Allostery Inhibition of BACE1 by Psychotic and Meroterpenoid Drugs in Alzheimer’s Disease Therapy

In over a century since its discovery, Alzheimer’s disease (AD) has continued to be a global health concern due to its incurable nature and overwhelming increase among older people. In this paper, we give an overview of the efforts of researchers towards identifying potent BACE1 exosite-binding antibodies and allosteric inhibitors. Herein, we apply computer-aided drug design (CADD) methods to unravel the interactions of some proposed psychotic and meroterpenoid BACE1 allosteric site inhibitors. This study is aimed at validating the allosteric potentials of these selected compounds targeted at BACE1 inhibition. Molecular docking, molecular dynamic (MD) simulations, and post-MD analyses are carried out on these selected compounds, which have been experimentally proven to exhibit allosteric inhibition on BACE1. The SwissDock software enabled us to identify more than five druggable pockets on the BACE1 structural surface using docking. Besides the active site region, a melatonin derivative (compound 1) previously proposed as a BACE1 allostery inhibitor showed appreciable stability at eight different subsites on BACE1. Refinement with molecular dynamic (MD) simulations shows that the identified non-catalytic sites are potential allostery sites for compound 1. The allostery and binding mechanism of the selected potent inhibitors show that the smaller the molecule, the easier the attachment to several enzyme regions. This finding hereby establishes that most of these selected compounds failed to exhibit strong allosteric binding with BACE1 except for compound 1. We hereby suggest that further studies and additional identification/validation of other BACE1 allosteric compounds be done. Furthermore, this additional allosteric site investigation will help in reducing the associated challenges with designing BACE1 inhibitors while exploring the opportunities in the design of allosteric BACE1 inhibitors.

