Investigation of flap flexibility of β-secretase using molecular dynamic simulations

Flap motif and its dynamics were extensively reported in aspartate proteases, e.g. HIV proteases and plasmepsins. Herein, we report the first account of flap dynamics amongst different conformations of β-secretase using molecular dynamics simulation. Various parameters were proposed and a selected few were picked which could appropriately describe the flap motion. Three systems were studied, namely Free (BACEFree) and two ligand-bound conformations, which belonged to space groups P6122 (BACEBound1) and C2221 (BACEBound2), respectively and four parameters (distance between the flaps tip residue, Thr72 and Ser325, d1; dihedral angle, ϕ (Thr72-Asp32-Asp228-Ser325); TriCα angles, θ1 (Thr72-Asp32-Ser325), and θ2 (Thr72-Asp228-Ser325)) were proposed to understand the change in dynamics of flap domain and the extent of flap opening and closing. Analysis of, θ2, d1, θ1 and ϕ confirmed that the BACEFree adopted semi-open, open and closed conformations with slight twisting during flap opening. However, BACEBound1 (P6122) showed an adaptation to open conformation due to lack of hydrogen bond interaction between the ligand and flap tip residue. A slight flap twisting, ϕ (lateral twisting) was observed for BACEBound1 during flap opening which correlates with the opening of BACEFree. Contradictory to the BACEBound1, the BACEBound2 locked the flap in a closed conformation throughout the simulation due to formation of a stable hydrogen bond interaction between the flap tip residue and ligand. Analyses of all three systems highlight that d1, θ2 and ϕ can be precisely used to describe the extent of flap opening and closing concurrently with snapshots along the molecular dynamics trajectory across several conformations of β-secretase.


