Prediction of Drug Potencies of BACE1 Inhibitors: A Molecular Dynamics Simulation and MM_GB(PB)SA Scoring

: Alzheimer’s disease (AD) is a progressive neurodegenerative brain disorder. One of the important therapeutic approaches of AD is the inhibition of β -site APP cleaving enzyme-1 (BACE1). This enzyme plays a central role in the synthesis of the pathogenic β -amyloid peptides (A β ) in Alzheimer’s disease. A group of potent BACE1 inhibitors with known X-ray structures (PDB ID 5i3X, 5i3Y, 5iE1, 5i3V, 5i3W, 4LC7, 3TPP) were studied by molecular dynamics simulation and binding energy calculation employing MM_GB(PB)SA. The calculated binding energies gave Kd values of 0.139 µ M, 1.39 nM, 4.39 mM, 24.3 nM, 1.39 mM, 29.13 mM, and 193.07 nM, respectively. These inhibitors showed potent inhibitory activities in enzymatic and cell assays. The K d values are compared with experimental values and the structures are discussed in view of the energy contributions to binding. Drug likeness of these inhibitors is also discussed. Accommodation of ligands in the catalytic site of BACE1 is discussed depending on the type of fragment involved in each structure. Molecular dynamics (MD) simulations and energy studies were used to explore the recognition of the selected BACE1 inhibitors by Asp32, Asp228, and the hydrophobic ﬂap. The results show that selective BACE1 inhibition may be due to the formation of strong electrostatic interactions with Asp32 and Asp228 and a large number of hydrogen bonds, in addition to π – π and van der Waals interactions with the amino acid residues located inside the catalytic cavity. Interactions with the ligands show a similar binding mode with BACE1. These results help to rationalize the design of selective BACE1 inhibitors.