The early symptoms of AD include disruption of daily life activities due to memory loss; difficulty in understanding visual images and completing known tasks; challenges with identification, speaking, and writing new words; disorganization; and constant mood swings [12,13]. Available AD treatments have focused on decreasing the level of Aβ aggregation and the accumulation of hyperphosphorylated tau. An alternative approach is using drugs to reverse the symptoms, prevent tissue/cell damage, reduce fats/swollenness, inoculate, and for hormonal treatment [14]. Currently, some tau-related therapies are under clinical trial. For an efficient AD treatment approach, early treatment assessment strategies and diagnostic biochemical markers are crucial [15]. The current research on AD treatment is three-fold. This includes identifying populations with higher vulnerability (persons > 60 years) and helping them with the preventions intended to modify the risk factors. Furthermore, taking advantage of the modern neuroimaging and cerebrospinal fluid examinations and subsequent data analysis in early diagnosis of the persons at  [1,11]. The cleavage mechanism by γ-secretase is non-amyloidogenic, while β-secretase or α-secretase cleavage results in the amyloidogenic pathway. CTFα and CTFβ are alpha and beta carboxy-terminal fragments, AICD represents the amyloid precursor protein intracellular domain.
The early symptoms of AD include disruption of daily life activities due to memory loss; difficulty in understanding visual images and completing known tasks; challenges with identification, speaking, and writing new words; disorganization; and constant mood swings [12,13]. Available AD treatments have focused on decreasing the level of Aβ aggregation and the accumulation of hyperphosphorylated tau. An alternative approach is using drugs to reverse the symptoms, prevent tissue/cell damage, reduce fats/swollenness, inoculate, and for hormonal treatment [14]. Currently, some tau-related therapies are under clinical trial. For an efficient AD treatment approach, early treatment assessment strategies and diagnostic biochemical markers are crucial [15]. The current research on AD treatment is three-fold. This includes identifying populations with higher vulnerability (persons > 60 years) and helping them with the preventions intended to modify the risk factors. Furthermore, taking advantage of the modern neuroimaging and cerebrospinal fluid examinations and subsequent data analysis in early diagnosis of the persons at the preclinical stage. And finally, the identification of disease-modifying compounds to inhibit Aβ and the NFTs accumulation [16].
Since its discovery in 1999, BACE1 has had a primary focus on the reduction of β-amyloid plaques in the brain. Although there is no approved BACE1 inhibitor [17], evidence of its inhibition towards AD management remains feasible. Its inhibition is one of the therapies targeted toward AD management [18,19]. Researchers have made relentless efforts to design potential BACE1 inhibitors with desirable properties in the reduction of Aβ accumulation in the brain. Six potential BACE1 drugs such as Atabecestat, Umibecestat, LY3202626, Elenbecestat, Lanabecestat, and Verubecestat, showed initial benefits in reducing the amyloid-beta plaques responsible for AD-related symptoms [18][19][20]. Several complications from these BACE1 drugs at phase II/III clinical trials informed discontinuation [21,22]. There is a need to understand why these BACE1 drugs failed. Good knowledge of the inhibiting potentials and properties of BACE1 inhibitors serves as a prerequisite for the development of molecules with high selectivity and specificity [18][19][20]. Also, more attention on BACE1 non-active sites inhibition might be a promising approach in targeting them in AD management.
Molecules 2022, 27, x FOR PEER REVIEW 3 of 18 the preclinical stage. And finally, the identification of disease-modifying compounds to inhibit Aβ and the NFTs accumulation [16]. Since its discovery in 1999, BACE1 has had a primary focus on the reduction of βamyloid plaques in the brain. Although there is no approved BACE1 inhibitor [17], evidence of its inhibition towards AD management remains feasible. Its inhibition is one of the therapies targeted toward AD management [18,19]. Researchers have made relentless efforts to design potential BACE1 inhibitors with desirable properties in the reduction of Aβ accumulation in the brain. Six potential BACE1 drugs such as Atabecestat, Umibecestat, LY3202626, Elenbecestat, Lanabecestat, and Verubecestat, showed initial benefits in reducing the amyloid-beta plaques responsible for AD-related symptoms [18][19][20]. Several complications from these BACE1 drugs at phase II/III clinical trials informed discontinuation [21,22]. There is a need to understand why these BACE1 drugs failed. Good knowledge of the inhibiting potentials and properties of BACE1 inhibitors serves as a prerequisite for the development of molecules with high selectivity and specificity [18][19][20]. Also, more attention on BACE1 non-active sites inhibition might be a promising approach in targeting them in AD management.
Allosteric inhibition is a form of noncompetitive inhibition that occurs at another protein site different from the active site [26]. In allosteric inhibition, the inhibitor binds to a distant site away from the active site, rendering the active site unfit for the substrate to bind [26,27]. We lately reviewed different studies from researchers on BACE1 allosteric sites inhibition and exosite-binding antibody development [20]. The survey facilitates the documentation of some inhibitors with potential propensities to the BACE1 allosteric sites ( Figure 3). Most of these compounds identified using experimental approaches are yet to be explored at the molecular level to evaluate if they indeed interact at the allosteric regions of BACE1. However, some phytocompounds, including 3,5,7,3 ,4 -pentamethoxyflavone (PMF) [28], nor-rubrofusarin 6-O-β-D-glucoside [29], loganin [30], and gamma-linolenic acid [31] have been proposed through docking and kinetics study as allosteric inhibitors. Besides, Rombouts et al. [32] applied a fragment-based virtual screening approach to identify BACE1 inhibitors, which bind at the other subsites without interacting with Asp32 and Asp228. Integrated analysis, including nuclear magnetic resonance (NMR), fluorescence resonance energy transfer (FRET) assay, and ThermoFluor (TF), produced six hits [32]. Refinement and analysis showed that four of these compounds [32] competitively bind with OM99-2 [33]. X-ray atomistic interaction revealed that one of these hits occupied the S1 and S3 subsites without interacting with the catalytic Asp32 and Asp228 ( Figure 3). Compound 12 with protein data bank (PDB) [34] code 5MXD [32] showed an IC 50 value of 0.5 mM, which is significant. Although the interaction pose occurred close to the BACE1 active site region, we assume this mechanism as a potential allosteric inhibition. Allosteric inhibition is a form of noncompetitive inhibition that occurs at another protein site different from the active site [26]. In allosteric inhibition, the inhibitor binds to a distant site away from the active site, rendering the active site unfit for the substrate to bind [26,27]. We lately reviewed different studies from researchers on BACE1 allosteric sites inhibition and exosite-binding antibody development [20]. The survey facilitates the documentation of some inhibitors with potential propensities to the BACE1 allosteric sites ( Figure 3). Most of these compounds identified using experimental approaches are yet to be explored at the molecular level to evaluate if they indeed interact at the allosteric regions of BACE1. However, some phytocompounds, including 3,5,7,3′,4′-pentamethoxyflavone (PMF) [28], nor-rubrofusarin 6-O-β-D-glucoside [29], loganin [30], and gamma-linolenic acid [31] have been proposed through docking and kinetics study as allosteric inhibitors. Besides, Rombouts et al. [32] applied a fragment-based virtual screening approach to identify BACE1 inhibitors, which bind at the other subsites without interacting with Asp32 and Asp228. Integrated analysis, including nuclear magnetic resonance (NMR), fluorescence resonance energy transfer (FRET) assay, and ThermoFluor (TF), produced six hits [32]. Refinement and analysis showed that four of these compounds [32] competitively bind with OM99-2 [33]. X-ray atomistic interaction revealed that one of these hits occupied the S1 and S3 subsites without interacting with the catalytic Asp32 and Asp228 ( Figure 3). Compound 12 with protein data bank (PDB) [34] code 5MXD [32] showed an IC50 value of 0.5 mM, which is significant. Although the interaction pose occurred close to the BACE1 active site region, we assume this mechanism as a potential allosteric inhibition. Figure 3. Zoomed-in 3D snapshot of BACE1 crystal structure complexed with compound 12 (PDB code: 5MXD) devoid of Asp32/228 interaction [32]. The compound binds at the flap region of BACE1 and we generated this image with the Discovery Studio R2017 [35].
Herein, we apply molecular modeling methods, including molecular docking and molecular dynamics (MD) simulations to unravel the interaction of some potential allosteric inhibitors with BACE1. Readers may be puzzled by the uniqueness of the present work since there exist investigations on BACE1 allosterism. Indeed, there are propositions on BACE1 allosteric or secondary sites inhibition, there is still a knowledge gap at the molecular level using computational methods, like MD simulations. Besides, Pietro et al.'s (2017) study of conformational ensemble and binding mode analysis of some multisite inhibitors using MD and docking method [36], studies that explore multisite targeting of BACE1 towards drug design for AD are still scarce. In this study, we selected compounds ( Figure 4) proposed by different authors [28][29][30][31][37][38][39][40][41] as potent BACE1 allosteric inhibitors for some computational experiments. Docking analysis revealed that these compounds bind at allosteric sites but are highly stable at the active site. Our observation aligns with the kinetic experiments and suggestions available in the literature [28][29][30][31][37][38][39][40][41] for the various compounds. Only compound 1 (a melatonin derivative) with 88% C Figure 3. Zoomed-in 3D snapshot of BACE1 crystal structure complexed with compound 12 (PDB code: 5MXD) devoid of Asp32/228 interaction [32]. The compound binds at the flap region of BACE1 and we generated this image with the Discovery Studio R2017 [35].
Herein, we apply molecular modeling methods, including molecular docking and molecular dynamics (MD) simulations to unravel the interaction of some potential allosteric inhibitors with BACE1. Readers may be puzzled by the uniqueness of the present work since there exist investigations on BACE1 allosterism. Indeed, there are propositions on BACE1 allosteric or secondary sites inhibition, there is still a knowledge gap at the molecular level using computational methods, like MD simulations. Besides, Pietro et al.'s (2017) study of conformational ensemble and binding mode analysis of some multisite inhibitors using MD and docking method [36], studies that explore multisite targeting of BACE1 towards drug design for AD are still scarce. In this study, we selected compounds ( Figure 4) proposed by different authors [28][29][30][31][37][38][39][40][41] as potent BACE1 allosteric inhibitors for some computational experiments. Docking analysis revealed that these compounds bind at allosteric sites but are highly stable at the active site. Our observation aligns with the kinetic experiments and suggestions available in the literature [28][29][30][31][37][38][39][40][41] for the various compounds. Only compound 1 (a melatonin derivative) with 88% inhibitory activity at 5 µM BACE1 concentration [37] showed appreciable stability at six different allostery pockets. We further refined the stability of the representative molecule (compound 1) at the BACE1 binding region to establish if these sites are transient or long-lived using several metrics. Understanding how potent small inhibitors modulate and inhibit BACE1 at the molecular level would enable us to manipulate synthetic enzymes or design drugs for AD treatment.
Molecules 2022, 27, x FOR PEER REVIEW 5 of 18 inhibitory activity at 5 μM BACE1 concentration [37] showed appreciable stability at six different allostery pockets. We further refined the stability of the representative molecule (compound 1) at the BACE1 binding region to establish if these sites are transient or longlived using several metrics. Understanding how potent small inhibitors modulate and inhibit BACE1 at the molecular level would enable us to manipulate synthetic enzymes or design drugs for AD treatment.

