Preprint
Article

Identification Mechanism of BACE1 on Inhibitors Probed by Using Multiple Separate Molecular Dynamics Simulations and Comparative Calculations of Binding Free Energies

Altmetrics

Downloads

107

Views

39

Comments

0

A peer-reviewed article of this preprint also exists.

Submitted:

31 May 2023

Posted:

01 June 2023

You are already at the latest version

Alerts
Abstract
The β-amyloid cleaving enzyme 1 (BACE1) is regarded as an important target of drug design toward treatment of Alzheimer's disease (AD). In this study, three separate molecular dynamics (MD) simulations and calculations of binding free energies were carried out to comparatively probe identification mechanism of BACE1 on three inhibitors 60W, 954 and 60X. The analyses on MD trajectories indicate that the presence of three inhibitors influences structural stability, flexibility and internal dynamics of BACE1. Binding free energies calculated by using solvated interaction energy (SIE) and molecular mechanics generalized Born surface area (MM-GBSA) methods reveal that the hydrophobic interactions provide decisive forces for inhibitor-BACE1 binding. The calculations of residue-based free energy decomposition suggest that the sidechains of residues L91, D93, S96, V130, Q134, W137, F169 and I179 play key roles in inhibitor-BACE1 binding, which provides a direction for future drug design toward treatment of AD.
Keywords: 
Subject: Biology and Life Sciences  -   Biophysics

1. Introduction

Alzheimer's disease (AD) is the most common chronic neurodegenerative disease, characterized by early clinical manifestations of memory impairment [1,2]. As the disease progresses, patients endure various neuro-psychiatric symptoms, including language disorders, emotional disturbances, motor impairments, visuospatial skill impairments, and even personality and behavioral changes, which seriously affect social, occupational, and daily life, leading to gradual loss of body functions and ultimate death [3]. The previous reports indicated that the development and deterioration of the AD are closely related to two major pathological features, involving amyloid plaques containing amyloid-β (Aβ) peptide and neurofibrillary tangles (NFTs) rich in tau [4,5,6]. The formation of Aβ peptide with a 38−43 amino acid length and subsequent aggregation is the culprit responsible for AD, as it blocks the transmission between neurons, ultimately leading to neuronal death [7,8]. The toxic amyloid-β peptides are yielded through procession of the β-amyloid precursor protein (APP) with β-amyloid cleaving enzyme 1 (BACE1) and γ-secretase [9,10]. Therefore, effective inhibition of the activity of BACE1 is a useful approach toward the treatment of AD.
In recent years, the pathway for inhibiting the activity of BACE1 through the design of small molecules has been paid widespread attentions [11,12,13,14,15,16]. Several potential BACE1 inhibitors have been designed by different groups and tested in clinical trials [17,18,19]. Unfortunately, the efficiency of inhibitors is greatly limited due to their side effects. To relieve this issue, many works, involved in insights into the factors affecting inhibitor-BACE1 binding, have been performed to probe efficient schemes for design of potent inhibitors [17,20,21,22,23,24]. Currently, it is still a great challenge to design the potent BACE1 inhibitors toward treatment of AD. Although some BACE1 inhibitors are proposed, molecular mechanism with regard to inhibition of the BACE1 activity is not insufficient, which exerts a heavy limit on development of clinically available BACE1 inhibitors. Therefore, it is highly necessary to explore binding mechanisms of inhibitors to BACE1 at atomic levels for designing efficient BACE1 inhibitors toward treatment of the AD.
Multiple computational technologies, such as molecular dynamics (MD) simulations [25,26,27,28,29,30], calculations of binding free energies [31,32,33,34] and principal component analysis (PCA) [35,36,37,38] etc., play significant roles in probing atomic-level binding mechanism of inhibitors to targets. Conventional MD simulations (cMD) are usually used to obtain conformational samplings of inhibitor-target complexes, but multiple separate MD (MSMD) simulations recently adopted by different work groups can rationally improve sampling efficiency of conformations [39,40,41,42,43,44,45]. Binding free energy calculations are applied to evaluate binding ability and modes of inhibitors to targets, which are involved in molecular mechanics Poisson Boltzmann/generalized Born surface area (MM-PB/GBSA) [46,47,48], solvated interaction energy (SIE) [49], thermodynamics integration (TI) [50,51,52] and free energy perturbation (FEP) [53,54,55,56]. Although FEP and TI methods can provide more accurate results, they are extremely time-consuming and require sufficient statistical samplings. Compared to the FEP and TI methods, MM-PB/GBSA and SIE methods can obtain fast and rational results in predictions of binding free energies. Interestingly, MD simulations and binding free energy calculations have been utilized to investigate inhibitor-BACE1 binding mechanism [16,57,58,59,60]. For instance, Chen et al. applied Gaussian accelerated molecular dynamics simulations and MM-GBSA calculations to study effect of pH-dependent protonation on inhibitor-BACE1 binding and their results revealed that pH-dependent protonation strongly affected the structural flexibility and correlated motions of BACE1 [61]. Hatmal and the coworkers combined MD simulations and ligand-receptor contacts analysis to develop valid pharmacophore models and their work rationally guided pharmacophore design [62]. Despite these successful studies, it is still of high significance to deeply probe molecular mechanism of inhibitor-BACE1 binding for design of clinically available drugs for treatment of the AD.
In this work, to rationally probe atomic-level binding mechanism of inhibitors to BACE1, three inhibitors 60W, 954 and 60X [63], indicated by using their identity document (ID) in the protein data bank (PDB), were selected to realize our aims. The topology structure of inhibitor-BACE1 complex and binding pocket were respectively depicted in Figure 1A,B. The structures of 60W, 954 and 60X were separately displayed at Figure 1C, 1D and 1E. As shown in Figure 1C–E, three inhibitors share similar molecular scaffold and have small structural difference. Binding ability of 60W, 954 and 60X to BACE1 is scaled by the Ki values of 1, 45 and 48 nM. Insights into effect of tiny structural difference on conformational changes of BACE1 will be of importance for design of potent inhibitors. To achieve our goal, MSMD simulations, MM-GBSA, SIE, PCA , dynamics cross-correlation maps (DCCMs) and free energy landscapes (FELs) were combined together to perform this current study and we expected three following tasks can be reached: (1) the changes in conformations and free energy profiles can be revealed through PCA and constructions of FELs, (2) binding free energies were estimated by using the MM-GBSA and SIE methods to evaluate the inhibitor-BACE1 binding ability and (3) hot interaction hotspots of inhibitors with BACE1 were identified through calculations of residue-based free energy decomposition. This work is also expected to provide useful information and theoretical guidance for design of efficient inhibitors toward BACE1.

2. Results and Discussion

2.1. Dynamics Equilibrium and Structural Fluctuation

To check structural stability of BACE1 during three separate MD simulations, root-mean-square deviations (RMSDs) of backbone atoms in BACE1 relative to the initially optimized structures were calculated based on three separate MD trajectories and the results were provided in supporting information Figure S1. It is noted that all four systems have reached the equilibrium after 300 ns of three separate MD simulations. The greater frequency distributions of the apo, 60W-, 954- and 60X-bound BACE1 are respectively located at 2.11, 2.29, 1.93 and 2.06 Å, moreover the distribution shape of the 954- and 60X-bound BACE1 moves toward the left relative to that of the apo BACE1 while the distribution shape of the 60W-bound BACE1 moves toward the right (Figure 2A). Thus, binding of 954 and 60X weakens the structural fluctuation of BACE1 compared to the apo state while the presence of 60W increases the structural fluctuation of BACE1. The RMSDs of heavy atoms from three inhibitors relative to the initially optimized structures were also computed to check their stability in binding pocket of BACE1 and the frequency distributions were depicted at Figure 2B. The RMSD of 60W is distributed at two peak values of 0.68 and 1.77 Å while the RMSDs of 954 and 60X are respectively populated at the peak values of 3.05 and 1.04 Å (Figure 2B). These results indicate that 60W and 60X have lower structure fluctuations in binding pocket of BACE1 than 954, implying that 60W and 60X have higher structure stability in binding pocket of BACE1 than 954.
To examine inhibitor’s binding-mediated impacts on structural flexibility of BACE1, root-mean-fluctuations (RMSFs) were computed by using the coordinates of the atoms Cα kept at the SMT (Figure 2C). It is found that the apo and bound states of BACE1 share similar flexible and rigid regions. The RMSF difference between the apo and bound states was also calculated by utilizing the equation R M S F = R M S F b o u n d R M S F a p o , in which R M S F , R M S F b o u n d and R M S F a p o represent the RMSF difference, the RMSFs of the bound and apo states (Figure 2D). The presence of the three inhibitors weakens the structural flexibility of the regions D1 (residues 110-123) and D4 (residues 178-202), which makes these two regions be more rigid than the apo BACE1 (Figure S2). However, binding of the three inhibitors strengthens the structural flexibility of the regions D2 (residues 146-153) and D3 (residues 164-178) relative to the apo BACE1 (Figure 2D and Figure S2). The presence of 60W and 60X totally enhances the structural flexibility of the region D5 (residues 217-245) but binding of 954 weakens that of this region relative to the apo BACE1 (Figure 2D). By comparison with the apo BACE1, binding of 60W increase the structural flexibility of the region D6 (residues 363-387) while the presence of 954 and 60X slightly weakens that of this region (Figure 2D and Figure S2)
To access stability of secondary structures for BACE1, the combination of the program CPPTRAJ and the DSSP second structure analysis [64] were used to probe the stability of the second structure of BACE1 from the apo and bound states in three separate MD simulations. The time evolutions of the secondary structures for the apo, 60W-, 954- and 60X-bound BACE1 were displayed in Figure 3A–C, S3, S4 and S5, individually. It is observed that the secondary structure of the apo BACE1 is stable through three separate MD simulations (Figure 3A,B). By comparison with the apo BACE1, the secondary structures of the 60W-, 954- and 60X-bound BACE1 do not generate obvious changes (Figures S3–S5), indicating that binding of inhibitors hardly affects the stability of the secondary structures of BACE1. To understand influences of inhibitor binding on the structure compact extents of BACE1, the gyrations of BACE1 in the apo and bound states were estimated based on the SMT and their frequency distributions were provided at Figure 3D. The gyration of the apo BACE1 is distributed at two peaks of 20.87 and 21.02 while the ones of the 60W-, 954- and 60X-associated BACE1 are separately populated at the single peaks of 21.17, 20.99 and 20.99 Å, furthermore the distribution shapes of gyrations for BACE1 in three bound states totally move toward the right by referencing to the apo BACE1 (Figure 3D), implying that the presence of the three inhibitors slightly leads to a more incompact structures of BACE1.
Based on the current analyses, binding of the three inhibitors affects structural fluctuations and compact extents of BACE1 but hardly changes the stability of secondary structures for BACE1. The difference in structures of 60W, 954 and 60X also leads to their different stability in binding pocket of BACE1. Meanwhile, binding of inhibitors generates obvious effect on structural flexibility of BACE1, which also implies possible hot spots of inhibitor-BACE1 interactions. These current results basically agree with that of a previous work [16].

2.2. Conformational Changes of BACE1 and Free Energy Profiles