Introduction
Alzheimer's disease (AD) is a devastating neurodegenerative disorder that affects more than 30 million people worldwide, and is also the common cause of dementia in older people (Kandalepas et al., 2013;Rajapaksha, Eimer, Bozza, & Vassar, 2011). This type of disorder is characterised by the presence of insoluble amyloid plaques and neurofibrillary tangles in the brain (Cole & Vassar, 2007;Hampel & Shen, 2009;Kandalepas et al., 2013;Rajapaksha et al., 2011). The progression of the illness generally causes the functions of the brain to decline, initially affecting the cognitive abilities and memories and resulting in incapacitation and eventually death. Studies have shown that β-Secretase plays an important role in the development of AD (Gong et al., 2010;Luo & Yan, 2010;Vassar, 2004;Zhou et al., 2010). This has led to β-secretase being one of the prominent targets for pharmacological intervention. β-secretase, or BACE1 (β-site of Amyloid precursor protein Cleaving Enzyme), is an aspartyl protease found in humans (Barman, Schürer, & Prabhakar, 2011;Chakraborty, Kumar, & Basu, 2011;Gorfe & Caflisch, 2005;Shimizu et al., 2008;Xu et al., 2012). It is a membrane-anchored aspartic protease, with three distinct domains: (a) an N-terminal ectodomain, (b) a single transmembrane domain and (c) a cytosolic C-terminus. The active site of BACE1 is made up of two aspartate residues: Asp32 and Asp228 ( Figure 1). The R group of both aspartates coordinates a single water molecule between them, allowing for a nucleophilic attack to occur on the carbonyls. Over the active site, a hairpin loop termed 'flap' (residues 67-77) is present, which partially covers the cleft (this is a common feature in most of the aspartates). While the active site remains inactive, the flap stays in its open conformation (Chakraborty et al., 2011;Gorfe & Caflisch, 2005;Shimizu et al., 2008;Xu et al., 2012). However, the flap is stabilised while closed over its substrate or some other inhibitor. BACE1 also consists of a 10s loop (residues 9-14) ( Figure 1) that is located in the S3 pocket. Conformational changes in the loop are believed to control substrate binding in the active site, i.e. when the 10s loop takes an open conformation, it allows for greater binding between the substrate and the S3 pocket (Hong & Tang, 2004). The 10s loop also contains glycine residue (Gly11) with which the substrate can form a hydrogen bond, allowing for further stabilisation of the 10s loop, as well as the overall β-secretase-substrate interaction.
There are a considerable number of reported non-structural proteins that display unique and specific movements, these being vital for outlining their accurate functions (Ishima & Louis, 2008). Determining the correct parameters that can be used to precisely describe these motions of proteins is very important, in understanding the functions of an enzyme as well as its impact on drug binding and resistance. For instance, flap dynamics is a distinctive motion observed amongst aspartate proteases, e.g. HIV protease (Hornak, Okur, Rizzo, & Simmerling, 2006), cathepsin and plasmepsins (Karubiu, Bhakat, McGillewie, & Soliman, 2015). Flaps have been shown to regulate entry to the active site of proteases by providing access for substrate and/or inhibitor binding (Cai, Yilmaz, Myint, Ishima, & Schiffer, 2012;Ishima & Louis, 2008). Flap opening and closure in HIV protease (Hornak et al., 2006) have been well studied using molecular dynamics (Galiano, Bonora, & Fanucci, 2007;Ishima, Freedberg, Wang, Louis, & Torchia, 1999). In general, different parameters have been proposed to describe the flap motion in the case of the HIV protease (Galiano et al., 2007).
The distance between Ile50-Ile50′ (interflap distance) is one of the most commonly used parameters to define flap motions (Hornak et al., 2006). However, this parameter does not adequately describe the curling behaviour observed, which leads to the introduction of other parameters such as measuring dihedral angle to better define HIV protease flap dynamics (Bandyopadhyay & Meher, 2006). The flap dynamics in HIV protease is one such example, which highlights the importance of defining the appropriate parameters that best describe specific dynamics associated with the inhibitor binding to the binding cavity. In one of our recent studies, we proposed similar parameters in order to understand the extent of flap opening and closing in the case of plasmepsin II (Plm II) (Karubiu et al., 2015). The study highlighted that in the case of Plm II, the distance between flap tip residues Val78 and Leu292, and the angle between Val78, Asp214 and Leu292, will describe the opening and closing phenomena of flaps to a great extent.
Exploring the motion of the flap in order to understand mechanistic events associated with binding the BACE1 inhibitors is essential to design more potent inhibitors. Crystallographic studies revealed that the binding cavity of BACE1 is a flexible pocket, thus designing more potent inhibitors will require a better understanding of its flexibility in response to different inhibitors. Experimental parameters previously defined are inadequate to examine and predict the binding modes and flap dynamics of different inhibitors (Cole & Vassar, 2008;Klaver et al., 2010;Vassar, 2002). Furthermore, the exact parameters to precisely describe flap opening and closing are not well defined in the literature, neither experimentally nor computationally. This prompted us to report on the first detailed computational study that highlights the flap dynamics for BACE1 and various proposed parameters. We believe that this report serves as a benchmark for the flap flexibility phenomenon in BACE1 as well as potentially for other proteases.

System preparation
The Free crystal structure (PDB ID: 3TPL) (Xu et al., 2012) and two ligand bound conformations of BACE1 bound with a potent inhibitor BSIIV (PDB ID: 3TPP, 3TPR) (Xu et al., 2012) were retrieved from the RSCB Protein Data Bank. The missing residues were added using a graphical user interface of molecular modelling tool Chimera (Pettersen et al., 2004), and a ligand interaction map was generated using the web version of PoseView (Stierand & Rarey, 2010).

Molecular dynamic simulations & analysis
Prior to molecular dynamics simulation, the correct protonation states of all three systems were added using H++ server (Anandakrishnan, Aguilar, & Onufriev, 2012) and flipped residues were adjusted accordingly. The systems were then subjected to pre-dynamics preparation using typical parameters as described in our previous reports (Bhakat, Martin, & Soliman, 2014;Karubiu, Bhakat, & Soliman, 2014). A detailed Figure 1. Crystal structures of isolated BACE1, highlighting the active site aspartate residues: Asp32 and Asp228 (red spacefill representation), hairpin loop termed 'flap' (orange) and the 10s loop (green).
description of system preparation is also mentioned in the supplementary material. Finally, all three systems were subjected to an all-atom molecular dynamics simulation using GPU version of PMEMD engine provided with Amber 14 (Roe & Cheatham, 2013). The trajectories were saved at every 1 ps and analysed using CPPTRAJ module (Roe & Cheatham, 2013) integrated with Amber 14. The systems were visualised using the graphical user interface of UCSF Chimera (Pettersen et al., 2004), and plots were generated using the GUI of Microcal Origin data analysis software, version 6 (www. originlab.com).