System Preparation
BACE1 X-ray 3D structures are available as several complexes in the PDB repositories. We selected 6PZ4 [42] to maintain consistency with previous studies [1,43]. Investigations have shown favorable outcomes with BACE1 mono-protonation of Asp32 [1,44], hence, we protonated Asp32 using PROPKA [45] at pH 7. Protein structure refinement entails using the Maestro protein preparation wizard package [46] to add hydrogens, assign bond orders, remove water and non-ligand molecules. We added missing residues in MODELLER 9.19 [47] implemented in the UCSF Chimera [48]. The ligands were prepared for docking in GaussView 6.0.16 [49], optimized in Avogadro [50], and saved in Mol2 format for the docking study.

Docking
Docking is a molecular modeling method that involves predicting the preferred orientation of one molecule when bound to another molecule [51]. It requires computational software compatible to dock small compounds into a macromolecule including protein structures [52]. The conformations of the docked ligand-enzyme complex prediction involve assessing the different poses of the ligand within the receptor's binding site. The scoring of the various poses enables predicting the molecular mechanisms, the binding free energy of the complex, and nonbonded interaction properties [53]. Over the years, the combination of docking and further binding free energy predictions remained useful in designing protease inhibitors such as HIV and BACE1 proteases [54,55].
We used SwissDock online docking software [56,57] to search for the regions where the selected compounds bind to the BACE1. SwissDock automatically prepares the ligand and enzyme structures before docking using a webserver. The implemented docking algorithm is the CHARMM force field [58] to offer an accurate docking approach. Therefore, ligands and the protein files (in Mol2 and PDB) automatically convert to CHARMM format after upload. Although prediction of ligand poses are premised on experimentally identified binding sites, the SwissDock is advantageous because it enables flexible