Introduction
Alzheimer's disease (AD) is a progressive, neurodegenerative disease of the brain. AD and the associated dementia are connected to amyloid plaque accumulated in the brain. The β-Site Amyloid Precursor Protein cleaving enzyme 1 (BACE1) is an aspartic protease enzyme fixed to the cell membrane; it acts to produce β-amyloid (Aβ) in the signaling pathways in Alzheimer's disease (AD). Excessive accumulation Aβ is believed to induce pathological changes and causes dementia in brains of AD patients.
The enzyme BACE1 initiates the cleavage of amyloid precursor protein (APP) at the β-secretase site, then Aβ is released as a result of further cleavage of the BACE1-cleaved C-terminal APP fragment [1]. Blocking BACE1 proteolytic activity will suppress Aβ generation and reduce the formation of amyloid. Research has been directed towards potent BACE1 inhibitors for AD therapy. Recent breakthroughs in developing BACE1 inhibitors which can penetrate brain cells make the targeting of amyloid deposition-mediated pathology as a therapy more obtainable. Various strategies that have successfully led to the discovery of BACE1 inhibitor drugs have reached the stage of clinical trials. The crucial catalytic aspartate (Asp) dyad, Asp32 and Asp228, is located at the interface of the two lobes [2]. A hairpin loop "flap" in the N terminal lobe partially covers the cleft in a perpendicular orientation and contains Valine 69 Tyrosine 71 and Threonine 72 (colored blue in Figure 1). The conformational changes in the flap control the substrate access to the active site. The first BACE1 substrate analogues inhibitors to mimic the APP-cleavage sequence which contains a non-cleavable peptide bond. This showed high potency but resulted in poor oral bioavailability and low brain penetration, which prevents therapeutic utility [3,4]. Amidine-containing compounds that form optimal interactions with Asp32 and Asp228 enhance the search for BACE1 inhibitors [1]. These Aspbinding amidine and guanidine inhibitors have been studied and the cyclohexyl groups were found to bind the S1 and the lipophilic S1′ pockets (Figure 1). The crucial catalytic aspartate (Asp) dyad, Asp32 and Asp228, is located at the interface of the two lobes [2]. A hairpin loop "flap" in the N terminal lobe partially covers the cleft in a perpendicular orientation and contains Valine 69 Tyrosine 71 and Threonine 72 (colored blue in Figure 1). The conformational changes in the flap control the substrate access to the active site. The first BACE1 substrate analogues inhibitors to mimic the APP-cleavage sequence which contains a non-cleavable peptide bond. This showed high potency but resulted in poor oral bioavailability and low brain penetration, which prevents therapeutic utility [3,4]. Amidine-containing compounds that form optimal interactions with Asp32 and Asp228 enhance the search for BACE1 inhibitors [1]. These Asp-binding amidine and guanidine inhibitors have been studied and the cyclohexyl groups were found to bind the S1 and the lipophilic S1 pockets ( Figure 1). Preparing ligand, receptor, and complex files for Amber [22]: The Antechamber [23] package in Amber Tools [24] was used to create topology and coordinate files for the simulations of ligands. Antechamber is designed to be used with the "general AMBER force field (GAFF)" [23], for organic molecules. This force field has been specifically designed to cover most pharmaceutical (organic) molecules and is compatible with the traditional AMBER force fields in such a way that the two can be mixed during a simulation. Hydrogens were added to the ligand (using reduce) then the ligand. frcmod and library files were prepared for amber, and the tleap editor was used to load the complex or combine the ligand and receptor. The complexes were solvated in a TIP3P [25] cubic water box with water molecules extending 15Å from the complex surface to the water box boundary, and Na + or Clions were added to neutralize the system depending on the charge. The structure of the complex was checked for errors and then converted to topology and coordinate files. The particle mesh Ewald method [26] was used for treating long range electrostatics, and a 9Å cutoff was set for long range interactions. The force field energy of each structure was minimized by progressively relaxing the system before starting the MD simulations. Minimizations were performed employing steepest descent followed by conjugate gradient minimizations (1000 cycles in tandem).
After relaxation of the system it was heated to 300 K by applying harmonic restraint (10 Kcal/Å 2 .mol) on solute. This was followed by an unrestrained 2 ns MD simulation at 300 K and 1 atm to equilibrate the system and adjust the density.
The SHAKE algorithm [27] was used to constrain hydrogen atoms to enable a longer time step (2 fs) in the simulation. A Langevin thermostat [28] with 2 ps −1 collision frequency and weak coupling barostat with 2 ps of relaxation time were employed. Production MD simulations were carried out for 150 ns and resulted in a converged trajectory evident in the RMSD behavior which showed good stability within 1.5Å. Trajectories were collected at 2 ps intervals. These trajectories were used to calculate the binding free energy using MMPBSA.py script [29]; 50 or 100 frames were used in the calculations. Loss in flexibility upon binding expressed as entropy change (T∆S) was calculated by normal modes using the same snapshots that were used for calculating ∆G binding. Then, the absolute free energy of binding was calculated (Equation (5)). The binding energy of the complex was calculated using the MM/GBSA method.