Results and discussion
Recently, Xu et al. conducted a study where they established the conformations of Free (PDB: 3TPL, BACE Free ) and two BSIIV-BACE1 complexes belonging to space group P6 1 22 (PDB: 3TPR) and C222 1 (PDB: 3TPP), respectively. The space group P6 1 22 complex (BACE Bound1 ) was obtained by soaking the inhibitor (BSIIV) into a crystal of the free enzyme belonging to the same space group. However, for the C222 1 space group, crystals of the complex (BACE Bound2 ) were obtained by co-crystallizing the inhibitor (BSIIV) and the enzyme (Xu et al., 2012). Upon analysis (superimposition) of the free enzyme (BACE Free ) and the two complexes, it was observed that the flap (residues 67-75) had assumed a different conformation in each of the enzymes ( Figure 2). Furthermore, the inhibitor in the two complexes was bound at almost identical loci, while the residues around the ligand were found to have different conformations ( Figure 3). The conformation of Free and complexes belonging to space group P6 1 22 and C222 1 led to the adaptation of different flap conformations. The crystal structures of BACE1 complexes which belong to P6 1 22 space group remained in an open conformation due to lack of hydrogen bond interaction with the inhibitor and flap tip residues. The complex with C222 1 space group adapted a closed conformation (Xu et al., 2012), as direct hydrogen bonding interactions between main-chain atoms of flap residues and the bound inhibitor locked the flap in a closed conformation ( Figure 3).
Furthermore, the report provided by Xu et al. (2012) did not compare different parameters in order to understand the opening and closing phenomena of the flap dynamics. It is also worth mentioning that to date, no computational study has proposed the parameters needed to understand the extent of flap motion as well as the opening and closing phenomena in the case of BACE1. In this study, four parameters proposed to understand the opening and closing phenomena of the flap and its inhibitor uptake mechanism from a dynamic perspective. The choice of residues played an important role in defining the correct parameters to understand flap dynamics. Thr72 is often described as the flap tip residue, which, when combined with Ser325 (opposite loop of Thr72), helps to understand the extent of the flap movement. While the catalytic aspartyl residues Asp32 and Asp228 not only assist in understanding the position of Thr72 and Ser325 relative to catalytic dyad, they also provide an estimate of flap opening-closing during MD simulation ( Figure 4).
The movement of the protein from one configurational space to the next, is something much more complex than the variation of two, three or four single geometrical parameters. A proper sampling from one conformation to the other is conceptually the right way to study this kind of protein motions. In this way, a Potential Mean Force of these parameters will offer much more information than the variation along the MD. However, this study sets the path for future efforts of enhanced sampling in order to see the total conformational space travelled by the protein. Thus, a better future approach will be to use enhanced sampling techniques such as replica exchange and metadynamics to visit the full energy landscape of BACE systems in order to find all possible flap conformations and understand the complex biological motions involved in flap movement.