System Preparation
BACE1 X-ray 3D structures are available as several complexes in the PDB repositories. We selected 6PZ4 [42] to maintain consistency with previous studies [1,43]. Investigations have shown favorable outcomes with BACE1 mono-protonation of Asp32 [1,44], hence, we protonated Asp32 using PROPKA [45] at pH 7. Protein structure refinement entails using the Maestro protein preparation wizard package [46] to add hydrogens, assign bond orders, remove water and non-ligand molecules. We added missing residues in MODELLER 9.19 [47] implemented in the UCSF Chimera [48]. The ligands were prepared for docking in GaussView 6.0.16 [49], optimized in Avogadro [50], and saved in Mol2 format for the docking study.

Docking
Docking is a molecular modeling method that involves predicting the preferred orientation of one molecule when bound to another molecule [51]. It requires computational software compatible to dock small compounds into a macromolecule including protein structures [52]. The conformations of the docked ligand-enzyme complex prediction involve assessing the different poses of the ligand within the receptor's binding site. The scoring of the various poses enables predicting the molecular mechanisms, the binding free energy of the complex, and nonbonded interaction properties [53]. Over the years, the combination of docking and further binding free energy predictions remained useful in designing protease inhibitors such as HIV and BACE1 proteases [54,55].
We used SwissDock online docking software [56,57] to search for the regions where the selected compounds bind to the BACE1. SwissDock automatically prepares the ligand and enzyme structures before docking using a webserver. The implemented docking algorithm is the CHARMM force field [58] to offer an accurate docking approach. Therefore, ligands and the protein files (in Mol2 and PDB) automatically convert to CHARMM format after upload. Although prediction of ligand poses are premised on experimentally identified binding sites, the SwissDock is advantageous because it enables flexible docking [56,57]. This unconstrained docking facilitates several ligand conformational orientations within the enzyme. The automated process also reduces human errors by using the web interface in generating alternative input files and parameters while interpreting the docking outcome [56]. SwissDock webserver is incorporated with the Apache web server, PHP (opensource technology) [59], which has dual Xeon E440 2.83 GHz at 1.7 Å and 16 GB. The docking procedures include minimizing at 100 steps steepest descent to remove any clashes inherent from adjusting the protein structure. Also, it uses a 5 kcal/mol.Å −2 constraint to restrict nonphysical movements of heavy atoms during minimization.

Molecular Dynamic Simulations
MD MD simulation is an approach in computational modeling to explore the conformational dynamics of molecules [60]. Molecular simulation provides a similitude of real-time possibilities to predict protein or molecule behavior with or without interacting with other molecules [61]. MD simulation is crucial in drug design; it facilitates binding and catalytic mechanism prediction after identifying a molecule with a plausible propensity for a target [62,63]. Many researchers have indicated MD as crucial to refine, validate, and improve docking evaluation [61,64]. There are several commercial, non-commercial, and web-based MD simulation programs for molecular modeling studies [63] in which the AMBER suite [65] is a prominent one.
We performed MD simulations on the AMBER 18 program package integrated with graphics processing units (GPU) [66] using the particle mesh Ewald method (PMEMD) [67] package with the Sander module. The long-range electrostatic interactions cut-off is at 12 Å. The simulation pre-step involves generating a partial atomic charge for the ligand using the General Amber Force Field (GAFF) [68] of the ANTECHAMBER module. The GAFF is a simple harmonic function developed as a complete and suitable force field for rational drug design [69]. Further procedures are topology and coordinate preparation, enzyme-ligand coupling, solvation, and neutralization. Explicit solvation was done with the Leap module of the AMBER 18 package in a TIP3P orthorhombic water box at 10 Å to any edge. The pre-MD production steps are partial minimization, total minimization, heating, and equilibration.
The partial minimization was at 10,000 steepest descent steps with 10 kcal/mol.Å −2 harmonic restraint on all heavy atoms to relax the system and remove potential atom clashes. We also run a full minimization with another steepest and conjugate gradient descents at 5000 steps each without constraints. Systems heating was from 0 K to 300 K for 300 ps using Langevin dynamics, 1 ps collision frequency, and a 5 kcal/mol.Å −2 applied harmonic restraints at a constant volume. We subsequently equilibrated for 500 ps at 300 K under constant pressure and temperature (NPT) ensemble. The final step is MD simulation for 120 ns at 2 fs time step without any restraints at 300 K and 1 atm in the NPT ensemble switching on the Langevin temperature scaling [70] and Berendsen barostat [71] algorithm for temperature and pressure, respectively. The molecular dynamic simulation also involves applying the SHAKE algorithm [72] to constrain hydrogen atoms.