To clarify inhibitor-mediated changes in correlated motions of BACE1, DCCMs were computed by means of the coordinates of the atoms Cα in BACE1 and the results were exhibited at Figure 4. The color bar coded by different colors was employed to embody the contents of the correlated motions between residues of BACE1. For the apo BACE1, several strongly correlated motions are observed: (1) the region R1 describes strong PCMs of the N-terminal from BACE1 relative to itself (Figure 4A), (2) the region R2 reflects strong ACMs of residues 206-265 relative to residues 115-158 and the region R3 characterizes strong ACMs between residues 278-331 and 146-200, (3) the region R4 embodies strong ACMs of residues 347-418 relative to the N-terminal of BACE1 and (4) the region R5 describes strong PCMs of residues 266-328 relative to themselves (Figure 4A). Compared to the apo BACE1, binding of the three inhibitors highly weakens the PCMs occurring at the regions R1 and R5 (Figure 4B–D). Meanwhile, binding of 60W, 954 and 60X also reduces the ACMs of the regions R2 and R3 by referencing to the apo BACE1 (Figure 4B–D). Differently, binding of 60W and 954 obviously abates the ACMs of the region R4 by comparison with the apo BACE1 (Figure 4B,C) but the presence of 60X slightly strengthens the ACMs of this region (Figure 4D). The aforementioned regions corresponding to evident alterations of correlated motions are possibly involved in hot interaction spots of inhibitors with BACE1.
To reveal impacts of inhibitor binding on concerted movements of structural domains in BACE1, PCA was carried out by using the CPPTRAJ in Amber 20. The first eigenvector from the PCA was visualized by means of the software VMD [65] and the results were depicted Figure 5. It is observed that inhibitor binding yields evident influences on the collective motions of the α helix α1 and the loop L2 (Figure 5). In the apo BACE1, the α1 and L2 generate a parallel concerted motion with the same direction (Figure 5A). Compared to the apo BACE1, binding of 60W enhances the concerted motion of α1 and L2 and obviously alters the direction of the concerted motion for the L2 (Figure 5B). By comparison with the apo BACE1, the presence of 954 not only leads to a completely opposite motion direction of the α1 and L2 but also inhibits the motion amplitude of the L2, which produces a tendency of the L2 away from the α1 (Figure 5C). By referencing to the apo BACE1, binding of 60X not only changes the concerted motion direction of the α1 and L2 but also apparently weakens the concerted motion of the L2 (Figure 5D). Binding of 60W and 954 slightly inhibits the concerted movement of the loop L1 relative to the apo BACE1 but the presence of 60X slightly strengthens the concerted motion of this loop (Figure 5B–D). In addition, binding of the three inhibitors weakens concerted motion of the loop L3 compared to the apo BACE1 (Figure 5B–D).
To unveil free energy profiles of the BACE1 conformation changes caused by inhibitor binding, FELs were created by using the projections (PC1 and PC2) of the SMT on the first two eigenvectors as reaction coordinates (RCs) and the presentative structures relating with free energy profiles were depicted in Figure 6 and S6. The projections of MD trajectories can rationally reflects conformational changes of BACE1, which is the cause for our selecting them as RCs.
In the case of the apo BACE1, three separate MD simulations capture three free energy valleys (EVs), including EV1, EV2 and EV3 (Figure S6A). According to the color bar, three EVs are located at the valley bottoms of the same depth (Figure S6A). Three presentative structures of the apo BACE1 in the EV1-EV3 were superimposed together (Figure S6B). The results suggest that the domains D4-D6 yield obvious deviations from each other. The structure domains D1 and D3 generate slight deviations and the D2 produces the slight sliding (Figure 6B)
Compared to the apo BACE1, binding of 60W and 954 only results in two EVs (Figure 6A,D), which is less than the energy states of the apo BACE1. This result implies that binding of 60W and 954 induces conformational arrangement of BACE1 relative to the apo BACE1. The representative structures of the 60W- and 954-bound BACE1 situated at the EV1 and EV2 were superimposed together to probe structural difference (Figure 6B,E). By comparison with the apo BACE1, binding of 60W and 954 reduces the structural deviations of the structural domains D1, D3, D4 and D5 (Figure 6B,E and S6B). Although binding of 60W weakens the structural deviation of the D2 and D6 relative to the apo BACE1, the association of 954 leads to great deviation of the D2 and D6 (Figure 6B,E and S6B). As shown in structural alignment of inhibitors 60W and 954 falling into the EV1 and EV2 (Figure 6C,F), 60W and 954 produce slightly sliding between two energy states, which possibly affects binding of these two inhibitors to BACE1. As for the 60X-bound BACE1, three EVs are detected through the entire MD simulations (Figure 6G). Although binding of 60X does not alter the number of the EVs relative to the apo BACE1, the presence of 60X enhances the energy barrier between the EV3 and two states EV1 and EV2 by referencing to the apo BACE1 (Figure 6G and S6A), which correspondingly increases the difficulty of the transitions between the EV3 and the EV1 and EV2. According to structural superimposition of the 60X-bound BACE1 trapped at the EV1-EV3 (Figure 6H), except for the structural domains D1, D3, D4 and D5, binding of 60X evidently increases the deviation of the D2 and D6 among three EVs compared to the apo BACE1. The structural alignment of 60X falling into the EV1-EV3 indicates that 60X yields slight deviations among three energetic states (Figure 6I).
Based on the aforementioned calculations of DCCMs, PCA and analyses of FELs, binding of inhibitors changes correlated motions between residues, affects concerted movements of structural domains and alters free energy profiles of BACE1. Some of the structural domains affected by inhibitor binding are located near the binding pocket, hence conformational changes caused by binding of inhibitors in turn alter the activity of BACE1. Several previous works also revealed similar results [20,61], which is in basic agreement with our current work.

3.3. Comparative Calculations of Binding Free Energies

To access binding ability of 60W, 954 and 60X to BACE1, the SIE method was applied to calculate binding affinities to the three inhibitors to BACE1 by using 500 snapshots extracted from the equilibrated section of three separated MD simulations, namely for the SMT, in a time interval of 1.8 ns. The calculated results were listed in Table 1. It is observed that the rank of binding affinities predicted by the SIE method is in consistence with that of the experimental values, which indicates that our current free energy analyses are rational and reliable.
According to Table 2, the components of binding affinities predicted by the SIE method mainly consist of the intermolecular Coulomb interactions ( E C ), the van der Waals ones ( E v d W ), the reaction energy ( G R ) and the energy changes in the molecular surface area upon binding( γ × M S A ). The energy contributions favoring binding of inhibitors are those from the van der Waals interactions between binding partners (-42.50 to -54.53 kcal·mol-1), the intermolecular Coulomb interactions (-12.55 to -16.18 kcal·mol-1) and the energy contributions relating with the changes in the molecular surface (-8.05 to -10.63 kcal·mol-1). The reaction energies fluctuate at a range of 21.27 to 24.77 kcal·mol-1 and this component provide an unfavorable force for the inhibitor bindings, which is also revealed by the previous work [49,66,67]. On the basis of Table 2, the unfavorable reaction energies of three inhibitor-BACE1 complexes are partially compensated by the favorable intermolecular Coulomb interaction. Meanwhile, the intermolecular van der Waals interactions also contribute partial compensation to this unfavorable effect. Among three inhibitors, 60W shows the strongest binding ability to BACE1 (-8.82 kcal·mol-1) while 954 has the weakest binding ability to BACE1 (-7.30 kcal·mol-1), which suggests that small structure difference among the three inhibitors impacts their binding ability to BACE1.
To comparatively study binding strength of 60W, 954 and 60X to BACE1, MM-GBSA method was adopted to predict binding free energies of three inhibitor-BACE1 complexes based on 500 snapshots extracted from the equilibrated section of three separated MD simulations, namely for the SMT, in a time interval of 1.8 ns. Because of expensive time in the entropy calculation, 100 snapshots taken from the above mentioned 500 snapshots were employed to perform the calculation of the entropy contributions to inhibitor-BACE1 binding. The MM-GBSA calculations are possibly involved in multiple generalized Born (GB) models. To understand influences of different GB models on the predicted results, four GB models, indicated by IGB=1, IGB=2, IGB=5 and IGB=66, were chosen to estimate binding free energies of the three inhibitors to BACE1. The empirical parameters involved in calculations of four GB models were given in Table 2, which includes two empirical parametersγand β together with the radii types. Binding free energies and their components computed by the MM-GBSA method were listed in Table 3.
Binding free energies are mainly composed of five components, including the van der Waals interactions ( E v d W ), the electrostatic interactions ( E e l e ), the polar solvation free energy ( G g b ), the non-polar solvation free energy ( G s u r f ) and the entropy contributions ( T S ), which is shown in Table 3. From free energy components, E v d W , E e l e and G s u r f are favorable for inhibitor-BACE1 binding but G g b and T S impair inhibitor-BACE1 associations (Table 3). The hydrophobic interactions ( G h y d r o ) formed by the sum of E v d W and G s u r f are favorable for inhibitor-BACE1 binding. The polar interactions ( G p o l ) formed by the sum of E e l e and G g b provide different-type force for inhibitor-BACE1 association. In details, G p o l predicted by the models of IGB=1, IGB=2 and IGB=66 is unfavorable for the inhibitor-BACE1 binding while that predicted by the model of IGB=5 contributes favorable forces to the inhibitor-BACE1 associations (Table 3). The sum of four components E v d W , E e l e , and G g b constructs the enthalpy contributions (∆H) to the inhibitor-BACE1 binding. Based on Table 3, the GB models used for calculations of MM-GBSA only produce evident impacts on polar solvation free energies and the selection of the empirical parameters γ and β obviously affects calculations of non-polar solvation free energies. By comparison on four GB models, the GB model of IGB=5 leads to the weakest polar solvation free energies for all inhibitors but that of IGB=66 yields the strongest polar solvation free energies (Table 3). Correspondingly, the GB model of IGB=5 generates the strongest enthalpy contributions to inhibitor-BACE1 association but that of IGB=66 produces the weakest enthalpy contributions to inhibitor-BACE1 binding. As a result, the selection of the GB models brings a vital impact on predictions of inhibitor-BACE1 binding free energies.
For our used GB models, binding free energies of 60W, 954 and 60X to BACE1 estimated with the GB model of IGB=2 is mostly close to the experimental values. Differently, binding free energies of the three inhibitors to BACE1 calculated through the GB models of IGB=1, 5 and 66 highly deviate from the experimental results. Meanwhile, the rank for binding free energies of 60W, 954 and 60X in four GB models are also in good agreement with that for the experimental values, verifying that our current results are reliable and rational. Based on the aforementioned analyses, the results calculated by the GB model of IGB=2 are utilized to access binding difference of the three inhibitors to BACE1. The electrostatic interactions of 60W and 60X with BACE1 are respectively strengthened by 7.98 and 5.63 kcal/mol relative to that of 954 with BACE1 but unfavorable polar solvation free energies of the 60W- and 60X-BACE1 complexes are raised by 11.42 and 3.99 kcal/mol by comparison with that of the 954-bace1 complex. On the whole, the polar interaction of 60W with BACE1 is increased by 3.44 kcal/mol by referencing to that of 954 with BACE1 while the polar interaction of 60X with BACE1 is reduced by 1.64 kcal/mol. The hydrophobic interaction of 60W with BACE1 strengthened by 13.3 kcal/mol compared to that of 954 with BACE1, but the hydrophobic interaction of 60X with BACE1 hardly changes relative to that of 954 with BACE1. As a result, the enthalpy contributions to the 60W- and 60X-BACE1 binding is improved by 9.86 and 1.63 kcal/mol by referencing to that of the 954-BACE1 binding. In addition, the unfavorable entropy contributions to the 60W- and 60X-BACE1 binding is increased by 3.92 and 0.34 kcal/mol relative to the 954-BACE1 binding. In summary, the binding ability of 60W and 60X to BACE1 is strengthened by 5.94 and 1.29 kcal/mol compared to that of 954 to BACE1 (Table 3). Therefore, although structural difference of the three inhibitors is tiny, their binding ability to BACE1 produces the bigger difference according to our current calculations, which should owe to the conformational changes caused by their binding.
By combination of the SIE and MM-GBSA calculations, it is found that hydrophobic interactions provide a key contribution to inhibitor-BACE1 binding, which is in good agreement with the previous reports [16]. Thus, the rational optimization on inhibitor-BACE1 hydrophobic interactions is of high significance for successful design of clinically available inhibitors toward BACE1. Based on this issue, more attentions should be paid to hydrophobic interactions of inhibitors with BACE1.