Molecular dynamic simulations and post-dynamic analysis of the flap motion
The positions of Asp228 and its companion dyad member Asp32 are quite fixed in all of the reported crystal structures. Therefore any change in the distance between them and the flap must be a consequence of either a conformational change in the flap or its movement. Therefore conventional studies have defined the open/close conformation by measuring the distances between the tip of the flap and the one of the asp dyad. For example, a closed conformation in which the flap moves closer to the two-Asp dyad has been observed in the crystal structures of complexes of BACE1 with inhibitors Hong et al., 2000;Hong, Turner, Koelsch, Ghosh, & Jordan, 2002;Turner, Hong, Koelsch, Ghosh, & Tang, 2005). In the apo structures, the flap was commonly observed to adopt an open conformation (in which the flap moves away to the two-Asp dyad) (Hong & Tang, 2004;Patel, Vuillard, Cleasby, Murray, & Yon, 2004;Shimizu et al., 2008;Xu et al., 2012). Therefore, it was imperative first to confirm that the results from this study were in agreement with previous findings Figure 3. Interaction maps of complexes belonging to space group P6 1 22 (A) and C222 1 (B), respectively. In the case of space group C222 1 , the inhibitor forms a critical hydrogen bond interaction with a flap residue of Gln73, which helps its locked (closed) conformation. Figure 4. Representation of the parameters used to understand the opening-closing phenomena of BACE1, distance (d 1 ) between the flap tip residues (Thr72-Ser325) and TriCα angles, θ 1 (Thr72-Asp32-Ser325), θ 2 (Thr72-Asp228-Ser325) and dihedral angle, ϕ (Thr72-Asp32-Asp228-Ser325). before we introduce the parameters. Therefore, we compared the experimental studies and the snapshots of the molecular dynamics trajectories (100 ns simulation) for all the systems. From the visual analysis of the snapshots, it was observed that BACE Bound1 transitions between a semi-open and close conformation throughout the duration of the simulation. Furthermore, inspection of the pattern of hydrogen bond interactions using LIG-PLOT showed that no direct hydrogen bonds formed between BSIIV and the flap in BACE Bound1 which could be the reason why the flap remains in the open conformation. In the case of BACE Bound2 , the flap tip residue Gln73 formed a key hydrogen bond interaction with the inhibitor, which led to a closed conformation of BACE Bound2 during the simulation period. The findings of this study are in agreement with the results reported by Xu et al. (2012), that in the crystal structure of complex (BACE Bound1 ) belonging to the P6122 space group the flap remains in the semi-open conformation, as the bound BSIIV does not form a stable hydrogen bond interaction with any of the flap tip residues. However, previous reports used only one parameter to define open/closed flap conformation (measuring the distances between the tip of the flap and the one of the asp dyad), but in dynamic state one parameter is insufficient to describe the complex motion of the flap. Thus, we embarked on a route where we have proposed additional parameters to measure flap dynamics, which we believe best describes not only flap and flexible region motions, but also the characteristic 'twisting' motion during the opening and closing of BACE1 enzyme ( Figure 5).
The per-residue C-α fluctuation (RMSF) of BACE Free , BACE Bound1 and BACE Bound2 highlighted the motion of its backbone during the simulation time ( Figure 6). BACE Free displayed a high overall conformational flexibility with an average RMSF of 0.95 Å, whereas the BACE Bound1 and BACE Bound2 displayed an average RMSF of 0.87 Å and 0.83 Å, respectively. A slight increase in average fluctuation (~0.04 Å) in the case of BACE Bound1 , as compared to BACE Bound2 , substantiates the fact that the lack of hydrogen bond formation between the inhibitor and the flap residues led to an increase in conformational flexibility.
The increase flap flexibility of BACE Bound1 was evidenced by the C-α fluctuation of flap tip residue, with Thr72 being 1.29 Å for BACE Bound1 , whereas in case of Figure 5. Transformation of flap position during simulation time for BACE Free (Black), BACE Bound1 (red) and BACE Bound2 (green). The starting and simulated structures are represented by cornflower blue and orange, respectively. BACE Bound2 , this fluctuation was significantly less (0.67 Å) due to the formation of a stable hydrogen bond interaction between the flap tip residue and the inhibitor. The BACE Free conformation showed an increase in flap flexibility, with a Thr72 C-α fluctuation of 1.58 Å, which is~0.30 Å and~0.90 Å greater than BACE Bound1 and BACE Bound2, respectively.
Furthermore, b-factor analysis (Figure 7) of BACE Free , BACE Bound1 (P6 1 22 space group) and BACE Bound2 (C222 1 space group) highlighted an increased flexibility at the hairpin region (containing flap tip residue, Thr72) and the flexible region (containing Ser325) in BACE Free as compared to BACE Bound1 . This was observed by the increased fluctuation of Thr72 and Ser325 during simulation that led to a fully open conformation in BACE Free . However, BACE Bound1 displayed significantly high fluctuations at the hairpin loop region/flap due to the lack of interaction between the flap tip residues and the inhibitor. However, it exhibited a higher flexibility at the 'flexible region' compared to that of BACE Bound2 . In the case of BACE Bound2 , the tight hydrogen bond interaction sealed the flap and flexible region in a stable conformation that led to lower the b-factor occupancy than for that of BACE Bound1 .
The trends in visual analysis of the snapshots, RMSF and b-factor fluctuations confirmed that BACE Free can adapt to an open conformation during simulation time in order to swallow the incoming inhibitor in its active site. The BACE Bound1 remained in an open conformation during simulation due to lack of hydrogen bond interaction and to steric hindrance. However, an opposite behaviour was observed in the case of BACE Bound2 , where the strong hydrogen bond interaction between the ligand and the flap residue induces a closed conformation thought simulation time.
3.1.1. Calculation of distance, d 1 , between flap tip residues and dihedral angle, phi (ϕ) The distance (d 1 ) between the flap tip residues and dihedral angle (ϕ), taking into account the flap tip residues and the catalytic dyads, emerged as the conventional strategy to understand the extent of the flap opening and twisting phenomena in HIV protease and plasmepsin (Karubiu et al., 2015). In the case of BACE1, the distance between C-α atoms of Thr72 and Ser325 was considered as d 1 . Furthermore, the dihedral angle (ϕ) amongst residues Thr72-Asp32-Asp228-Ser325 was taken into account in order to understand the flap twisting of BACE1. It is interesting to note that BACE Free displayed semi-open, open and closed conformations during the simulation, with d 1 values of 15.88 Å (~4.5 ns), 17.23 Å (~42 ns) and 5.36 Å (~76 ns), respectively (Figure 8(A) and 9(A)). This phenomena of BACE1, as described by d 1 , is in accordance with a similar phenomenon observed in HIV protease and plasmepsin (Karubiu et al., 2015). BACE Bound1 adapted a semi-open conformation during a very short period of simulation, with a d 1 of 15.61 Å at~42 ns (Figure 8(A) and 9(B)). The open flap conformation of BACE Bound1 could be due to the lack of hydrogen bond interaction between the inhibitor and the flap tip residue, and to steric hindrance as a result of the P6 1 22 space group. However, the value of d 1 remained stable in BACE Bound2 except small conformation in semi-open conformation, with an average value of 12.33 Å during simulation due to the critical hydrogen bond interaction between the flap tip residue and the inhibitor (Figure 8(A) and 9(C)).
Calculation of ϕ highlighted that the major twisting phenomena did not occur in all three systems. A slight flap twisting (lateral twisting) was observed in BACE Free (~28-42 ns) and BACE Bound1 (~20-42 ns) (Figure 8(B)). However, only some degree of drift in ϕ was observed in both cases, which confirmed that parameterising flap twisting played a lesser role in opening and closing the flaps. Similar to d 1 , the values of ϕ remained constant in the case of BACE Bound2 , with an average value of −42.24° (Figure 8(B)).
3.1.2. Calculation of TriCα Angles, θ 1 and θ 2 Calculating the fluctuation of TriCα angles, θ 1 and θ 2 has emerged as a useful strategy to understand the extent of the flap opening in the case of plasmepsin (Karubiu et al., 2015). Calculating θ 1 (Thr72-Asp32-Ser325) and θ 2 (Thr72-Asp228-Ser325) proved to be a useful parameter in order to understand the opening and closing phenomena of BACE1. Calculating θ 1 led to identifying the flap opening and closing confirmations. In the case of BACE Free , the semi-open and open states correspond to a value of 60.82°(~4.5 ns) and 57.76°(~42 ns), whereas Figure 6. The C-α RMSF of BACE Free (black), BACE Bound1 (red) and BACE Bound2 (green), respectively. The flap tip region with increased flexibility is highlighted in the crystal structure. the closing has a value of 18.09°(~77 ns) (Figure 8(C)). The trend of θ 1 is similar to that of d 1 , but the extent of the semi-open and open states is not described properly due the presence of Asp32 in close proximity to the flexible 10s loop. For BACE Bound1 , the flap opening is described by a θ 1 of 52.94°(~42 ns) (Figure 8(C)), which is in agreement with the fluctuation of d 1 . The fluctuation of θ 1 , in the case of BACE Bound2 , remained at a steady state, with an average value of 41.09° (  Figure 8(C)), this possibly being due to the formation of a stable hydrogen bond interaction between the flap tip residue and the inhibitor. The θ 2 values for BACE Free were~73.88°(~4.4 ns), 78.32°(~42 ns) and~29.61°(~7 8 ns), which corresponded to the semi-open, open and closed conformations, respectively. The fluctuation of the θ 2 values showed a trend similar to that of d 1 (Figure 8(D)). The opening of BACE Bound1 was described by a θ 2 value of 64.25°at~42 ns, while the BACE Bound2 displayed a steady average fluctuation of 56.82°with a trend similar to that of d 1 (Figure 8(D)).
The proper representation of three systems using the angle created amongst residues Thr72-Asp228-Ser325 (θ 2 ) highlighted the fact that θ 2 can be used in combination with d 1 to understand the extent of the flap opening and closing phenomena in BACE1.