Post-Simulation Analysis
We measured system stability through root-mean-square deviation (RMSD) calculation. The RMSD trajectory of the protein backbone alpha carbon (Cα) generated with the CPPTRAJ [73] module uses Equation (1) for its estimation. The standard deviation of the interatomic distance between Cα backbone atoms of two amino acids v and w at n points in Equation (1) represents v i as Cα coordinates in v at the time i, and w i is the coordinates of Cα atom in w at the time i.
The radius of gyration (RoG) is the moment inertia of atoms from their center of mass. RoG is often applicable in quantifying the molecular rigidity of a system [74]. The RoG  (2)) is the square root of the inertia moment of atoms, where n is the number of atoms, r i represents the atomic position, and r m signifies the mean position of all atoms.
We estimate the root-mean-square fluctuation (RMSF) to predict the conformational changes on a per-residue basis. Equation (3) shows the RMSF equation whereby x i(j) represents the i-th Cα atom position in the j-th model structure, and (x i ) denotes the averaged location of the i-th Cα backbone atom in all models.
Further analyses include secondary structure prediction with the definition secondary structure of protein (DSSP) approach of Kabsch and Sander [75] implemented in the CPPTRAJ program. The DSSP approach involves analyzing the most likely secondary structure assignment through the 3D structure of a protein. It entails interpreting atom position in a protein and calculating the hydrogen bond (HB) energy between all atoms. DSSP algorithm ignores any hydrogen available in the input structure then computes the significant hydrogen positions after placing them at 1.0 Å from the backbone nitrogen in the reverse direction from the backbone carbonyl bond [75]. The assignment completes by using the top two HBs in N and C=O to predict the most likely secondary structure class for each residue in the protein.

Allosteric Sites and Prediction
Blind molecular docking with SwissDock generated different binding poses at various sites on the BACE1 protein structure. All the compounds (Figure 4) show favorable interaction with BACE1 widely around the active site region and sparsely at other domains ( Figure 5). This observation indicates that these compounds are selective for the BACE1 active site over other subsites. We notice a significant overlay of LY2811376 poses on the BACE1 flap tip as the only identified allostery site. The observation is consistent with previous reports [38,39] on LY2811376 modulatory and inhibitory potentials. All the predicted sites are feasible subsites as suggested in a survey of ligand-binding modes of co-crystallized BACE1 structures [76] and molecular modeling study [36].
Zoomed-in profile of all compounds' interactions at the BACE1 active site is available in the supporting information (SI) captioned as Figure S1. Gamma-linolenic and sargahydroquinoic acids with the most favored binding affinity (Table 1) at BACE1 active site show two distinct hydrogen bond (HB) interactions ( Figure S1) via Thr33/Gly34 and Gly34/Gln73 (Table 2), respectively. Electrostatic interaction from hydrogen bonding is imperative for inhibitor stability at the active site of an enzyme. The electrostatic effect represents an approach to predict catalysis [77]. Compound 1 shows an interaction at eight different subsites on the BACE1 scaffold ( Figure 6), thereby indicating the molecule as a potent allostery modulator and multitarget directed ligand (MTDL). The compound is relatively sizeable compared to others (Figure 4) with unique functional groups, including amide. These physicochemical properties make compound 1 flexible with the propensity to attach at several regions on the BACE1 scaffold through different molecular interactions, including electrostatic and van der Waals ( Figure 6). Molecules 2022, 27, x FOR PEER REVIEW 8 of 18 Zoomed-in profile of all compounds' interactions at the BACE1 active site is available in the supporting information (SI) captioned as Figure S1. Gamma-linolenic and sargahydroquinoic acids with the most favored binding affinity (Table 1) at BACE1 active site show two distinct hydrogen bond (HB) interactions ( Figure S1) via Thr33/Gly34 and Gly34/Gln73 (Table 2), respectively. Electrostatic interaction from hydrogen bonding is imperative for inhibitor stability at the active site of an enzyme. The electrostatic effect represents an approach to predict catalysis [77]. Compound 1 shows an interaction at eight different subsites on the BACE1 scaffold ( Figure 6), thereby indicating the molecule as a potent allostery modulator and multitarget directed ligand (MTDL). The compound is relatively sizeable compared to others (Figure 4) with unique functional groups, including amide. These physicochemical properties make compound 1 flexible with the propensity to attach at several regions on the BACE1 scaffold through different molecular interactions, including electrostatic and van der Waals ( Figure 6).