The Generalized Born/Surface Area Model
The MM/PBSA [30] and MM/GBSA methods [31] have been used to estimate ligand-binding affinities in many systems, giving correlation coefficients compared with experiments of R 2 in the range of 0.3 to 0.9, depending on the protein, with MM/GBSA giving better results in this case. The results strongly depend on details in the method, especially the continuum-solvation method, the charges, the dielectric constant, the sampling method, and the entropies. The methods often overestimate differences between sets of ligands.
The (MM/PB(GB)SA [30] method uses representative snapshots from an ensemble of conformations to calculate free energy change between the bound and unbound states of receptor and ligand, Equations (1) and (2)). Before using MM-GBSA [32] the system equilibration was verified by considering temperature, density, total energy, and root mean squared deviation of coordinates (RMSD). An RMSD value relative to the crystal structure of 1.5Å was deemed acceptable. Extensive analysis of each trajectory was performed to make sure the energies calculated are reliable depending on the snapshots [33]. To estimate the total solvation free energy of a molecule, ∆G solv , one typically assumes that it can be decomposed into the "electrostatic" and "non-electrostatic" parts: where ∆G nonel is the free energy of solvating a molecule from which all charges have been removed (i.e., partial charges of every atom are set to zero), and ∆G el is the free energy after removing all charges in vacuum, and then adding them back in the presence of a continuum solvent environment. Generally speaking, ∆G nonel comes from the combined effect of two types of interaction: the favorable van der Waals attraction between the solute and solvent molecules, and the unfavorable cost of breaking the structure of the solvent (water) around the solute. In the current Amber code, this is taken to be proportional to the total solvent accessible surface area (SA) of the molecule, with a proportionality constant derived from experimental solvation energies of small non-polar molecules, and uses a fast LCPO algorithm to compute an analytical approximation to the solvent accessible area of the molecule. Within Amber GB models, each atom in a molecule is represented as a sphere of radius Ri with a charge qi at its center; the interior of the atom is assumed to be filled uniformly with a material of dielectric constant 1. The molecule is surrounded by a solvent of a high dielectric constant 80 for water at 300 K) (Equation (3)). The GB model approximates ∆G el and the nonpolar energy is usually estimated using the solvent-accessible surface area (SASA) (Equation (7)) [8]: ∆G SOL(PB/GB) = ∆G PB/GB + ∆G SA(PB/GB) where ∆E MM is total gas phase energy (sum of ∆E internal + ∆E electrostatic + ∆E vdw ). ∆G SOL(PB/GB) is the sum of nonpolar and polar contributions to solvation calculated by PB or GB. T∆S is conformational entropy upon binding computed by normal-mode analysis on a set of conformational snapshots taken from MD simulations. ∆E internal is internal energy arising from bond, angle, and dihedral terms in the MM force field. ∆E electrostatic is electrostatic energy as calculated by the molecular mechanics (MM) force field. ∆E vdw is the van der Waals contribution from MM. ∆G PB/GB is the nonpolar contribution to the solvation free energy calculated by an empirical model. The nonpolar solvation free energy is typically given by an empirical formula that is proportional to the solvent accessible surface area of the solute: ∆ G SA = γ ·SASA + b, where γ is the surface tension constant and b is a correction constant (γ = 0.00542 kcal·mol −1 ·Å −2 and b = 0.92 kcal/mol in the AMBER package). ∆G SA/GB is the electrostatic contribution to the solvation free energy calculated by the PB or GB method.
One-thousand 2 ps spaced snapshots of each complex were generated from the MD trajectories, and all water molecules and counter-ions were removed before MM-PBSA/GBSA calculations. Coordinates were extracted by using the extract-coords.mmpbsa script and the ∆G values were calculated using the "MMPBSA.py" script [29].

Data Analysis
The RMSDs, dynamic cross-correlation analysis, and principal component analysis (PCA) were processed using the CPPTRAJ module in Amber 18 package [34]. The principal component analysis (PCA) was performed to help in sampling [35,36].
System stability under MD simulations (see Figure 2a-d).
Computation 2020, 8, x FOR PEER REVIEW 6 of 25 constant and b is a correction constant (γ = 0.00542 kcal·mol −1 ·Å −2 and b = 0.92 kcal/mol in the AMBER package). ΔGSA/GB is the electrostatic contribution to the solvation free energy calculated by the PB or GB method. One-thousand 2 ps spaced snapshots of each complex were generated from the MD trajectories, and all water molecules and counter-ions were removed before MM-PBSA/GBSA calculations. Coordinates were extracted by using the extract-coords.mmpbsa script and the ΔG values were calculated using the "MMPBSA.py" script [29].

Data Analysis
The RMSDs, dynamic cross-correlation analysis, and principal component analysis (PCA) were processed using the CPPTRAJ module in Amber 18 package [34]. The principal component analysis (PCA) was performed to help in sampling [35,36].
System stability under MD simulations (see Figure 2a-d). Before starting MD analysis, the root mean square deviation (RMSD) evolution of the protein backbone Cα for each complex was monitored throughout the 200 ns MD simulations to ensure stability of the systems. As shown in Figure 2, the RMSD evolution for Cα of BACE1 bound with inhibitors exhibited relatively small fluctuations at the start of simulation, then was stable and changes were within 1.0 Å. Accordingly, the RMSD evolution of the heavy atoms of the inhibitors maintained relative stability (RMSD fluctuation <1 Å) during the 200 ns simulation. Pairwise RMSD for specific snapshots was computed using pytraj in Amber. The RMSD to the experimental structure reference was computed, then, pairwise RMSD for the first 5000 snapshots, skipping every 10 frames, was computed ( Figure 2d).