2.4. Analyses of Inhibitor-BACE1 Interaction Networks

To obtain atomic-level insights into interaction modes of inhibitors with BACE1, residue-based free energy decomposition method was applied to estimate inhibitor-residue interaction spectrum of three inhibitor-BACE1 complexes (Figure 7). The contributions from the sidechains and backbones of residues to inhibitor-BACE1 associations were provided in Table 4. The hydrogen bonding interactions (HBIs) between inhibitors and residues of BACE1 were analyzed by using the program CPPTRAJ and the results were listed in Table 5. The geometric information with regard to inhibitor-residue interactions was depicted in Figure 8. Meanwhile, the frequency distributions of the distances relating with inhibitor-residue interactions were also calculated and the results were displayed in Figure 9.
For the 60W-BACE1 complex, 60W produces the interactions stronger than -1.0 kca/mol with six residues of BACE1, including L91, D93, S96, V130, Q134 and I179 (Figure 7A and 7D). The three residues D93, S96 and V130 are situated near the hydrophobic rings R1 and R2 of 60W (Figure 8A). Hence, D93 forms the CH-O interactions with these two rings, S96 generates the CH-π and CH-O interactions with the ring R1 and V130 yields the CH-π interaction with the ring R1 of 60W (Figure 8A). According to Table 4, the energetic contributions of S96 and V130 to the 60W-BACE1 binding mostly arise from the sidechains of these two residues. Additionally, the carbonyl of D93 generates four HBIs with the ring R2 of 60W and their occupancy is higher than 46.7% (Table 5 and Figure 8B), meanwhile the favorable 60W-D93 interaction mainly comes from the electrostatic interaction of the sidechain of D93 (Table 4). On the whole, D93, S96 and V130 provide energy contributions of -2.04, -1.13 and -1.75 kca/mol to the 60W-BACE1 binding, respectively (Figure 7A,D and Table 4). The distances for the mass centers of the sidechains of V130 and S96 away from that of the ring R1 are respectively distributed at 4.03 and 4.03 Å (Figure 9A), which verifies the interactions of these two residues with 60W. The distance between the mass center of the carbonyl of D93 and that of the ring R2 in 60W is populated at 6.09 Å, which agrees with the weak CH-O interaction of D93 with 60W (Figure 9A). The residues Q134 and I179 are next to the ring R3 of 60W and these two residues form the CH-π interactions with the ring R3 of 60W (Figure 7A,D and Figure 8A). As shown in Table 4, the van der Waals interactions of the sidechains from Q134 and I179 with the ring R3 of 60W contributes the most forces to the 60W-BACE1 association. The distance of the carbon atom from Q134 and the mass center of the alkyl group in I179 away from mass center of the ring R3 in 60W are situated at 4.03 and 3.66 Å, respectively (Figure 9A). As a result, the two residues Q134 and I179 separately provide the interaction energies of -2.58 and -2.26 kcal/mol for binding of 60W to BACE1 (Figure 7A,D). The interaction energy of L91 with 60W is -1.72 kcal/mol (Figure 7A,D), which structurally stems from the CH-π interaction between the alkyl group of L91 and the ring R4 of 60W (Figure 8A). More interestingly, the energy contribution of L91 is mainly provided by the van der Waals interactions between the sidechain of L91 and the ring R4 of 60W (Table 4). The distance between the mass center of the alkyl group from L91 and that of the ring R4 in 60W is distributed at 3.84 Å (Figure 9A), which demonstrates the existence of the CH-π interaction between 60W and L91.
With respect to the 954-BACE1 compound, five residues are involved in interactions stronger than -1.0 kcal/mol with the inhibitor 954 and these residues include L91, Y132, W137, F169 and I179 (Figure 7B,D). The interaction energies of Y132 and W137 with 954 are -1.88 and -1.04 kcal/mol, individually, which structurally agrees with the π-π interactions of the phenyl group in Y132 with the ring R2 of 954 and the hydrophobic ring of W137 with the ring R1 of 954 (Figure 8C). The distances of the mass centers for the hydrophobic rings of Y132 and W137 away from that of the rings R2 and R1 from 954 are located at 4.87 and 5.09 Å (Figure 9B), which further supports the interactions of these two residues with 954. Based on Table 4, the energy contributions of Y132 and W137 to the 954-BACE1 binding are mostly provided by the van der Waals interactions of the sidechains from Y132 and W137 with 954. The hydrophobic groups L91, F169 and I179 are located near the ring R4 of 954 (Figure 8C). Therefore the alkyl group of L91, the phenyl group of F169 and the alkyl group of I179 tend to generate the CH-π, π-π and CH-π interactions with the R4 of 954. The distances of the mass centers of the sidechains in L91, F169 and I179 away from that of the ring R4 in 954 are respectively populated at 4.86, 6.42 and 4.24 Å (Figure 9B), which verifies the hydrophobic interactions of these three residues with 954. On the whole, L91, F169 and I179 contribute the interaction energies of -1.14, -1.44 and -1.91 kcal/mol to the 954-BACE1 binding (Figure 7B,D and Table 4). More importantly, the interaction energies of L91, F169 and I179 with 954 mostly origin from the van der Waals interactions of the sidechains in these three residues with 954 (Table 4). In addition, the carbonyl group of D93 forms four HBIs with the ring R3 of 954 and the occupancy of these four hydrogen bonds are higher than 46.9% (Table 5 and Figure 8D). However, D93 only provides an energy contribution of -0.77 kcal/mol (Figure 7B,D), which mainly stems from the electrostatic interaction between the sidechain of D93 and 954 (Table 4).
With regard to the 60X-BACE1 complex, 60X yields the interactions stronger than -1.0 kcal/mol with five residues D93, V130, Y132, W137 and I179 in BACE1 (Figure 7C,D). The hydrophobic groups of I179 and Y132 are adjacent to the ring R2 of 60X (Figure 8E), hence the alkyl group of I179 forms the CH-π interactions with the ring R2 of 60X and the phenyl group of Y132 generates the π-π interaction with the ring R2 of 60X (Figure 8E). The distances of the mass centers for the hydrophobic groups of Y132 and I179 away from that of the ring R2 in 60X are individually populated at 4.52 and 5.13 Å (Figure 9C), which further supports the interactions of Y132 and I179 with 60X. Y132 and I179 separately contribute the interaction energies of -2.0 and -2.31 kcal/mol to the 60X-BACE1 binding (Figure 7C,D), furthermore they mostly come from the van der Waals interactions of the sidechains in Y132 and I179 with the ring R2 of 60X (Table 4). Two residues V130 and W137 produce the interactions of -1.28 and -1.29 kcal/mol with 60X (Figure 7C,D), which is in good agreement with the CH-π interaction of the alkyl group from V130 and the π-π interaction of the hydrophobic ring of W137 with the ring R3 of 60X (Figure 8E). The distances of the mass centers for the sidechains of V130 and W137 away from that of the ring R3 in 60X are separately distributed at 4.32 and 6.13 Å (Figure 9C), implying the existence of the interactions of V130 and W137 with 60X. More importantly, the energy contributions of V130 and W137 to the 60X-BACE1 association are mainly provided by the van der Waals interactions of the sidechains of V130 and W137 with 60X (Table 4). Besides, the carbonyl group of the residue D93 not only produces the CH-O interactions with the ring R1 of 60X but also forms four HBIs with the occupancy higher than 48.7% with the ring R2 of 60X (Figure 8F and Table 5). The distance between the mass center of the carbonyl group of D93 and that of the ring R2 in 60X is distributed at 4.72 Å (Figure 9C), implying the existence of the CH-O interactions of D93 with 60X.
Based on the aforementioned description, three inhibitors form hydrophobic interactions with L91, S96, V130, Q134, W137, F169 and I179 and the energy contributions of these residues to inhibitor binding mostly come from the interactions of their sidechains with inhibitors. The residue D93 produces four HBIs with inhibitors and these HBIs are formed between the carbonyl group (the sidechain) of D93 and inhibitors. It is concluded that the sidechains of the above mentioned residues play key roles in binding of inhibitors to BACE1. Therefore, it is of high significance to rationally optimize the interactions of inhibitors with the sidechains of key residues in BACE1 for design of efficient inhibitors toward BACE1.

3. Materials and Methods

3.1. Construction of Initial Systems

The initial structures of 60W-, 954- and 60X-BACE1 complexes were obtained from the PDB. The ID 5HDU, 5HDZ and 5HE7 respectively correspond to the 60W-, 954- and 60X-BACE1 complexes. The apo BACE1 without associations of inhibitors were obtained by cutting 60W from the crystal structure 5HDU. The missing residues in three crystal structures were repaired by using the program Modeller [68]. All crystal water and non-inhibitor molecules were deleted from the initial model. The protonated states of residues from BACE1 were checked through the program H++ 3.0 [69]. Then, the following tasks were accomplished with the help of the module Leap in Amber 20 [70,71]: (1) the force field parameters of BACE1 were assigned by employing the ff19SB force field [72], (2) three disulfide bonds were established between C176 and C380, C238 and C403, and C290 and C340, respectively, (3) an octahedral TIP3P water molecule periodic box with a buffer of 10.0 Å was constructed to solve four BACE1-related systems, (4) counter ions were added within the systems in 0.15 M salt environment to neutralize each system, in which the force parameters of sodium ions (Na+) and chloride ions (Cl-) arise from the work of Joung and Cheatham [73,74]. The molecular structures of the three inhibitors 60W, 954 and 60X were optimized at a semiempirical AM1 level, and then, the BCC charges [75,76] were given to each atom of the inhibitors using the Antechamber module in Amber [77]. The general Amber force field (GAFF2) [78,79] was adopted to generate the force field parameters of three inhibitors 60W, 954 and 60X.

3.2. MD Simulations