LY2811376
Gamma-linolenic acid Sargahydroquinoic acid Anisoperidone Compound 1 Figure 5. Binding modes of the allostery inhibitors on BACE1 through docking. LY2811376 converges mainly at the active site and partially on top of the flap and gamma-linolenic acid binds favorably at the primary site and two other secondary sites as predicted previously [31]. Sargahydroquinoic acid and anisoperidone interact at the active site and extend to the flap tip, while compound 1 binds at several regions.      Table 1 shows the full fitness energy and the predicted binding scores for all the potent molecules' interactions at the active site. The binding poses scored using their full fitness and clustered show interaction energy values within −7.51 and −8.64 kcal/mol. Although with the lowest binding score, compound 1 has the most favored full fitness of its moieties in the BACE1 active site. Its binding energy to other binding regions within the enzyme ranges from −6.12 to 7.10 kcal/mol with favorable fitness scores ( Figure 6). The interaction of Thr232 and flap region residues buttress previous molecular docking prediction of compound 1 binding pose [37]. To give a clearer picture of the stability of compound 1 interaction to BACE1 at the different binding subsites, we simulate each complex ( Figure 6) using the amber force field. The molecular docking approach enables us to unravel druggable regions for potential multisite targeting in BACE1 towards AD treatment.

Protein Structural Changes
We assess changes to the BACE1 backbone structure when bound to compound 1 using RMSD calculation of the protein Cα atoms. In molecular modeling, RMSD scoring enables conformational dynamics and system stability prediction. Figure 7 depicts the RMSD plots per time in which allosteric sites 2 and 3 ( Figure 6) are unavailable because their MD simulation failed. The data indicate the relative stability of the protein's primary structure over the 120 nanoseconds production run. Each complex shows appreciable overlap with relative stability around 90-120 ns, indicating that all the protein structures converged without large-scale conformational transitions. RMS deviation higher than 3.5 Å is a potential indication of significant BACE1 conformational switch in domains such as flap opening and closing [43]. Evolution of the alpha carbon backbone atom root-mean-square deviation (RMSD/angstrom) over 120 ns molecular simulations for compound 1 binding to various regions on BACE1 structure. Allos. site represents the allosteric site and the definition for each site including the interaction residues is available in Figure 6.
Compound 1 binds to the active site with an average RMSD of 2.162 ± 0.206 Å, while allosteric sites 1, 4-8 show mean RMSD scores within 1.766 and 2.072 Å (Table 3). This result denotes that compound 1 binding at the active site increases the protein backbone dynamics compared to the six potential allosteric sites. Note that the minimum RMSD value is 0 Å, and allostery site 4 shows the lowest outcome in all the parameters (Table 3), while the free BACE1 structure has an average RMSD value of 3.75 ± 0.44 Å [43].  Evolution of the alpha carbon backbone atom root-mean-square deviation (RMSD/angstrom) over 120 ns molecular simulations for compound 1 binding to various regions on BACE1 structure. Allos. site represents the allosteric site and the definition for each site including the interaction residues is available in Figure 6.
Compound 1 binds to the active site with an average RMSD of 2.162 ± 0.206 Å, while allosteric sites 1, 4-8 show mean RMSD scores within 1.766 and 2.072 Å (Table 3). This result denotes that compound 1 binding at the active site increases the protein backbone dynamics compared to the six potential allosteric sites. Note that the minimum RMSD value is 0 Å, and allostery site 4 shows the lowest outcome in all the parameters (Table 3), while the free BACE1 structure has an average RMSD value of 3.75 ± 0.44 Å [43]. To estimate the effect of compound 1 interaction on BACE1 intrinsic arrangement during the simulation, we estimate the protein secondary structure using the DSSP [75] method. Figure 8 shows that the ligand binding at the various spaces induces small changes to the protein structure. For instance, the binding of compound 1 at allosteric site 1 diminishes the alpha arrangement to anti around residues 335-338. BACE1 has a few highly folded alpha helixes (pink color, Figure 5). Alteration of the helix to anti (Figure 8) might impact the protein's structural stability significantly. The DSSP map of allosteric sites 4-8 (SI, Figure S2) is like the active site ( Figure 8). changes to the protein structure. For instance, the binding of compound 1 at allosteric site 1 diminishes the alpha arrangement to anti around residues 335-338. BACE1 has a few highly folded alpha helixes (pink color, Figure 5). Alteration of the helix to anti (Figure 8) might impact the protein's structural stability significantly. The DSSP map of allosteric sites 4-8 (SI, Figure S2) is like the active site ( Figure 8).