Radius of gyration and principal componenet analysis (PCA)
The R g is defined as the moment of inertia of the C-α atoms from its centre of mass. The tendency of semiopen, open and closed phenomena led to a similar trend of R g (Figure 8(D)) to that of d 1 and θ 2 , thus confirming the trend in flap dynamics in the case of BACE Free . The lack of interaction between the inhibitor and the flap in BACE Bound1 led to an instable fluctuation of R g in the cause of simulation time. The stable hydrogen bond interaction in c BACE Bound2 led to a stable fluctuation of R g . The trend of R g in all three systems was similar to that of d 1 and θ 2 , and thus highlights the flap dynamics   Figure 9. It is interesting to note that this plot helps in identification of the region of these conformations during simulation time as often a particular value of certain parameters don't describe the range of certain ensembles. Figure 10 highlights that BACE Free opened with a later twist in the flap residues (which includes Thr72) and flexible region (Ser325) during the simulation period, which can be highlighted by parameters d 1 and θ 2 , where ϕ described the extent of lateral shift to some extent. The BACE Bound1 adapted an open conformation during a short period of the simulation time due to a lack of ligand-protein hydrogen bond interaction at the flap region. However, the presence of hydrogen bond interaction between the flap residue and the ligand locked the BACE Bound2 in a closed conformation during simulation, which is in accordance with the steady fluctuation of d 1 , θ 2 and ϕ. Thus, combining the distance of d 1 , the tri C-α angle, θ 2 and the dihedral angle, ϕ, will provide a better representation of the flap opening and closing, with a slight lateral twisting in case BACE Free and BACE Bound1 .