To remove possible high energy contacts between atoms formed during the initial process of four simulated systems and relieve the instability of systems, two-step energy minimizations were implemented before a real MSMD simulation, composed of a 50000-step steepest descent optimization and a 50000-step conjugate gradient one. Subsequently, all systems were provided a slow heating process from 0 to 300 k in 1 ns in the canonical ensemble (NVT), in which all non-hydrogen atoms in BACE1 and inhibitors were constrained in a weak harmonic restriction of 2 kcalmol−1‧Å2. Then, a 2-ns equilibrium phase was executed on four BACE1-related systems at 300 K under the isothermal−isobaric ensemble (NPT) to further optimize systems. Finally, 600-ns MD simulations were conducted on each system to deeply relax the systems. The above mentioned simulations processes were repeatedly performed for three times, and in which the initial atomic velocities were produced by means of the Maxwell distribution. As a result, three separate MD simulations were completed. Through the aforementioned simulation stages, the Langevin thermostat [80] was utilized to adjust the system temperature and meanwhile the collision frequency was set as 2 ps-1. The shake algorithm [81] was applied to constrain all chemical bonds involved in hydrogen atoms. The long-range electrostatic interactions between atoms were estimated with the particle mesh Ewald algorithm [82] with a cutoff value of 9 Å. At the same time, this cutoff was also employed to calculate van der Waals interactions between atoms. The equilibrium phases of three separate MD trajectories were connected into a single MD trajectory (SMT) to facilitate the post-processing analysis. For this work, all simulations were performed by employing the program pmemd.cuda inlayed in Amber 20 [83,84].

3.3. Calculations of Solvated Interaction Energy

The SIE method can be used to fast and rationally predict binding free energies of inhibitors to targets. The SIE function to calculate inhibitor-BACE1 binding free energy was expressed as the following equation (1)
G b i n d ρ , D i n , α , γ , C = α × E c D i n + G R + E v d W + γ × M S A ρ + C
in which E c and E v d W indicate the intermolecular Coulomb and van der Waals interaction energies between atoms in the bound state, individually, and they were calculated using the Amber molecular mechanics force field ff19SB. The component G R represents the alteration in the reaction field energy caused by binding of inhibitors to BACE1 and was calculated by solving the Poisson equation using the boundary element method BRI BEM [85,86] together with a variable-radius solvent probe [87]. The term γ × M S A is used to reflect the change of free energies relating with the molecular surface area caused by inhibitor binding. The parameters ρ , D i n , γ and C represent Amber van der Waals radii linear scaling coefficient, the solute interior dielectric constant, the molecular surface area coefficient and a constant, respectively. The parameter α indicates the global proportionality coefficient related to the loss of conformational entropy upon binding [88]. The optimized values of the aforementioned parameters are α = 0.1048 , ρ = 1.1 , D i n = 2.25 , γ = 0.0129 kcal·mol-1·Å-2 and C = 2.89 kcal·mol-1 [49]. The SIE calculations were implemented by means of the program Sietraj [67].

3.4. MM-GBSA Calculations

MM-PB/GBSA methods have been widely used to calculate inhibitor-target binding free energies. Hou’s group have been performed a series of works to check the performance of these two methods [89,90,91]. According to their tests, MM-GBSA method was selected to calculate inhibito-BACE1 binding free energies with the following equation 2
G b i n d = E e l e + E v d w + G g b + G s u r f T S
where E e l e and E v d w respectively represent the electrostatic and van der Waals interactions of inhibitors with BACE1, while G p o l and G n o n p o l separately indicate the polar and nonpolar contributions to solvent free energy of the inhibitor-BACE1 complexes. The two components E e l e and E v d w were obtained from the Amber force field ff19SB force. The term G s u r f was estimated by using the empirical equation G s u r f = γ × S A S A + β , in which S A S A denotes the solvent accessible surface area. The G p o l was computed by using the generalized Born (GB) model [92]. For this current study, we selected different GB models that are individually represented by IGB=1, 2, 5 and 66 [92,93] to calculate the G g b so that we can examine impacts of different GB model on calculations of binding free energies. The type of GB models and the corresponding parameters, including radii types, γ and β, were provided in Table 1. The last component -T∆S represents the contribution of the entropic change to binding free energies and was estimated through the mmpbsa_py_nabnmode program in Amber 20 [94]. In our current work, 500 snapshots were extracted from the SMT to compute binding free energies. Since the entropy calculation is too expensive, 100 snapshots picked from the above mentioned 500 snapshots were utilized to calculate the entropy contributions to inhibitor-BACE1 binding.

3.5. Principal Component Analysis

PCA is helpful for our insights into concerted motions of structure domains in BACE1. Hence PCA was executed to clarify how binding of inhibitor impacts concerted motions of BACE1. In this work, PCA was realized through the diagonalization on a covariance matrix C constructed with the coordinates of the Cα atoms in BACE1 based on the equation (3)
C = < ( q i < q i > ) ( q j < q j > ) T >
where, the terms q i and q j are the Cartesian coordinates of the ith and jth Cαatoms in BACE1, respectively, while the terms < q i > and < q j > are their averaged positions on conformational ensembles recorded at the SMT. In general, this average is computed by performing superimposition of the SMT with a referenced structure to abolish overall translations and rotations by using a least-square fitting procedure [95]. The eigenvalues and eigenvectors stemming from the PCA are usually applied to respectively embody the fluctuation amplitude along an eigenvector and concerted motions of the structural domains. We completed the PCA through the program CPPTRAJ [96] inlayed in Amber 20 in this study.

3.6. Dynamics Cross-correlation Map

DCCMs are an efficient approach to explore internal dynamics of targets [97,98,99]. To clarify influence of inhibitor binding on internal dynamics of BACE1, DCCMs were calculated by using the coordinates of the atoms Cα in BACE1 saved at the SMT through the following equation (4).
C i j = < r i · r j > ( < r i 2 > < r j 2 > ) 1 / 2
in which two components r i and r j are the displacement of the Cα atoms i and j relative to their corresponding averaged positions. The angle brackets is an indicator of ensemble averages on the snapshots kept at the SMT. The element values ( C i j ) of DCCMs fluctuate at a range of -1 to 1. The positive and negative C i j values respectively correspond to the positively correlated movements (PCMs) and the anti-correlated motions (ACMs) between the Cα atoms i and j of BACE1. The color-coded bars were employed to characterize the extent of correlated motions between residues of BACE1. In this study, the calculations of DCCMs were finished by using the program in Amber 20.

4. Conclusions

BACE1 plays an important role in the production of the toxic amyloid-β peptides that cause ADs. Insights into inhibitor-BACE1 binding mechanism and conformational changes of BACE1 due to inhibitor binding are significant for development of efficient drugs targeting BACE1. Three separate MD simulations with a total simulation time of 1.8 μs, each for running 600 ns, were respectively conducted on the apo, 60W-, 954- and 60X-binding BACE1. The analyses of RMSDs and RMSFs verify that inhibitor binding not only affects structural stability but also changes the structural flexibility of BACE1. The results from the calculations of DCCMs and PCA indicate that inhibitor binding alters correlated motions and dynamics behavior of BACE1. Binding free energies calculated through the SIE and MM-GBSA methods comparatively reveal that hydrophobic interactions drive inhibitor-BACE1 binding, which should be paid special attentions in future drug design toward BACE1. The inhibitor-residue spectrum from calculations of residue-based free energy decomposition shows that the sidechains of residues L91, D93, S96, V130, Q134, W137, F169 and I179 provide main energy contributions for binding of inhibitors to BACE1. This study is also expected to contribute useful information to development of potent inhibitors toward BACE1.

Supplementary Materials

The following supporting information can be downloaded at the website of this paper posted on Preprints.org. Figure S1. Root-mean-square deviations (RMSDs) of backbone atoms in HSP90 in three independent simulations; Figure S2. Structural domains of HSP90 with obvious changes of RMSF, in which the D1, D2 and D3 display different domains with colors; Figure S3. The function of eigenvalues as eigenvector indexes from principal component analysis, which is used to describe structural fluctuations of HSP90 along the eigenvectors; Figure S4. Free energy landscapes and the repre-sentative structures of the apo HSP90.

Author Contributions

Conceptualization, W.H., J.C. and B.W.; methodology, W.H., J.C. and B.W.; software, Y.W. and F.Y.; validation, Y.W., F.Y. and D.Y..; formal analysis, Y.W., Y.Z. and Y.Z.; inves-tigation,Y.Z..; data curation, Y.W. and D.Y., Y.Z..; writing—original draft preparation, Y.W.; writing—review and editing, Y.W. and D.Y.; visualization, Y.Z..; supervision, Y.W. and D.Y.; project administration, W.H., J.C. and B.W.; funding acquisition, W.H. and J.C. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by high-caliber talent of tuojiang scholar from Shandong Jiaotong University (No. TJXZ202203 and TJXZ202204) and Natural Science Foundation of Shan-dong Province Grant (ZR2019MA040, ZR2021MA069 and ZR2020ME231).

Institutional Review Board Statement

Not applicable.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Sample Availability

Samples of the compounds are not available from the authors.