Allostery Dynamics
We use the RMSF metric to probe how compound 1 dynamically perturbs the motion of the entire protein structure. The computation includes all atomic flexibility (both backbone and sidechain atoms) to identify the effect of the studied molecule at the various regions on BACE1 per residue. RMSF is a crucial concept to evaluate protein dynamics and individual residue flexibility [78]. The RMSF projections for compound 1 binding (Figure 9) are fascinating, with a high spectrum separation between the active site and others (potential allosteric sites). This outcome signifies that compound 1 binding to allostery sites decreases BACE1 flexibility, whereas residue mobility is high when bound to the active site. The fluctuation pattern and projection for active site binding (Figure 9) are akin to our previous simulation of apo BACE1 RMSF [43], denoting that compound 1 binding at the active region is most likely inconsequential on the protein mobility, thus supporting Panyatip et al. [37] in silico predictions. Eighteen residues show very high fluctuations (>20 Å) with the highest RMSF of 23.422 Å (see maximum value in Table 4) in Ser58 for the active site system.

Allostery Dynamics
We use the RMSF metric to probe how compound 1 dynamically perturbs the motion of the entire protein structure. The computation includes all atomic flexibility (both backbone and sidechain atoms) to identify the effect of the studied molecule at the various regions on BACE1 per residue. RMSF is a crucial concept to evaluate protein dynamics and individual residue flexibility [78]. The RMSF projections for compound 1 binding (Figure 9) are fascinating, with a high spectrum separation between the active site and others (potential allosteric sites). This outcome signifies that compound 1 binding to allostery sites decreases BACE1 flexibility, whereas residue mobility is high when bound to the active site. The fluctuation pattern and projection for active site binding (Figure 9) are akin to our previous simulation of apo BACE1 RMSF [43], denoting that compound 1 binding at the active region is most likely inconsequential on the protein mobility, thus supporting Panyatip et al. [37] in silico predictions. Eighteen residues show very high fluctuations (>20 Å) with the highest RMSF of 23.422 Å (see maximum value in Table 4) in Ser58 for the active site system.
Molecules 2022, 27, x FOR PEER REVIEW 13 of 18 Figure 9. Per residue (a total of 385) fluctuation through all-atom root-mean-square fluctuation (RMSF/angstrom) scoring for compound 1 binding to several regions on BACE1. Allos. site represents the allosteric site and the definition for each site including the interacting residues is available in Figure 6. Table 4. Statistics of the RMSF (in Å) scores per residue. Std. represents standard and Allos. site represents the allosteric site and the definition for each site including the interaction residues is available in Figure 6. We propose that the allostery mechanism of compound 1 on BACE1 is both modulatory and inhibitory. The analysis reflects that the melatonin derivative can allosterically restrict BACE1 motility. It acts as an inhibitor at the potential allosteric sites and as a modulator at the active site. This phenomenon indicates that compound 1 is a potential MTDL targeting multiple subsites on BACE1. The mechanism would distort BACE1 structural arrangement and availability for natural substrate binding, lowering its level in AD conditions. Irrespective of the targeted subsite, compound 1 exerts relatively low flexibility across all the 385 residues. The average RMSF values are within 1.095 and 1.348 Å (Table  4), with allosteric site 7 showing slightly higher projections (above 5 Å) at residues 164-167.

Rigid Body Motion
The mass-weighted radius of gyration (RoG) is a moment of atom inertia from their center of mass [79]. The calculation setup for RoG analysis includes all the atoms to notice Figure 9. Per residue (a total of 385) fluctuation through all-atom root-mean-square fluctuation (RMSF/angstrom) scoring for compound 1 binding to several regions on BACE1. Allos. site represents the allosteric site and the definition for each site including the interacting residues is available in Figure 6. Table 4. Statistics of the RMSF (in Å) scores per residue. Std. represents standard and Allos. site represents the allosteric site and the definition for each site including the interaction residues is available in Figure 6. We propose that the allostery mechanism of compound 1 on BACE1 is both modulatory and inhibitory. The analysis reflects that the melatonin derivative can allosterically restrict BACE1 motility. It acts as an inhibitor at the potential allosteric sites and as a modulator at the active site. This phenomenon indicates that compound 1 is a potential MTDL targeting multiple subsites on BACE1. The mechanism would distort BACE1 structural arrangement and availability for natural substrate binding, lowering its level in AD conditions. Irrespective of the targeted subsite, compound 1 exerts relatively low flexibility across all the 385 residues. The average RMSF values are within 1.095 and 1.348 Å (Table 4), with allosteric site 7 showing slightly higher projections (above 5 Å) at residues 164-167.