Prediction of Binding Mode and Key Interactions of Inhibitors to BACE1
MD simulations were performed to elucidate the key interactions of inhibitors responsible for inhibitory activity against β-amyloid (Aβ) accumulation. The MD simulations were performed to evaluate the favored binding modes and key interactions of BACE1 with various inhibitors ( Figure  3). Before starting MD analysis, the root mean square deviation (RMSD) evolution of the protein backbone Cα for each complex was monitored throughout the 200 ns MD simulations to ensure stability of the systems. As shown in Figure 2, the RMSD evolution for Cα of BACE1 bound with inhibitors exhibited relatively small fluctuations at the start of simulation, then was stable and changes were within 1.0 Å. Accordingly, the RMSD evolution of the heavy atoms of the inhibitors maintained relative stability (RMSD fluctuation <1 Å) during the 200 ns simulation. Pairwise RMSD for specific snapshots was computed using pytraj in Amber. The RMSD to the experimental structure reference was computed, then, pairwise RMSD for the first 5000 snapshots, skipping every 10 frames, was computed ( Figure 2d).

Prediction of Binding Mode and Key Interactions of Inhibitors to BACE1
MD simulations were performed to elucidate the key interactions of inhibitors responsible for inhibitory activity against β-amyloid (Aβ) accumulation. The MD simulations were performed to evaluate the favored binding modes and key interactions of BACE1 with various inhibitors (Figure 3).   The binding energies of inhibitors with BACE1 are shown in Tables 1 and 2 and in Figures 4 and 5;  inhibitors under study bind Asp32 and Asp228 (Tables 3-5, and Figure 6), except for 4LC7 which binds Asp93 and Asp289 (See later). Table 2. The different components of binding free energy (kcal/mol) between the inhibitors-BACE1 complex evaluated using the MM/GBSA method.         (Table 3). Colors: ß sheets Yellow; coils green; flap blue; helices red; Nitrogen atoms blue; inhibitor cyan. All inhibitors occupy similar binding pockets and more importantly form hydrogen bond interactions with the catalytic dyad of Asp32 and Asp228. The active site of BACE1 is mostly hydrophobic, with no charged residues within a distance of 8 Å of the Asp dyad; the aspartate  (Table 3). Colors: ß sheets Yellow; coils green; flap blue; helices red; Nitrogen atoms blue; inhibitor cyan.  All inhibitors occupy similar binding pockets and more importantly form hydrogen bond interactions with the catalytic dyad of Asp32 and Asp228. The active site of BACE1 is mostly hydrophobic, with no charged residues within a distance of 8 Å of the Asp dyad; the aspartate residues form bonds with the amine and the nitrogen of the pyridine ring; see Figure 1 structure 8 and Figures 7 and 8.  The hydrophobic interactions Tyr71, Val69, Gly13, Gly230. Phe108, Leu30, and Ile118 are common in all 68J, 68K, 68L, and 68M inhibitors, and all display hydrophobic contacts with residues. The nitrogen containing heterocycles are often referred to as the aspartyl binding motif; see Figure 8.  Figures 7 and 8, where the terminal CR3 forms hydrophobic interactions in the S2′ pocket which contains D228. The correlation coefficient of binding energy (∆H) for these 4 inhibitors with vdW energy is 0.95 and Esurface is 0.63. The 2aminopyridine fragment forms hydrogen bonds with Asp32 (2.6Å) and a weaker interaction with Asp 228 (4.9 Å). Table 4. Correlation coefficients (R 2 ) of ΔH with contributing energies (from Table 2) for groups of inhibitors. Figure 2)    The hydrophobic interactions Tyr71, Val69, Gly13, Gly230. Phe108, Leu30, and Ile118 are common in all 68J, 68K, 68L, and 68M inhibitors, and all display hydrophobic contacts with residues. The nitrogen containing heterocycles are often referred to as the aspartyl binding motif; see Figure 8.  Figures 7 and 8, where the terminal CR3 forms hydrophobic interactions in the S2′ pocket which contains D228. The correlation coefficient of binding energy (∆H) for these 4 inhibitors with vdW energy is 0.95 and Esurface is 0.63. The 2aminopyridine fragment forms hydrogen bonds with Asp32 (2.6Å) and a weaker interaction with Asp 228 (4.9 Å). Table 4. Correlation coefficients (R 2 ) of ΔH with contributing energies (from Table 2) for groups of inhibitors. Figure 2)   The hydrophobic interactions Tyr71, Val69, Gly13, Gly230. Phe108, Leu30, and Ile118 are common in all 68J, 68K, 68L, and 68M inhibitors, and all display hydrophobic contacts with residues. The nitrogen containing heterocycles are often referred to as the aspartyl binding motif; see Figure 8.