References

  1. Goedert, M.; Spillantini, M.G. A Century of Alzheimer's Disease. Science 2006, 314, 777–781. [Google Scholar] [CrossRef]
  2. 2017 Alzheimer's disease facts and figures. Alzheimer's & Dementia 2017, 13, 325–373.
  3. Mendiola-Precoma, J.; Berumen, L. C.; Padilla, K.; Garcia-Alcocer, G. Therapies for Prevention and Treatment of Alzheimer’s Disease. BioMed Research International 2016, 2016, 2589276. [Google Scholar] [CrossRef]
  4. Armstrong, R.A. The molecular biology of senile plaques and neurofibrillary tangles in Alzheimer’s disease. Folia Neu ropathol. 2009, 47, 289–299. [Google Scholar]
  5. Vassar, R. BACE1 inhibitor drugs in clinical trials for Alzheimer’s disease. Alzheimer's Res. Ther. 2014, 6, 89. [Google Scholar] [CrossRef]
  6. Sadleir, K. R.; Kandalepas, P. C.; Buggia-Prévot, V.; Nicholson, D. A.; Thinakaran, G.; Vassar, R. Presynaptic dystrophic neurites surrounding amyloid plaques are sites of microtubule disruption, BACE1 elevation, and increased Aβ generation in Alzheimer’s disease. Acta Neuropathol. 2016, 132, 235–256. [Google Scholar] [CrossRef] [PubMed]
  7. Fobare, W. F.; Solvibile, W. R.; Robichaud, A. J.; Malamas, M. S.; Manas, E.; Turner, J.; Hu, Y.; Wagner, E.; Chopra, R.; Cowling, R. Thiophene substituted acylguanidines as BACE1 inhibitors. Bioorg. Med. Chem. Lett. 2007, 17, 5353–5356. [Google Scholar] [CrossRef]
  8. Malamas, M. S.; Barnes, K.; Johnson, M.; Hui, Y.; Zhou, P.; Turner, J.; Hu, Y.; Wagner, E.; Fan, K.; Chopra, R. Di-substituted pyridinyl aminohydantoins as potent and highly selective human β-secretase (BACE1) inhibitors. Bioorg. Med. Chem. 2010, 18, 630–639. [Google Scholar] [CrossRef] [PubMed]
  9. Jordan, J. B.; Whittington, D. A.; Bartberger, M. D.; Sickmier, E. A.; Chen, K.; Cheng, Y.; Judd, T. Fragment-Linking Approach Using 19F NMR Spectroscopy To Obtain Highly Potent and Selective Inhibitors of β-Secretase. J. Med. Chem. 2016, 59, 3732–3749. [Google Scholar] [CrossRef]
  10. Vassar, R.; Kovacs, D. M.; Yan, R.; Wong, P. C. The β-Secretase Enzyme BACE in Health and Alzheimer's Disease: Regulation, Cell Biology, Function, and Therapeutic Potential. J. Neurosci. 2009, 29, 12787–12794. [Google Scholar] [CrossRef]
  11. Zou, Y.; Li, L.; Chen, W.; Chen, T.; Ma, L.; Wang, X.; Xiong, B.; Xu, Y.; Shen, J. Virtual Screening and Structure-Based Discovery of Indole Acylguanidines as Potent β-secretase (BACE1) Inhibitors. Molecules 2013, 18, 5706–5722. [Google Scholar] [CrossRef]
  12. Malamas, M. S.; Erdei, J.; Gunawan, I.; Turner, J.; Hu, Y.; Wagner, E.; Fan, K.; Chopra, R.; Olland, A.; Bard, J. Design and Synthesis of 5,5′-Disubstituted Aminohydantoins as Potent and Selective Human β-Secretase (BACE1) Inhibitors. J. Med. Chem. 2010, 53, 1146–1158. [Google Scholar] [CrossRef]
  13. Xu, Y.; Li, M.-j.; Greenblatt, H.; Chen, W.; Paz, A.; Dym, O.; Peleg, Y.; Chen, T.; Shen, X.; He, J. Flexibility of the flap in the active site of BACE1 as revealed by crystal structures and molecular dynamics simulations. Acta Crystallogr. D 2012, 68, 13–25. [Google Scholar] [CrossRef]
  14. Ruderisch, N.; Schlatter, D.; Kuglstatter, A.; Guba, W.; Huber, S.; Cusulin, C.; Benz, J.; Rufer, A. C.; Hoernschemeyer, J.; Schweit zer, C. Potent and Selective BACE-1 Peptide Inhibitors Lower Brain Aβ Levels Mediated by Brain Shuttle Transport. eBioMedicine 2017, 24, 76–92. [Google Scholar] [CrossRef]
  15. Fujimoto, K.; Matsuoka, E.; Asada, N.; Tadano, G.; Yamamoto, T.; Nakahara, K.; Fuchino, K.; Ito, H.; Kanegawa, N.; Moechars, D. Structure-Based Design of Selective β-Site Amyloid Precursor Protein Cleaving Enzyme 1 (BACE1) Inhibitors: Targeting the Flap to Gain Selectivity over BACE2. J. Med. Chem. 2019, 62, 5080–5095. [Google Scholar] [CrossRef] [PubMed]
  16. Chen, J.; Wang, J.; Yin, B.; Pang, L.; Wang, W.; Zhu, W. Molecular Mechanism of Binding Selectivity of Inhibitors toward BACE1 and BACE2 Revealed by Multiple Short Molecular Dynamics Simulations and Free-Energy Predictions. ACS Chem. Neurosci. 2019, 10, 4303–4318. [Google Scholar] [CrossRef] [PubMed]
  17. Johansson, P.; Kaspersson, K.; Gurrell, I. K.; Bäck, E.; Eketjäll, S.; Scott, C. W.; Cebers, G.; Thorne, P.; McKenzie, M. J.; Beaton, H. Toward β-Secretase-1 Inhibitors with Improved Isoform Selectivity. J. Med. Chem. 2018, 61, 3491–3502. [Google Scholar] [CrossRef] [PubMed]
  18. Oehlrich, D.; Prokopcova, H.; Gijsen, H. J. M. The evolution of amidine-based brain penetrant BACE1 inhibitors. Bioorg. Med. Chem. Lett. 2014, 24, 2033–2045. [Google Scholar] [CrossRef]
  19. Cebers, G.; Alexander, R. C.; Haeberlein, S. B.; Han, D.; Goldwater, R.; Ereshefsky, L.; Olsson, T.; Ye, N.; Rosen, L.; Russell, M. AZD3293: Pharmacokinetic and&nbsp;Pharmacodynamic Effects in&nbsp;Healthy&nbsp;Subjects and Patients with&nbsp;Alzheimer’s Disease. J. Alzheimer's Dis. 2017, 55, 1039–1053. [Google Scholar]
  20. Chen, J.; Yin, B.; Wang, W.; Sun, H. Effects of Disulfide Bonds on Binding of Inhibitors to β-Amyloid Cleaving Enzyme 1 De coded by Multiple Replica Accelerated Molecular Dynamics Simulations. ACS Chem. Neurosci. 2020, 11, 1811–1826. [Google Scholar] [CrossRef]
  21. Cheng, Y.; Judd, T. C.; Bartberger, M. D.; Brown, J.; Chen, K.; Fremeau, R. T., Jr.; Hickman, D.; Hitchcock, S. A.; Jordan, B.; Li, V. From Fragment Screening to In Vivo Efficacy: Optimization of a Series of 2-Aminoquinolines as Potent Inhibitors of Beta-Site Amyloid Precursor Protein Cleaving Enzyme 1 (BACE1). J. Med. Chem. 2011, 54, 5836–5857. [Google Scholar] [CrossRef] [PubMed]
  22. Butler, C. R.; Ogilvie, K.; Martinez-Alsina, L.; Barreiro, G.; Beck, E. M.; Nolan, C. E.; Atchison, K.; Benvenuti, E.; Buzon, L.; Doran, S. Aminomethyl-Derived Beta Secretase (BACE1) Inhibitors: Engaging Gly230 without an Anilide Functionality. J. Med. Chem. 2017, 60, 386–402. [Google Scholar] [CrossRef] [PubMed]
  23. Koriyama, Y.; Hori, A.; Ito, H.; Yonezawa, S.; Baba, Y.; Tanimoto, N.; Ueno, T.; Yamamoto, S.; Yamamoto, T.; Asada, N. Discovery of Atabecestat (JNJ-54861911): A Thiazine-Based β-Amyloid Precursor Protein Cleaving Enzyme 1 Inhibitor Advanced to the Phase 2b/3 EARLY Clinical Trial. J. Med. Chem. 2021, 64, 1873–1888. [Google Scholar] [CrossRef]
  24. Rueeger, H.; Lueoend, R.; Machauer, R.; Veenstra, S. J.; Holzer, P.; Hurth, K.; Voegtle, M.; Frederiksen, M.; Rondeau, J.-M.; Tintelnot-Blomley, M. Synthesis of the Potent, Selective, and Efficacious β-Secretase (BACE1) Inhibitor NB-360. J. Med. Chem. 2021, 64, 4677–4696. [Google Scholar] [CrossRef]
  25. Bao, H. Y.; Wang, W.; Sun, H. B.; Chen, J. Z. Binding modes of GDP, GTP and GNP to NRAS deciphered by using Gaussian accelerated molecular dynamics simulations. SAR QSAR in Environ. Res. 2023, 34, 65–89. [Google Scholar] [CrossRef] [PubMed]
  26. Xue, W.; Wang, P.; Tu, G.; Yang, F.; Zheng, G.; Li, X.; Li, X.; Chen, Y.; Yao, X.; Zhu, F. Computational identification of the binding mechanism of a triple reuptake inhibitor amitifadine for the treatment of major depressive disorder. Phys. Chem. Chem. Phys. 2018, 20, 6606–6616. [Google Scholar] [CrossRef]
  27. Wang, J.; Arantes, P. R.; Bhattarai, A.; Hsu, R. V.; Pawnikar, S.; Huang, Y.-m. M.; Palermo, G.; Miao, Y. Gaussian accelerated molecular dynamics: Principles and applications. WIREs Comput. Mol. Sci. 2021, 11, e1521. [Google Scholar] [CrossRef]
  28. Sun, Z.; Gong, Z.; Xia, F.; He, X. Ion dynamics and selectivity of Nav channels from molecular dynamics simulation. Chem. Phys. 2021, 548, 111245. [Google Scholar] [CrossRef]
  29. Shao, Q.; Xiong, M.; Li, J.; Hu, H.; Su, H.; Xu, Y. Unraveling the catalytic mechanism of SARS-CoV-2 papain-like protease with allosteric modulation of C270 mutation using multiscale computational approaches. Chem. Sci. 2023. [Google Scholar] [CrossRef]
  30. Chen, Z.; Zhang, X.; Peng, C.; Wang, J.; Xu, Z.; Chen, K.; Shi, J.; Zhu, W. D3Pockets: A Method and Web Server for Systematic Analysis of Protein Pocket Dynamics. J. Chem. Inf. Model. 2019, 59, 3353–3358. [Google Scholar] [CrossRef]
  31. Xue, W.; Yang, F.; Wang, P.; Zheng, G.; Chen, Y.; Yao, X.; Zhu, F. What Contributes to Serotonin–Norepinephrine Reuptake Inhibitors’ Dual-Targeting Mechanism? The Key Role of Transmembrane Domain 6 in Human Serotonin and Norepinephrine Transporters Revealed by Molecular Dynamics Simulation. ACS Chem. Neurosci. 2018, 9, 1128–1140. [Google Scholar] [CrossRef] [PubMed]
  32. Sun, Z.; He, Q.; Gong, Z.; Kalhor, P.; Huai, Z.; Liu, Z. A General Picture of Cucurbit[8]uril Host&ndash;Guest Binding: Recali brating Bonded Interactions. Molecules 2023, 28, 3124. [Google Scholar] [PubMed]
  33. Chen, J.; Zeng, Q.; Wang, W.; Sun, H.; Hu, G. Decoding the Identification Mechanism of an SAM-III Riboswitch on Ligands through Multiple Independent Gaussian-Accelerated Molecular Dynamics Simulations. J. Chem. Inf. Model. 2022, 62, 6118–6132. [Google Scholar] [CrossRef] [PubMed]
  34. Hou, T.; Yu, R. Molecular Dynamics and Free Energy Studies on the Wild-type and Double Mutant HIV-1 Protease Complexed with Amprenavir and Two Amprenavir-Related Inhibitors:  Mechanism for Binding and Drug Resistance. J. Med. Chem. 2007, 50, 1177–1188. [Google Scholar] [CrossRef] [PubMed]
  35. Amadei, A.; Linssen, A. B. M.; Berendsen, H.J.C. Essential dynamics of proteins. Proteins 1993, 17, 412–425. [Google Scholar] [CrossRef]
  36. Levy, R. M.; Srinivasan, A. R.; Olson, W. K.; McCammon, J.A. Quasi-harmonic method for studying very low frequency modes in proteins. Biopolymers 1984, 23, 1099–1112. [Google Scholar] [CrossRef]
  37. Bao, H.; Wang, W.; Sun, H.; Chen, J. Probing mutation-induced conformational transformation of the GTP/M-RAS complex through Gaussian accelerated molecular dynamics simulations. J. Enzym. Inhib. Med. Ch. 2023, 38, 2195995. [Google Scholar] [CrossRef]
  38. Bao, H.; Wang, W.; Sun, H.; Chen, J. The switch states of the GDP-bound HRAS affected by point mutations: a study from Gaussian accelerated molecular dynamics simulations and free energy landscapes. J. Biomol. Struct. Dyn. 2023. [CrossRef]
  39. Auffinger, P.; Westhof, E. RNA hydration: three nanoseconds of multiple molecular dynamics simulations of the solvated tRNAAsp anticodon hairpin. J. Mol. Biol. 1997, 269, 326–341. [Google Scholar] [CrossRef]
  40. Wang, L.; Wang, Y.; Yu, Y.; Liu, D.; Zhao, J.; Zhang, L. Deciphering Selectivity Mechanism of BRD9 and TAF1(2) toward Inhib itors Based on Multiple Short Molecular Dynamics Simulations and MM-GBSA Calculations. Molecules 2023, 28, 2583. [Google Scholar] [CrossRef]
  41. Chen, J.; Liu, X.; Zhang, S.; Chen, J.; Sun, H.; Zhang, L.; Zhang, Q. Molecular mechanism with regard to the binding selectivity of inhibitors toward FABP5 and FABP7 explored by multiple short molecular dynamics simulations and free energy analyses. Phys. Chem. Chem. Phys. 2020, 22, 2262–2275. [Google Scholar] [CrossRef] [PubMed]
  42. Suruzhon, M.; Bodnarchuk, M. S.; Ciancetta, A.; Viner, R.; Wall, I. D.; Essex, J.W. Sensitivity of Binding Free Energy Calculations to Initial Protein Crystal Structure. J. Chem. Theory Comput. 2021, 17, 1806–1821. [Google Scholar] [CrossRef] [PubMed]
  43. Wang, R.; Zheng, Q. Multiple Molecular Dynamics Simulations and Free-Energy Predictions Uncover the Susceptibility of Var iants of HIV-1 Protease against Inhibitors Darunavir and KNI-1657. Langmuir 2021, 37, 14407–14418. [Google Scholar] [CrossRef] [PubMed]
  44. Caves, L. S. D.; Evanseck, J. D.; Karplus, M. Locally accessible conformations of proteins: Multiple molecular dynamics simula tions of crambin. Protein Sci. 1998, 7, 649–666. [Google Scholar] [CrossRef] [PubMed]
  45. Knapp, B.; Ospina, L.; Deane, C. M. Avoiding False Positive Conclusions in Molecular Simulation: The Importance of Replicas. J. Chem. Theory Comput. 2018, 14, 6127–6138. [Google Scholar] [CrossRef]
  46. Wang, W.; Kollman, P. A. Free energy calculations on dimer stability of the HIV protease using molecular dynamics and a continuum solvent model11Edited by B. Honig. J. Mol. Biol. 2000, 303, 567–582. [Google Scholar] [CrossRef] [PubMed]
  47. Wang, W.; Kollman, P. A. Computational study of protein specificity: The molecular basis of HIV-1 protease drug resistance. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 14937–14942. [Google Scholar] [CrossRef]
  48. Wang, J.; Morin, P.; Wang, W.; Kollman, P. A. Use of MM-PBSA in Reproducing the Binding Free Energies to HIV-1 RT of TIBO Derivatives and Predicting the Binding Mode to HIV-1 RT of Efavirenz by Docking and MM-PBSA. J. Am. Chem. Soc. 2001, 123, 5221–5230. [Google Scholar] [CrossRef]
  49. Naïm, M.; Bhat, S.; Rankin, K. N.; Dennis, S.; Chowdhury, S. F.; Siddiqi, I.; Drabik, P.; Sulea, T.; Bayly, C. I.; Jakalian, A. Solvated Interaction Energy (SIE) for Scoring Protein−Ligand Binding Affinities. 1. Exploring the Parameter Space. J. Chem. Inf. Model. 2007, 47, 122–133. [Google Scholar] [CrossRef]
  50. Tzoupis, H.; Leonis, G.; Mavromoustakos, T.; Papadopoulos, M. G. A Comparative Molecular Dynamics, MM–PBSA and Ther modynamic Integration Study of Saquinavir Complexes with Wild-Type HIV-1 PR and L10I, G48V, L63P, A71V, G73S, V82A and I84V Single Mutants. J. Chem. Theory Comput. 2013, 9, 1754–1764. [Google Scholar] [CrossRef]
  51. Chen, J.; Wang, X.; Zhu, T.; Zhang, Q.; Zhang, J. Z. H. A Comparative Insight into Amprenavir Resistance of Mutations V32I, G48V, I50V, I54V, and I84V in HIV-1 Protease Based on Thermodynamic Integration and MM-PBSA Methods. J. Chem. Inf. Model. 2015, 55, 1903–1913. [Google Scholar] [CrossRef] [PubMed]
  52. Leonis, G.; Steinbrecher, T.; Papadopoulos, M. G. A Contribution to the Drug Resistance Mechanism of Darunavir, Amprenavir, Indinavir, and Saquinavir Complexes with HIV-1 Protease Due to Flap Mutation I50V: A Systematic MM–PBSA and Thermodynamic Integration Study. J. Chem. Inf. Model. 2013, 53, 2141–2153. [Google Scholar] [CrossRef] [PubMed]
  53. Aldeghi, M.; Heifetz, A.; Bodkin, M. J.; Knapp, S.; Biggin, P. C. Predictions of Ligand Selectivity from Absolute Binding Free Energy Calculations. J. Am. Chem. Soc. 2017, 139, 946–957. [Google Scholar] [CrossRef] [PubMed]
  54. Shirts, M. R.; Pitera, J. W.; Swope, W. C.; Pande, V. S. Extremely precise free energy calculations of amino acid side chain analogs: Comparison of common molecular mechanics force fields for proteins. J. Chem. Phys. 2003, 119, 5740–5761. [Google Scholar] [CrossRef]
  55. Chen, J.; Wang, X.; Pang, L.; Zhang, J. Z. H.; Zhu, T. Effect of mutations on binding of ligands to guanine riboswitch probed by free energy perturbation and molecular dynamics simulations. Nucleic Acids Res. 2019, 47, 6618–6631. [Google Scholar] [CrossRef] [PubMed]
  56. Aldeghi, M.; Heifetz, A.; Bodkin, M. J.; Knapp, S.; Biggin, P. C. Accurate calculation of the absolute free energy of binding for drug molecules. Chem. Sci. 2016, 7, 207–218. [Google Scholar] [CrossRef] [PubMed]
  57. Saravanan, K.; Sugarthi, S.; Suganya, S.; Kumaradhas, P. Probing the intermolecular interactions, binding affinity, charge den sity distribution and dynamics of silibinin in dual targets AChE and BACE1: QTAIM and molecular dynamics perspective. J. Biomol. Struct. Dyn. 2022, 40, 12880–12894. [Google Scholar] [CrossRef] [PubMed]
  58. Ellis, C. R.; Tsai, C.-C.; Hou, X.; Shen, J. Constant pH Molecular Dynamics Reveals pH-Modulated Binding of Two Small-Mol ecule BACE1 Inhibitors. J. Phys. Chem. Lett. 2016, 7, 944–949. [Google Scholar] [CrossRef]
  59. Bao, L.-Q.; Baecker, D.; Mai Dung, D. T.; Phuong Nhung, N.; Thi Thuan, N.; Nguyen, P. L.; Phuong Dung, P. T.; Huong, T. T. L.; Rasulev, B.; Casanola-Martin, G. M. Development of Activity Rules and Chemical Fragment Design for In Silico Discovery of AChE and BACE1 Dual Inhibitors against Alzheimer's Disease. Molecules 2023, 28, 3588. [Google Scholar] [CrossRef]
  60. Hernández-Rodríguez, M.; Correa-Basurto, J.; Gutiérrez, A.; Vitorica, J.; Rosales-Hernández, M. C. Asp32 and Asp228 determine the selective inhibition of BACE1 as shown by docking and molecular dynamics simulations. Eur. J. Med. Chem. 2016, 124, 1142–1154. [Google Scholar] [CrossRef]
  61. Chen, J.; Zhang, S.; Wang, W.; Sun, H.; Zhang, Q.; Liu, X. Binding of Inhibitors to BACE1 Affected by pH-Dependent Protonation: An Exploration from Multiple Replica Gaussian Accelerated Molecular Dynamics and MM-GBSA Calculations. ACS Chem. Neurosci. 2021, 12, 2591–2607. [Google Scholar] [CrossRef]
  62. Hatmal, M. m. M.; Jaber, S.; Taha, M. O. Combining molecular dynamics simulation and ligand-receptor contacts analysis as a new approach for pharmacophore modeling: beta-secretase 1 and check point kinase 1 as case studies. J. Comput. Aid. Mol. Des. 2016, 30, 1149–1163. [Google Scholar] [CrossRef] [PubMed]
  63. Mandal, M.; Wu, Y.; Misiaszek, J.; Li, G.; Buevich, A.; Caldwell, J. P.; Liu, X.; Mazzola, R. D.; Orth, P.; Strickland, C. Struc ture-Based Design of an Iminoheterocyclic β-Site Amyloid Precursor Protein Cleaving Enzyme (BACE) Inhibitor that Lowers Central Aβ in Nonhuman Primates. J. Med. Chem. 2016, 59, 3231–3248. [Google Scholar] [CrossRef] [PubMed]
  64. Kabsch, W.; Sander, C. Dictionary of protein secondary structure: Pattern recognition of hydrogen-bonded and geometrical features. Biopolymers 1983, 22, 2577–2637. [Google Scholar] [CrossRef]
  65. Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graph. 1996, 14, 33–38. [Google Scholar] [CrossRef] [PubMed]
  66. Wang, Y.-T.; Su, Z.-Y.; Hsieh, C.-H.; Chen, C.-L. Predictions of Binding for Dopamine D2 Receptor Antagonists by the SIE Method. J. Chem. Inf. Model. 2009, 49, 2369–2375. [Google Scholar] [CrossRef] [PubMed]
  67. Cui, Q.; Sulea, T.; Schrag, J. D.; Munger, C.; Hung, M.-N.; Naïm, M.; Cygler, M.; Purisima, E. O. Molecular Dynamics—Solvated Interaction Energy Studies of Protein–Protein Interactions: The MP1–p14 Scaffolding Complex. J. Mol. Biol. 2008, 379, 787–802. [Google Scholar] [CrossRef] [PubMed]
  68. Webb, B.; Sali, A. Comparative Protein Structure Modeling Using MODELLER. Current Protocols in Bioinformatics 2014, John Wiley & Sons, Inc., 2014; pp 2015.2016.2011-2015.2016.2032.
  69. Anandakrishnan, R.; Aguilar, B.; Onufriev, A. V. H++ 3.0: automating pK prediction and the preparation of biomolecular struc tures for atomistic molecular modeling and simulations. Nucleic Acids Res. 2012, 40, W537–W541. [Google Scholar] [CrossRef] [PubMed]
  70. Salomon-Ferrer, R.; Case, D. A.; Walker, R. C. An overview of the Amber biomolecular simulation package. WIREs Comput. Mol. Sci. 2013, 3, 198–210. [Google Scholar] [CrossRef]
  71. Case, D. A.; Cheatham III, T. E.; Darden, T.; Gohlke, H.; Luo, R.; Merz Jr., K. M.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. The Amber biomolecular simulation programs. J. Comput. Chem. 2005, 26, 1668–1688. [Google Scholar] [CrossRef]
  72. Tian, C.; Kasavajhala, K.; Belfon, K. A. A.; Raguette, L.; Huang, H.; Migues, A. N.; Bickel, J.; Wang, Y.; Pincay, J.; Wu, Q. ff19SB: Amino-Acid-Specific Protein Backbone Parameters Trained against Quantum Mechanics Energy Surfaces in Solution. J. Chem. Theory Comput. 2020, 16, 528–552. [Google Scholar] [CrossRef] [PubMed]
  73. Joung, I. S.; Cheatham, T. E., III. Determination of Alkali and Halide Monovalent Ion Parameters for Use in Explicitly Solvated Biomolecular Simulations. J. Phys. Chem. B 2008, 112, 9020–9041. [Google Scholar] [CrossRef] [PubMed]
  74. Joung, I. S.; Cheatham, T. E., III. Molecular Dynamics Simulations of the Dynamic and Energetic Properties of Alkali and Halide Ions Using Water-Model-Specific Ion Parameters. J. Phys. Chem. B 2009, 113, 13279–13290. [Google Scholar] [CrossRef] [PubMed]
  75. Jakalian, A.; Jack, D. B.; Bayly, C. I. Fast, efficient generation of high-quality atomic charges. AM1-BCC model: II. Parameteriza tion and validation. J. Comput. Chem. 2002, 23, 1623–1641. [Google Scholar] [CrossRef] [PubMed]
  76. Jakalian, A.; Bush, B. L.; Jack, D. B.; Bayly, C. I. Fast, efficient generation of high-quality atomic charges. AM1-BCC model: I. Method. J. Comput. Chem. 2000, 21, 132–146. [Google Scholar] [CrossRef]
  77. Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. Automatic atom type and bond type perception in molecular mechanical cal culations. J. Mol. Graph. Model. 2006, 25, 247–260. [Google Scholar] [CrossRef] [PubMed]
  78. Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. J. Comput. Chem. 2004, 25, 1157–1174. [Google Scholar] [CrossRef]
  79. He, X.; Man, V. H.; Yang, W.; Lee, T.-S.; Wang, J. A fast and high-quality charge model for the next generation general AMBER force field. J. Chem. Phys. 2020, 153. [Google Scholar] [CrossRef]
  80. Izaguirre, J. A.; Catarello, D. P.; Wozniak, J. M.; Skeel, R. D. Langevin stabilization of molecular dynamics. J. Chem. Phys. 2001, 114, 2090–2098. [Google Scholar] [CrossRef]
  81. Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23, 327–341. [Google Scholar] [CrossRef]
  82. Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. J. Chem. Phys. 1995, 103, 8577–8593. [Google Scholar] [CrossRef]
  83. Salomon-Ferrer, R.; Götz, A. W.; Poole, D.; Le Grand, S.; Walker, R. C. Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 2. Explicit Solvent Particle Mesh Ewald. J. Chem. Theory Comput. 2013, 9, 3878–3888. [Google Scholar] [CrossRef] [PubMed]
  84. Götz, A. W.; Williamson, M. J.; Xu, D.; Poole, D.; Le Grand, S.; Walker, R. C. Routine Microsecond Molecular Dynamics Simu lations with AMBER on GPUs. 1. Generalized Born. J. Chem. Theory Comput. 2012, 8, 1542–1555. [Google Scholar] [CrossRef] [PubMed]
  85. Purisima, E. O. Fast summation boundary element method for calculating solvation free energies of macromolecules. J. Comput. Chem. 1998, 19, 1494–1504. [Google Scholar] [CrossRef]
  86. Purisima, E. O.; Nilar, S. H. A simple yet accurate boundary element method for continuum dielectric calculations. J. Comput. Chem. 1995, 16, 681–689. [Google Scholar] [CrossRef]
  87. Bhat, S.; Purisima, E. O. Molecular surface generation using a variable-radius solvent probe. Proteins 2006, 62, 244–261. [Google Scholar] [CrossRef]
  88. Perdih, A.; Bren, U.; Solmajer, T. Binding free energy calculations of N-sulphonyl-glutamic acid inhibitors of MurD ligase. J. Mol. Model. 2009, 15, 983–996. [Google Scholar] [CrossRef]
  89. Sun, H.; Li, Y.; Shen, M.; Tian, S.; Xu, L.; Pan, P.; Guan, Y.; Hou, T. Assessing the performance of MM/PBSA and MM/GBSA methods. 5. Improved docking performance using high solute dielectric constant MM/GBSA and MM/PBSA rescoring. Phys. Chem. Chem. Phys. 2014, 16, 22035–22045. [Google Scholar] [CrossRef]
  90. Sun, H.; Li, Y.; Tian, S.; Xu, L.; Hou, T. Assessing the performance of MM/PBSA and MM/GBSA methods. 4. Accuracies of MM/PBSA and MM/GBSA methodologies evaluated by various simulation protocols using PDBbind data set. Phys. Chem. Chem. Phys. 2014, 16, 16719–16729. [Google Scholar] [CrossRef]
  91. Sun, H.; Duan, L.; Chen, F.; Liu, H.; Wang, Z.; Pan, P.; Zhu, F.; Zhang, J. Z. H.; Hou, T. Assessing the performance of MM/PBSA and MM/GBSA methods. 7. Entropy effects on the performance of end-point binding free energy calculation approaches. Phys. Chem. Chem. Phys. 2018, 20, 14450–14460. [Google Scholar] [CrossRef]
  92. Onufriev, A.; Bashford, D.; Case, D. A. Exploring protein native states and large-scale conformational changes with a modified generalized born model. Proteins 2004, 55, 383–394. [Google Scholar] [CrossRef] [PubMed]
  93. Tsui, V.; Case, D. A. Theory and applications of the generalized born solvation model in macromolecular simulations. Biopoly mers 2000, 56, 275–291. [Google Scholar] [CrossRef]
  94. Miller, B. R., III; McGee, T. D., Jr.; Swails, J. M.; Homeyer, N.; Gohlke, H.; Roitberg, A. E. MMPBSA.py: An Efficient Program for End-State Free Energy Calculations. J. Chem. Theory Comput. 2012, 8, 3314–3321. [Google Scholar] [CrossRef] [PubMed]
  95. McLachlan, A. D. Gene duplications in the structural evolution of chymotrypsin. J. Mol. Biol. 1979, 128, 49–79. [Google Scholar] [CrossRef]
  96. Roe, D. R.; Cheatham, T. E., III. PTRAJ and CPPTRAJ: Software for Processing and Analysis of Molecular Dynamics Trajectory Data. J. Chem. Theory Comput. 2013, 9, 3084–3095. [Google Scholar] [CrossRef]
  97. Ichiye, T.; Karplus, M. Collective motions in proteins: A covariance analysis of atomic fluctuations in molecular dynamics and normal mode simulations. Proteins 1991, 11, 205–217. [Google Scholar] [CrossRef]
  98. Chen, J.; Zhang, S.; Wang, W.; Pang, L.; Zhang, Q.; Liu, X. Mutation-Induced Impacts on the Switch Transformations of the GDP- and GTP-Bound K-Ras: Insights from Multiple Replica Gaussian Accelerated Molecular Dynamics and Free Energy Analysis. J. Chem. Inf. Model. 2021, 61, 1954–1969. [Google Scholar] [CrossRef]
  99. Yu, Z.; Su, H.; Chen, J.; Hu, G. Deciphering Conformational Changes of the GDP-Bound NRAS Induced by Mutations G13D, Q61R, and C118S through Gaussian Accelerated Molecular Dynamic Simulations. Molecules 2022, 27, 5596. [Google Scholar] [CrossRef]