Rigid Body Motion
The mass-weighted radius of gyration (RoG) is a moment of atom inertia from their center of mass [79]. The calculation setup for RoG analysis includes all the atoms to notice the effect of allostery binding on the BACE1 surface. Figure 10 shows the RoG evolution per time for each model system. Besides sites 7 and 8 with average RoG of 21.5 and 21.4 Å, the approximate average RoG value is 21.2 Å (Table 5). These values are close, denoting that the atoms in the BACE1 structure have similar moment inertia in the various models. However, allosteric site 7 facilitates slightly higher RoG ( Figure 10) like its RMSF (Figure 9). The data signifies less compactness in the allosteric site 7 BACE1 model.
Also, allostery site 7 shows the highest data in all the parameters, while site 5 shows the least mean RoG score ( Table 5). The RoG metric enables us to predict how BACE1 moves as a rigid body when complexed with the melatonin derivative at the active region and allosteric sites. The outcome shows that compound 1 exerts compactness on the BACE1, restricting protein motions during the simulation. The similarity of the RoG projections and the same mean values depicts compound 1 as a potential active site and allosteric inhibitor. This outcome denotes that the computationally predicted sites are potential target sites for drug binding in BACE1. the effect of allostery binding on the BACE1 surface. Figure 10 shows the RoG evolution per time for each model system. Besides sites 7 and 8 with average RoG of 21.5 and 21.4 Å, the approximate average RoG value is 21.2 Å (Table 5). These values are close, denoting that the atoms in the BACE1 structure have similar moment inertia in the various models. However, allosteric site 7 facilitates slightly higher RoG ( Figure 10) like its RMSF ( Figure  9). The data signifies less compactness in the allosteric site 7 BACE1 model.   Figure 10. Time (in nanoseconds) evolution of the mass-weighted RoG-radius of gyration (in angstrom) induced by coupling compound 1 to the various regions on BACE1 structure. Allos. site represents the allosteric site and the definition for each site including the interacting residues is available in Figure 6.

Conclusions
Allosteric inhibition is a type of noncompetitive inhibition occurring at another protein site different from the active site. In this mechanism, the inhibitor binds to a site(s) other than the active site, thereby rendering the active site unfit for the substrate to bind. In allosteric inhibition, it is a case in which the moiety gets to the enzyme or protein first, which blocks the other from binding. As a sequel to our previous review study on some identified allosteric inhibitors and exosite-binding antibodies of BACE1 over the last 8 years (2013-2020), we, hereby, in this study applied CADD methods to further validate the claims as contained in the reviewed study. We selected six compounds that were reported to bind on sites different from the BACE1 active site. Some of these compounds are classified as meroterpenoids (Sargahydroquinoic acid and Gamma-Linolenic acid) and psychotic drugs (Anisoperidone). Good knowledge of the inhibiting potentials and properties of BACE1 inhibitors also serves as a prerequisite for the development of molecules with high selectivity and specificity. We show in this study that besides the active site, BACE1 targeting at other subsites is plausible towards AD therapy. From the results of the molecular docking, molecular dynamic (MD) simulations, and subsequent post-MD analyses, the identified regions correspond with some predicted subsites and the small molecules (like compound 1) bind through several other BACE1 residues. The allostery and binding mechanisms of the selected potent inhibitors show that the smaller the molecule, the easier the attachment to several enzyme regions. This finding hereby establishes that most of these selected compounds failed to exhibit strong allosteric binding with BACE1 except for compound 1. We hereby suggest further studies and additional identification/validation of other BACE1 allosteric compounds be done. Furthermore, this additional allosteric site investigation will help in reducing the associated challenges with designing BACE1 inhibitors, while exploring the opportunities in the design of allosteric BACE1 inhibitors. We further suggest that diverting attention from the conventional active site through a detailed computer-aided drug design approach would assist in designing a potential BACE1 inhibitor that might attain the approval stage for AD management.
Supplementary Materials: The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/molecules27144372/s1, Figure S1: The interaction profile of compounds at the BACE1 active site and Figure S2: The protein secondary structure plot for compound 1 interaction at the various allosteric sites of BACE1.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.