Inhibitors Number (from
Inhibitors 1, 2, 3, and 4 share fragments A and B in Figures 7 and 8, where the terminal CR3 forms hydrophobic interactions in the S 2 pocket which contains D228. The correlation coefficient of binding energy (∆H) for these 4 inhibitors with vdW energy is 0.95 and E surface is 0.63. The 2-aminopyridine fragment forms hydrogen bonds with Asp32 (2.6Å) and a weaker interaction with Asp 228 (4.9 Å).
The correlation with electrostatic energy is very small (Table 4) indicating a mostly hydrophobic interaction on this side. The phenyl rings in structure B (Figure 7) bind Tyr 71 (3.0 Å) and Val69 (4.0 Å). Inhibitors 1 and 2 have an extra fragment C which binds the S3 pocket and differ by one fragment (where fragment D is in inhibitor 1 and replaced by fragment E in inhibitor 2) which binds the S1 pocket and Gly34, where in inhibitor 1 it is D and in inhibitor 2 it is E. Fragment D in Figure 7 with its Лcloud provides stronger interaction than that of E. All differences arise from different vdW interactions, and the π-π stacking interaction between the phenyl-imino group and Phe108 added stability with the enzyme 5i3W-68L (inhibitor 5) binds Asp 32, Asp 228, Gly230, Tyr71, Leu30, and Gly13, see Figure 9E. This inhibitor shares fragment C in Figure 7 with inhibitors 1-4, which binds S1 and yields an experimental ∆G value of −12.5 kcal/mol and comparable vdw energy to other inhibitors 1-4, whereas the calculated value is −8.15(4) kcal/mol. The fragment The correlation with electrostatic energy is very small (Table 4) indicating a mostly hydrophobic interaction on this side. The phenyl rings in structure B (Figure 7) bind Tyr 71 (3.0 Å) and Val69 (4.0 Å). Inhibitors 1 and 2 have an extra fragment C which binds the S3 pocket and differ by one fragment (where fragment D is in inhibitor 1 and replaced by fragment E in inhibitor 2) which binds the S1 pocket and Gly34, where in inhibitor 1 it is D and in inhibitor 2 it is E. Fragment D in Figure 7 with its Л cloud provides stronger interaction than that of E. All differences arise from different vdW interactions, and the π-π stacking interaction between the phenyl-imino group and Phe108 added stability with the enzyme 5i3W-68L (inhibitor 5) binds Asp 32, Asp 228, Gly230, Tyr71, Leu30, and Gly13, See Figure 9E. This inhibitor shares fragment C in Figure 7 with inhibitors 1-4, which binds S1 and yields an experimental ∆G value of −12.5 kcal/mol and comparable vdw energy to other inhibitors 1-4, whereas the calculated value is −8.15(4) kcal/mol. The fragment O N NH 2 binds Asp32 and Leu30. The attachment of the phenyl ring could lead to a significant hydrophobic interaction, which would increase the probability of permeability into the brain. Thus, many BACE1 inhibitors were designed using phenyl -based analogs.
In BACE1 bound to inhibitor 7 (4LC7), shown in Figure 6G, and the heterocyclic pentatomic ring binds Asp93 and Asp289. This feature is shared with 5i3W ( Figure 9E) in which the same ring binds Asp32 and Asp228. Inhibitor 5 (5i3W) has an extra phenyl group that binds the hydrophobic pocket (near Tyr71) which enhanced its binding over 4LC7. Inhibitor 6 (3TPP) has a different structure but shares an aryl ring with other inhibitors and it showed enhanced binding ( Figure 6F Figure 9F).
binds Asp32 and Leu30. The attachment of the phenyl ring could lead to a significant hydrophobic interaction, which would increase the probability of permeability into the brain. Thus, many BACE1 inhibitors were designed using phenyl -based analogs.
In BACE1 bound to inhibitor 7 (4LC7), shown in Figure 6G, and the heterocyclic pentatomic ring binds ld lead to a significant hydrophobic ility into the brain. Thus, many BACE1 G, and the heterocyclic pentatomic ring red with 5i3W ( Figure 9E) in which the an extra phenyl group that binds the g over 4LC7. Inhibitor 6 (3TPP) has a itors and it showed enhanced binding d aryl group interacts with Gln73, the ring-NH binds the other end of Asn233 ther side of Thr231 all make hydrogen 9F).
binds Asp93 and Asp289. This feature is shared with 5i3W ( Figure 9E) in which the same ring binds Asp32 and Asp228. Inhibitor 5 (5i3W) has an extra phenyl group that binds the hydrophobic pocket (near Tyr71) which enhanced its binding over 4LC7. Inhibitor 6 (3TPP) has a different structure but shares an aryl ring with other inhibitors and it showed enhanced binding ( Figure 6F. The sulfate group binds Asn233 and the attached aryl group interacts with Gln73, the fragment The correlation with electrostatic energy is very small (Table 4) indicating a mostly hydrophobic interaction on this side. The phenyl rings in structure B (Figure 7) bind Tyr 71 (3.0 Å) and Val69 (4.0 Å). Inhibitors 1 and 2 have an extra fragment C which binds the S3 pocket and differ by one fragment (where fragment D is in inhibitor 1 and replaced by fragment E in inhibitor 2) which binds the S1 pocket and Gly34, where in inhibitor 1 it is D and in inhibitor 2 it is E. Fragment D in Figure 7 with its Л cloud provides stronger interaction than that of E. All differences arise from different vdW interactions, and the π-π stacking interaction between the phenyl-imino group and Phe108 added stability with the enzyme 5i3W-68L (inhibitor 5) binds Asp 32, Asp 228, Gly230, Tyr71, Leu30, and Gly13, See Figure 9E. This inhibitor shares fragment C in Figure 7 with inhibitors 1-4, which binds S1 and yields an experimental ∆G value of −12.5 kcal/mol and comparable vdw energy to other inhibitors 1-4, whereas the calculated value is −8.15(4) kcal/mol. The fragment O N NH 2 binds Asp32 and Leu30. The attachment of the phenyl ring could lead to a significant hydrophobic interaction, which would increase the probability of permeability into the brain. Thus, many BACE1 inhibitors were designed using phenyl -based analogs.
In BACE1 bound to inhibitor 7 (4LC7), shown in Figure 6G, and the heterocyclic pentatomic ring binds Asp93 and Asp289. This feature is shared with 5i3W ( Figure 9E) in which the same ring binds Asp32 and Asp228. Inhibitor 5 (5i3W) has an extra phenyl group that binds the hydrophobic pocket (near Tyr71) which enhanced its binding over 4LC7. Inhibitor 6 (3TPP) has a different structure but shares an aryl ring with other inhibitors and it showed enhanced binding ( Figure 6F The correlation with electrostatic energy is very small (Table 4) indicating a mostly hydrophobic interaction on this side. The phenyl rings in structure B (Figure 7) bind Tyr 71 (3.0 Å) and Val69 (4.0 Å). Inhibitors 1 and 2 have an extra fragment C which binds the S3 pocket and differ by one fragment (where fragment D is in inhibitor 1 and replaced by fragment E in inhibitor 2) which binds the S1 pocket and Gly34, where in inhibitor 1 it is D and in inhibitor 2 it is E. Fragment D in Figure 7 with its Л cloud provides stronger interaction than that of E. All differences arise from different vdW interactions, and the π-π stacking interaction between the phenyl-imino group and Phe108 added stability with the enzyme 5i3W-68L (inhibitor 5) binds Asp 32, Asp 228, Gly230, Tyr71, Leu30, and Gly13, See Figure 9E. This inhibitor shares fragment C in Figure 7 with inhibitors 1-4, which binds S1 and yields an experimental ∆G value of −12.5 kcal/mol and comparable vdw energy to other inhibitors 1-4, whereas the calculated value is −8.15(4) kcal/mol. The fragment O N NH 2 binds Asp32 and Leu30. The attachment of the phenyl ring could lead to a significant hydrophobic interaction, which would increase the probability of permeability into the brain. Thus, many BACE1 inhibitors were designed using phenyl -based analogs.
In BACE1 bound to inhibitor 7 (4LC7), shown in Figure 6G, and the heterocyclic pentatomic ring binds Asp93 and Asp289. This feature is shared with 5i3W ( Figure 9E) in which the same ring binds Asp32 and Asp228. Inhibitor 5 (5i3W) has an extra phenyl group that binds the hydrophobic pocket (near Tyr71) which enhanced its binding over 4LC7. Inhibitor 6 (3TPP) has a different structure but shares an aryl ring with other inhibitors and it showed enhanced binding ( Figure 6F   The aryl group on the opposite end makes hydrophobic interactions with Phe108, Gly13, Gln12, and Leu30. The oxygen of the peptide bond also interacts with Gln73. The sulfate fragment in 3TPP-5HA binds S2, as seen in the structure below. The aryl group on the opposite end makes hydrophobic interactions with Phe108, Gly13, Gln12, and Leu30. The oxygen of the peptide bond also interacts with Gln73. The sulfate fragment in 3TPP-5HA binds S2, as seen in the structure below

Drug Likeness
Lipinski's rule of five was used to evaluate drug likeness or determine if a compound with a certain pharmacological activity has properties that would make it a feasible orally active drug in humans (Table 6). The rule was based on the observation that most orally administered drugs are relatively small and moderately lipophilic. The rule predicts the absorption, distribution, metabolism, and excretion of the compound. Lipinski's rule states that, in general, an orally active drug has no more than one violation of the following criteria: -No more than five hydrogen bond donors (total H-N, H-O bonds); -No more than 10 hydrogen bond acceptors (all N+O atoms); -Molecular mass less than 500; -LogP value less than 5 (octanol-water partition coefficient); -Drug likeness improved LogP (−0.4 to 5.6), molecular weight 180 to 480, total atoms 20 to 70, including N and O, Veber's Rule: Good oral bioavailability, questioned the 500 molecular weight cutoff. Introduced polar surface area (PSA), no greater than 140 Å 2 , and 10 rotatable bonds or less (Table 7).

Drug Likeness
Lipinski's rule of five was used to evaluate drug likeness or determine if a compound with a certain pharmacological activity has properties that would make it a feasible orally active drug in humans (Table 6). The rule was based on the observation that most orally administered drugs are relatively small and moderately lipophilic. The rule predicts the absorption, distribution, metabolism, and excretion of the compound. Lipinski's rule states that, in general, an orally active drug has no more than one violation of the following criteria: Veber's Rule: Good oral bioavailability, questioned the 500 molecular weight cutoff. Introduced polar surface area (PSA), no greater than 140 Å 2 , and 10 rotatable bonds or less (Table 7).  PSA is a commonly used metric for the optimization of a drug's ability to permeate cells [39]. Molecules with a polar surface area of greater than 140 Å 2 tend to be poor at permeating cell membranes. For molecules to penetrate the blood-brain barrier; a PSA less than 90 Å 2 is usually needed [39]. Inspecting the properties of the seven inhibitors used (Table 6) shows that each can be a suitable drug.
Inhibitors 1, 2, 3, and 4, which share a hydrophobic mioty ( Figure 7) and the 2-aminopyridine fragment (Figure 7) in their structure, showed the best correlation between PSA (Table 8) and binding energy (∆H), E vdw , E surface , and E electrostatic ; Table 4. Furthermore, the vdw energy showed the best correlation with PSA for these inhibitors. Inhibitors 1 to 6 showed the best correlation with surface area energy (Figure 7). When structure 5 was added to the group, the correlation of PSA with E electrostatic improved due to the presence of hydrogen bond donors and acceptors in inhibitor 5, but the correlation with E surface was not changed. Analysis of energies involved in the binding of inhibitors to BACE1 will aid the design of new inhibitors with more efficacy. Ligand efficiency (LE) [40] is calculated by scaling affinity by molecular size (Table 9). LE was introduced as a metric for the molecular structure to normalize the affinity with respect to molecular size by scaling the standard free energy of binding (∆G • ) with the number of heavy atoms (N nH ), using the formula: LE (T, P, C) = −∆G/N nH Table 8. The areas of hydrophobic pockets in BACE1 for each inhibitor binding ( Figure S5). The calculated energies resulting from hydrophobicity using the formula −25 cal/Å 2 of surface area and comparing the estimated hydrophobic energy with that resulting from reported polar surface area (PSA) [21].  LE values vary with conditions, and a value of 0.3 or higher is considered favorable. LE decreases with an increase in the number of heavy atoms. There was no obvious trend followed in the inhibitors in this work due to variation in structure. This variation in the results of high energy cost for desolvation of ligands depends on the changes that took place. Ligand efficiency values of inhibitors were in the range of 0.09 to 0.41 (Table 9).

Proteins in
The drug-like properties when applying Lipinski's rule of five, Veber's Rule, and the MDDR Rule (No. rings ≥ 3, No. rigid bonds ≥ 18, No. rotatable bonds ≥ 6) changed depending on functional groups and molecular weights. There was a good correlation between the Gibbs free energy (∆G) calculated and the experimentally obtained values [21,41].

Conclusions
The parameters for successful drugs depend on the specificity and binding to the receptor; a 500 molecuar weight and a Kd value in the range of 1 to 10 nM is preferred for good absorption. The potency depends on the specificity of binding (Asp) and increased hydrophobic binding residues are preferred. However, this depends on the specificity, and a balance between specific binding and hydrphobicity should be maintained. The higher LE, the more promising the drug binding to a specific target.
The binding energy of drug to its target depends on a group of energies [42]; the first is the desolvation energy needed to break the hydrogen bonds between the inhibitor and solvent, then the energy released upon binding of the inhibitor to the receptor, and burying the inhibitor hydrophobic surface. Polar interactions and hydrophobic surface burial depend on the surface area (every 1Å 2 of surface area releases approximately 25 cal); see Tables 7 and 8. The disadvantage in the drugs under study is the limited surface area, of around 90 Å 2 , for the drugs to enter brain cells. Differences between calculated and actual ∆G values are due to imperfect H-bonds caused by steric and distance factors, which result in a higher E-cost for desolvation.
Research on the mechanism of AD has considered BACE1 as a key enzyme that participates in the formation of Aβ, and which broadly exists in the brains of AD patients. Compounds with peptide mimetic structures are effective in BACE1 inhibition according to experimental aspartic proteinase results in vitro. Nevertheless, these kinds of BACE1 inhibitors did not perform well in pre-clinical trials due to their excessive number of hydrogen bond donors and acceptors, which increase the polarity and further lead to a lack of permeability across the BBB (Blood Brain Barrier). Based on molecular dynamics and energy studies, the amino acid residues Asp228 and Asp32 in the BACE1 enzyme play an important role in the interactions between compounds and the enzyme. Furthermore, S1, S3, S2 and other pockets also exhibited a central role in binding with the BACE1 inhibitors. In the light of these studies, compounds with amino heterocycles were designed and synthesized. The presence of amino and aromatic rings maintained the inhibitory ability while decreasing the polarity of the structure.
MM/PBSA energies are calculated for snapshots obtained by MD simulations. Variations are normally solved by calculating only interaction energies, studying many snapshots, and using several independent simulations. It has been suggested that the calculations can be performed using only minimized structures, but such results may depend on the starting structures. Finally, MM/GBSA, compared with other ligand binding methods, showed reasonable accuracy.
MM/GBSA is a popular method for calculating absolute binding affinities with a modest computational effort. Energy results from six well-defined terms. However, questions remain about the dielectric constant, parameters for the non-polar energy, the radii used for the PB or GB calculations, whether to include the entropy term, and whether to perform MD simulations or minimizations. In practice, MM/GBSA typically provides results of intermediate quality, which are often better than docking and scoring but worse than FEP; for example, R 2 = 0.3 for the whole PDB bind database, but R 2 = 0.0-0.8 for individual proteins.
Funding: This Research was funded by Birzeit University.