Figure 1. Molecular structures: (A) inhibitor-BACE1 complex, in which BACE1 is shown in cartoon and surface modes, while the inhibitor is displayed in stick modes, (B) binding pocket of inhibitor to BACE1, (C), (D) and (E) correspond to 60W, 954 and 60X, respectively, in which three inhibitors are displayed in line modes. In this figure, the shaded areas indicate the groups of inhibitors that possibly produce hydrophobic interactions with BACE1.
Figure 1. Molecular structures: (A) inhibitor-BACE1 complex, in which BACE1 is shown in cartoon and surface modes, while the inhibitor is displayed in stick modes, (B) binding pocket of inhibitor to BACE1, (C), (D) and (E) correspond to 60W, 954 and 60X, respectively, in which three inhibitors are displayed in line modes. In this figure, the shaded areas indicate the groups of inhibitors that possibly produce hydrophobic interactions with BACE1.
Preprints 75229 g001
Figure 2. Frequency of RMSDs and RMSFs: (A) frequency distributions of RMSDs for the apo BACE1, 60W-, 954- and 60X-bound ones, (B) frequency of RMSDs for three inhibitors, (C) RMSFs of the apo and bound states of BACE1 and (D) the difference in RMSFs between the apo BACE1 and the inhibitor-bound ones.
Figure 2. Frequency of RMSDs and RMSFs: (A) frequency distributions of RMSDs for the apo BACE1, 60W-, 954- and 60X-bound ones, (B) frequency of RMSDs for three inhibitors, (C) RMSFs of the apo and bound states of BACE1 and (D) the difference in RMSFs between the apo BACE1 and the inhibitor-bound ones.
Preprints 75229 g002
Figure 3. Stability of secondary structures and structure compact extents of BACE1: (A) time evolution of secondary structure for the apo BACE1 in the simulation 1, (B) time evolution of secondary structure for the apo BACE1 in simulation 2, (C) time evolution of secondary structure for the apo BACE1 in the simulation 3 and (D) the frequency distribution of the BACE1 gyration.
Figure 3. Stability of secondary structures and structure compact extents of BACE1: (A) time evolution of secondary structure for the apo BACE1 in the simulation 1, (B) time evolution of secondary structure for the apo BACE1 in simulation 2, (C) time evolution of secondary structure for the apo BACE1 in the simulation 3 and (D) the frequency distribution of the BACE1 gyration.
Preprints 75229 g003
Figure 4. DCCMs calculated by using the coordinates of the atoms Cα from BACE1: (A) the apo BACE1, (B) the 60W-bound BACE1, (C) the 954-bound BACE1 and (D) the 60X-bound BACE1. In this figure, the color bar is used to reflect the extents of correlated motions between the regions of BACE1.
Figure 4. DCCMs calculated by using the coordinates of the atoms Cα from BACE1: (A) the apo BACE1, (B) the 60W-bound BACE1, (C) the 954-bound BACE1 and (D) the 60X-bound BACE1. In this figure, the color bar is used to reflect the extents of correlated motions between the regions of BACE1.
Preprints 75229 g004
Figure 5. Concerted motions of structural domains from BACE1 revealed by PCA: (A) the apo BACE1, (B) the 60W-bound BACE1, (C) the 954-bound BACE1 and (D) the 60X-bound BACE1.
Figure 5. Concerted motions of structural domains from BACE1 revealed by PCA: (A) the apo BACE1, (B) the 60W-bound BACE1, (C) the 954-bound BACE1 and (D) the 60X-bound BACE1.
Preprints 75229 g005
Figure 6. Free energy landscapes and the representative structures: (A) free energy landscape of the 60W-bound BACE1, (B) structural superimpositions of the 60W-bound BACE1 located at the EV1 and EV2, (C) structural alignment of 60W falling into the EV1 and EV2, (D) free energy landscape of the 954-bound BACE1, (E) structural alignment of the 954-bound BACE1 situated at the EV1 and EV2, (F) structural superimposition of 954 located at the EV1and EV2, (G) free energy landscape of the 60X-bound BACE1, (H) superimposition of the 60X-bound BACE1 and (I) structural alignment of 60X trapped at the EB1-EB3.
Figure 6. Free energy landscapes and the representative structures: (A) free energy landscape of the 60W-bound BACE1, (B) structural superimpositions of the 60W-bound BACE1 located at the EV1 and EV2, (C) structural alignment of 60W falling into the EV1 and EV2, (D) free energy landscape of the 954-bound BACE1, (E) structural alignment of the 954-bound BACE1 situated at the EV1 and EV2, (F) structural superimposition of 954 located at the EV1and EV2, (G) free energy landscape of the 60X-bound BACE1, (H) superimposition of the 60X-bound BACE1 and (I) structural alignment of 60X trapped at the EB1-EB3.
Preprints 75229 g006
Figure 7. Interactions of inhibitors with BACE1: (A) interaction spectrum of 60W with separate residues of BACE1, (B) interaction spectrum of 954 with each residue of BACE1, (C) interaction spectrum of 60X with separate residues of BACE1 and (D) key residues in inhibitor-BACE1 interactions.
Figure 7. Interactions of inhibitors with BACE1: (A) interaction spectrum of 60W with separate residues of BACE1, (B) interaction spectrum of 954 with each residue of BACE1, (C) interaction spectrum of 60X with separate residues of BACE1 and (D) key residues in inhibitor-BACE1 interactions.
Preprints 75229 g007
Figure 8. Geometric information of inhibitor-residue interactions: (A) the hydrophobic interact tions of 60W with residues, (B) the 60W-BACE1 HBIs, (C) the hydrophobic interactions between 954 and residues, (D) the 954-BACE1 HBIs, (E) the hydrophobic interactions of 60X with residues and (F) the 60X-BACE1 HBIs.
Figure 8. Geometric information of inhibitor-residue interactions: (A) the hydrophobic interact tions of 60W with residues, (B) the 60W-BACE1 HBIs, (C) the hydrophobic interactions between 954 and residues, (D) the 954-BACE1 HBIs, (E) the hydrophobic interactions of 60X with residues and (F) the 60X-BACE1 HBIs.
Preprints 75229 g008
Figure 9. Frequency distributions of the distances relating with key inhibitor-residue interactions: (A) the 60W-BACE1 complex, (B) the 954-BACE1 complex and (C) the 60X-BACE1 complex.
Figure 9. Frequency distributions of the distances relating with key inhibitor-residue interactions: (A) the 60W-BACE1 complex, (B) the 954-BACE1 complex and (C) the 60X-BACE1 complex.
Preprints 75229 g009
Table 1. Binding free energies of inhibitors to BACE1 calculated by using SIE methoda.
Table 1. Binding free energies of inhibitors to BACE1 calculated by using SIE methoda.
Components 60W-BACE1 954-BACE1 60X-BACE1
average StdErr average StdErr average StdErr
E v d W -54.53 0.50 -42.68 0.97 -42.50 0.43
E C -16.18 0.35 -12.55 0.42 -15.22 0.29
G R 24.77 0.37 21.27 0.57 22.80 0.37
γ × M S A -10.63 0.09 -8.14 0.18 -8.05 0.06
C -2.89 0.00 -2.89 0.00 -2.89 0.00
b G b i n d -8.82 0.07 -7.30 0.11 -7.39 0.05
c G e x p -12.3 -10.0 -11.4
aAll energy components are scaled in kcal·mol-1; b G b i n d = α × E c + G R + E v d W + γ × M S A + C ; c G e x p is obtained by using G = R T l n K i with the experimental value of K i [63].
Table 2. The parameters used in MM-GBSA calculations with different generalized Born model.
Table 2. The parameters used in MM-GBSA calculations with different generalized Born model.
Parameters IGB=1 IGB=2 IGB=5 IGB=66
a γ 0.0072 0.005 0.005 0.005
a β 0.00 0.00 0.00 0.00
b r a d i i mbondi mbondi2 mbondi2 bondi
Table 3. Binding free energies calculated by using MM-GBSA method with different GB modela.
Table 3. Binding free energies calculated by using MM-GBSA method with different GB modela.
Energy 60W 954 60X
IGB=1 IGB=2 IGB=5 IGB=66 IGB=1 IGB=2 IGB=5 IGB=66 IGB=1 IGB=2 IGB=5 IGB=66
E e l e -36.23 -36.23 -36.23 -36.23 -28.25 -28.25 -28.25 -28.25 -33.88 -33.88 -33.88 -33.88
E v d W -54.58 -54.58 -54.58 -54.58 -42.36 -42.36 -42.36 -42.36 -42.31 -42031 -42.31 -42.31
G g b 46.47 54.13 11.96 69.01 35.94 42.71 5.57 55.77 40.06 46.70 6.31 56.78
G s u r f -7.21 -5.00 -5.00 -5.00 -5.64 -3.92 -3.92 -3.92 -5.71 -3.96 -3.96 -3.96
b G p o l 10.24 17.9 -24.27 32.78 7.69 14.46 -22.98 27.52 6.18 12.82 -27.57 22.9
c G h y d r o -61.79 -59.58 -59.58 -59.58 -48 -46.28 -46.28 -46.28 -48.02 -46.27 -46.27 -46.27
d H -51.55 -41.68 -83.85 -26.8 -40.31 -31.82 -69.26 -18.76 -41.84 -33.45 -73.84 -23.37
T S 22.52 18.60 18.94
G b i n d -29.03 -19.16 -61.33 -4.28 -21.71 -13.22 -50.66 -0.16 -22.9 -14.51 -54.9 -4.43
e G e x p -12.3 -10.0 -11.4
a All free energy components are in scaled in kcal/mol; b G p o l = E e l e + G g b which is used to describe polar interactions of inhibitors with BACE1; c G h y d r o = E v d W + G s u r f which is utilized to signify hydrophobic interactions of inhibitors with BACE1; d H = G p o l + G h y d r o which is adopted to indicate the enthalpy effect during bindings of inhibitors to BACE1; e G e x p is obtained by using G = R T l n K i with the experimental value of K i [63].
Table 4. Contributions of the side chains and backbones to inhibitor-residue interactionsa.
Table 4. Contributions of the side chains and backbones to inhibitor-residue interactionsa.
Inhibitor Residue S v d W B v d W T v d W S e l e B e l e T e l e S g b B g b T g b G
60W L91 -1.35 -0.07 -1.42 0.06 -0.18 -0.12 -0.05 -0.00 -0.05 -1.72
D93 -0.59 -0.19 -0.78 -16.62 -0.45 -17.08 15.42 0.51 15.93 -2.04
S96 -1.69 -0.36 -2.05 1.39 -0.24 1.15 -0.47 0.37 -0.10 -1.13
V130 -1.43 -0.11 -1.54 -0.03 -0.09 -0.12 0.03 0.07 0.10 -1.75
Y132 -0.04 -0.07 -0.11 -0.02 0.15 0.13 0.03 -0.02 0.01 0.03
Q134 -2.96 -0.74 -3.70 -0.26 -0.76 -1.02 0.99 1.52 2.51 -2.58
W137 -1.34 -0.04 -1.38 -0.06 0.03 -0.03 0.59 -0.03 0.56 -0.93
F169 -0.69 -0.14 -0.83 0.11 -0.14 -0.03 0.20 0.25 0.45 -0.45
I179 -1.81 -0.21 -2.02 0.17 -0.05 0.12 -0.11 -0.12 -0.23 -2.26
954 L91 -0.95 -0.05 -1.00 0.07 -0.23 -0.16 -0.04 0.15 0.11 -1.14
D93 -0.37 -0.10 -0.47 -12.67 -0.32 -12.99 12.51 0.28 12.79 -0.77
S96 -1.47 -0.34 -1.81 1.09 0.22 1.31 -0.13 0.00 -0.13 -0.78
V130 -0.72 -0.09 -0.81 -0.02 -0.03 -0.05 0.04 0.13 0.17 -0.77
Y132 -2.32 -0.11 -2.43 -0.99 0.08 -0.91 1.70 0.01 1.71 -1.88
Q134 -0.23 -0.06 -0.29 -0.36 -0.04 -0.40 0.39 0.06 0.45 0.27
W137 -1.66 -0.04 -1.70 -0.38 0.03 -0.35 1.14 -0.02 1.12 -1.04
F169 -1.44 -0.32 -1.76 -0.46 -0.63 -1.09 0.57 0.96 1.52 -1.44
I179 -1.64 -0.13 -1.77 0.13 -0.02 0.11 -0.12 -0.04 -0.16 -1.91
60X L91 -0.75 -0.05 -0.79 0.10 -0.19 -0.09 -0.06 0.11 0.05 -0.90
D93 -0.30 -0.11 -0.41 -15.06 -0.44 -15.5 14.33 0.44 14.77 -1.23
S96 -1.71 -0.35 -2.06 1.34 -0.04 1.30 0.07 0.24 0.31 -0.61
V130 -1.14 -0.11 -1.25 -0.06 0.06 0.00 0.07 0.03 0.10 -1.28
Y132 -2.52 -0.13 -2.65 -0.61 -0.12 -0.73 1.47 0.18 1.65 -2.00
Q134 -0.12 -0.05. -0.17 -0.20 -0.04 -0.24 0.24 0.07 0.31 -0.13
W137 -1.55 -0.04 -1.59 -0.86 0.06 -0.80 1.25 -0.05 1.20 -1.29
F169 -1.04 -0.21 -1.25 -0.10 -0.30 -0.40 0.35 0.47 0.82 -0.90
I179 -1.91 -0.24 -2.15 0.19 -0.11 0.08 -0.11 -0.01 -0.12 -2.31
a All energy components are scaled in kcal/mol. S v d W and B v d W separately indicate contributions of the side chains and backbones to van der Waals interactions ( T v d W ) of inhibitors with residues. S e l e and B e l e respectively correspond to contributions of the side chains and backbones to electrostatic interactions ( T e l e ) of inhibitors with residues. S g b and B g b individually represent contributions of the side chains and backbones to inhibitor-residue polar solvation free energies.
Table 5. Hydrogen bonds formed between inhibitors and residues analyzed using the CPPTRAJ.
Table 5. Hydrogen bonds formed between inhibitors and residues analyzed using the CPPTRAJ.
compound aHydrogen bonds Distance(Å) Angle (°) bOccupancy (%)
60W-BACE1 A93-OD1…60W-H3-N1 3.0 157.1 91.1
A93-OD2…60W-H1-N 2.8 163.3 82.4
A93-OD2…60W-H3-N1 3.2 147.8 54.1
A93-OD1…60W-H1-N 3.1 143.8 46.7
954-BACE1 A93-OD2…954-H4-N2 3.1 150.1 76.7
A93-OD1…954-H4-N2 3.1 149.5 75.6
A93-OD1…954-H2-N1 2.8 159.6 49.2
A93-OD2…954-H2-N1 2.9 158.5 46.9
60X-BACE1 A93-OD2…60X-H5-N4 3.1 150.4 89.1
A93-OD1…60X-H5-N4 3.1 149.3 84.4
A93-OD1…60X-H1-N3 2.8 161.6 62.7
A93-OD2…60X-H1-N3 2.9 157.7 48.7
a Hydrogen bonding interaction are recognized by an acceptor···donor distance of < 3.5 Å and acceptor···H-donor angle of > 120°; b Occupancy (%) is defined as the percentage of simulation time that a specific hydrogen bond exists.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

Subscribe

© 2024 MDPI (Basel, Switzerland) unless otherwise stated