Conclusion
The binding cavity of BACE1 can adapt according to the size and shape of different inhibitors as well as the space group. The presence of P6 1 22 and C222 1 space groups changed the conformational landscape of BACE1. Hence, defining the appropriate parameters of the flap motion will assist in understanding the binding landscape of novel inhibitors in the active site cavity, as well as establish its effect on the flap motion. In some cases, the Arg235 residue can rotate to increase the binding space in order to incorporate a larger inhibitor. In addition, the inhibitor that formed a hydrogen bond interaction with the flap residues locked its motion. Understanding flap dynamics will therefore play a key role in the ligand uptake mechanism.
This study took into account the Free (BACE Free ) and two bound conformation (BACE Bound1 , BACE Bound2 ) systems. In BACE Free , a symmetrical flap opening with lateral twisting was observed, while BACE Bound1 described an asymmetrical flap opening (opening of hairpin loop), which might be due to the lack of hydrogen bond connection between inhibitor and the residue of the hairpin loop. Overall, d 1 , θ 2 and ϕ described the motion of the flap in all three systems, where d 1 and θ 2 accurately indicated the opening and closing phenomena of flap residues, and ϕ the later twisting that was observed in case of BACE Free . Thus, combining d 1 , θ 2 and ϕ and snapshots along the MD trajectory will help to understand the flap motion across different conformations of BACE1. The methods used in this study will not only help to parametrise the flap dynamics in the case of BACE-2, but also across different aspartyl proteases, thereby improving our understanding of the effect of the flap motion in ligand uptake and binding.