A PM7 dynamic residue-ligand interactions energy landscape of the BACE1 inhibitory pathway by hydroxyethylamine compounds. Part I: The flap closure process

BACE1 is an enzyme of scientific interest because it participates in the progression of Alzheimer’s disease. Hydroxyethylamines (HEAs) are a family of compounds which exhibit inhibitory activity toward BACE1 at a nanomolar level, favorable pharmacokinetic properties and oral bioavailability. The first step in the inhibition of BACE1 by HEAs consists of their entrance into the protease active site and the resultant conformational change in the protein, from Apo to closed form. These two conformations differ in the position of an antiparallel loop (called the flap) which covers the entrance to the catalytic site. For BACE1, closure of this flap is vital to its catalytic activity and to inhibition of the enzyme due to the new interactions thereby formed with the ligand. In the present study a dynamic energy landscape of residue-ligand interaction energies (ReLIE) measured for 112 amino acids in the BACE1 active site and its immediate vicinity during the closure of the flap induced by 8 HEAs of different inhibitory power is presented. A total of 6.272 million ReLIE calculations, based on the PM7 semiempirical method, provided a deep and quantitative view of the first step in the inhibition of the aspartyl protease. The information suggests that residues Asp93, Asp289, Thr292, Thr293, Asn294 and Arg296 are anchor points for the ligand, accounting for approximately 45% of the total protein-ligand interaction. Additionally, flap closure improved the BACE1-HEA interaction by around 25%. Furthermore, the inhibitory activity of HEAs could be related to the capacity of these ligands to form said anchor point interactions and maintain them over time: the lack of some of these anchor interactions delayed flap closure or impeded it completely, or even caused the flap to reopen. The methodology employed here could be used as a tool to evaluate future structural modifications which lead to improvements in the favorability and stability of BACE1-HEA ReLIEs, aiding in the design of better inhibitors.     

Determination of the protonation state for the catalytic dyad in β-secretase when bound to hydroxyethylamine transition state analogue inhibitors: A molecular dynamics simulation study

BACE1 is an aspartyl protease of pharmacological interest for its direct participation in Alzheimer’s disease (AD) through β-amyloid peptide production. Two aspartic acid residues are present in the BACE1 catalytic region which can adopt multiple protonation states depending on the chemical nature of its inhibitors, i.e., monoprotonated, diprotonated and di-deprotonated states. In the present study a series of protein-ligand molecular dynamics (MD) simulations was carried out to identify the most feasible protonation state adopted by the catalytic dyad in the presence of hydroxyethylamine transition state analogue inhibitors. The MD trajectories revealed that the di-deprotonated state is most prefered in the presence of hydroxyethilamine (HEA) family inhibitors. This appears as a result after evaluating, for all 9 protonation state configurations during the simulation time, the deviations of a set of distances and dihedral angles measured on the ligand, protein and protein-ligand complex with reference to an X-ray experimental BACE1/HEA crystallographic structure. These results will help to clarify the phenomena related to the HEAs inhibitory pathway, and improve HEAs databases’ virtual screening and ligand design processes targeting β-secretase protein.     

Modeling the protonation states of β-secretase binding pocket by molecular dynamics simulations and docking studies

β-secretase (BACE1) is an aspartyl protease that processes the β-amyloid peptide in the human brain in patients with Alzheimer’s disease. There are two catalytic aspartates (ASP32 and ASP228) in the active domain of BACE1. Although it is believed that the net charge of the Asp dyad is −1, the exact protonation state still remains a matter of debate. We carried out molecular dynamic (MD) simulations for the four protonation states of BACE1 proteins. We applied Glide docking studies to 21 BACE1 inhibitors against the MD extracted conformations. The dynamic results infer that the protein/ligand complex remains stable during the entire simulation course for HD32D228 model. The results show that the hydrogen bonds between the inhibitor and the Asp dyad are maintained in the 10,000th ps snapshot of HD32D228 model. Our results also reveal the significant loop residues in maintaining the active binding conformation in the HD32D228 model. Molecular docking results show that the HD32D228 model provided the best enrichment factor score, suggesting that this model was able to recognize the most active compounds. Our observations provide an evidence for the preference of the anionic state (HD32D228) in BACE1 binding site and are in accord with reported computational data. The protonation state study would provide significant information to assign the correct protonation state for structure-based drug design and docking studies targeting the BACE1 proteins as a tactic to develop potential AD inhibitors.     

Computational insights into the protonation states of catalytic dyad in BACE1–acyl guanidine based inhibitor complex

Developing small compound based drugs targeting the β-secretase (BACE) enzyme is one of the most promising strategies in treatment of the Alzheimer’s disease. As the enzyme shows the activity based on the acid-base reaction at a very narrow pH range, the protonation state of aspartic acids with the residue number 32 and 228 (Asp32 and Asp228), which forms the active site dyad, along with the protonation state of the ligand (substrate or inhibitor) play very critical role in interactions between the ligand and enzyme. Thus, understanding the nature of the protonation state of both enzyme’s active site dyad and ligand is crucial for drug design in Alzheimer’s disease field. Here we have investigated the protonation state of the Asp32 and Asp228 residues in the presence of a highly potent beta secretase inhibitor, containing acyl guanidine warhead that have recently been devised but not extensively studied. Our Quantum Mechanical, Molecular Dynamics and Docking studies on all the possible protonation states have suggested that the dyad residues are in di-deprotonated states in the presence of protonated inhibitor.     

Structural analysis of lead fullerene-based inhibitor bound to human immunodeficiency virus type 1 protease in solution from molecular dynamics simulations

Molecular dynamics (MD) simulations of the HIV-1 protease (HIVP) complexed with lead fullerene-based inhibitor (diphenyl C60 alcohol) in the three protonated states, unprotonated (Un-), monoprotonated (Mono-), and diprotonated (Di-) states at Asp25 and Asp25' were performed. As the X-ray structure of the investigated complex is not available, it was built up starting with the X-ray crystallographic structure of the HIVP complexed with non-peptide inhibitor (PDB code: 1AID) and that of the diphenyl C60 alcohol optimized using the integrated ONIOM molecular orbital calculations. The inhibitor was, then, introduced into the enzyme pocket using a molecular docking method. Change of the HIVP binding cavity for all three states were evaluated in terms of distance between the two catalytic residues, Asp25 and Asp25' as well as those between the catalytic residues and the flap regions. The torsional angles formed by the O-C-C-O of the two carboxyl groups of the catalytic dyad show the non-planar configuration with the most frequency at about -45 degrees for the Un-, 35 degrees and -95 degrees for the Mono- and 60 degrees for the Di-systems. At equilibrium, different orientations of the fullerene-based inhibitor in the three protonation states were observed. For the Di-state, the OH group of the inhibitor stably forms hydrogen bonds with the two aspartic residues. It turns to the flap region to form hydrogen bonding to the backbone N of Ile50' for the Un-state. In contrast, the OH group turns to locate between the catalytic and the flap region for the Mono-states. Beside the molecular orientation, the rotation of the OH group of the inhibitor in the Un-state was also detected. In terms of solvation, the carboxylate oxygens of the aspartic residues in the Un- and Mono-states were solvated by one to three water molecules while the OH group in these two states was coordinated by one water molecule. This is in contrast to the Di-state in which no water molecule is available in the radius of 5-6A around the oxygen atoms of the carboxylate groups of enzyme and of the OH group of the inhibitor. The simulated results lead to the conclusion that the active site of the HIVP complexed with the diphenyl C60 alcohol is the diprotonation states on Asp25 and Asp25'.     

Structural and dynamical properties of different protonated states of mutant HIV-1 protease complexed with the saquinavir inhibitor studied by molecular dynamics simulations

To understand the basis of drug resistance, particularly of the HIV-1 PR, three molecular dynamics (MD) simulations of HIV-1 PR mutant species, G48V, complexed with saquinavir (SQV) in explicit aqueous solution with three protonation states, diprotonation on Asp25 and Asp25' (Di-pro) and monoprotonation on each Asp residue (Mono-25 and Mono-25'). For all three states, H-bonds between saquinavir and HIV-1 PR were formed only in the two regions, flap and active site. It was found that conformation of P2 subsite of SQV in the Mono-25 state differs substantially from the other two states. The rotation about 177 degrees from the optimal structure of the wild type was observed, the hydrogen bond between P2 and the flap residue (Val48) was broken and indirect hydrogen bonds with the three residues (Asp29, Gly27, and Asp30) were found instead. In terms of complexation energies, interaction energy of -37.3 kcal/mol for the Mono-25 state is significantly lower than those of -30.7 and -10.7kcal/mol for the Mono-25' and Di-pro states, respectively. It was found also that protonation at the Asp25 leads to a better arrangement in the catalytic dyad, i.e., the Asp25-Asp25' interaction energy of -8.8 kcal/mol of the Mono-25 is significantly lower than that of -2.6kcal/mol for the Mono-25' state. The above data suggest us to conclude that interaction in the catalytic area should be used as criteria to enhance capability in drug designing and drug screening instead of using the total inhibitor/enzyme interaction.     

Flap opening mechanism of HIV-1 protease

The active site of aspartic proteases, such as HIV-1 protease (PR), is covered by one or more flaps, which restrict access to the active site. For HIV-1 PR, X-ray diffraction studies suggested that in the free enzyme the two flaps are packed onto each other loosely in a semi-open conformation, while molecular dynamics (MD) studies observed that the flaps can also separate into open conformations. In this study, the mechanism of flap opening and the structure and dynamics of HIV-1 PR with semi-open and open flap conformations were investigated using molecular dynamics simulations. The flaps showed complex dynamic behavior as two distinct mechanisms of flap opening and various stable flap conformations (semi-open, open and curled) were observed during the simulations. A network of weakly polar interactions between the flaps were proposed to be responsible for stabilizing the semi-open flap conformation. It is hypothesized that such interactions could be responsible for making flap opening a highly sensitive gating mechanism which control access to the active site.     

Computational design of V-CoCrFeMnNi high-entropy alloys: An atomistic simulation study

In the high-entropy alloy (HEA) community, many researchers have been trying to improve the strength of the CoCrFeMnNi HEA by generating a transformation-induced-plasticity (TRIP) effect and/or maximizing the solid solution hardening effect. Adding vanadium (V) to the CoCrFeMnNi HEAs could be an effective way to improve strength, because vanadium stabilizes the body-centered cubic (bcc) phase and its atomic size is larger than Co, Cr, Fe, Mn, and Ni. To design high strength V-added HEAs, we investigated the effect of vanadium on the critical resolved shear stress (CRSS) by utilizing an atomistic simulation, proposing an empirical equation to estimate the relative effect of alloying elements on the CRSS. For this, we first developed the Co-Cr-Fe-Mn-Ni-V hexanary interatomic potential by newly developing the Cr-V, Fe-V, and Mn-V binary interatomic potentials. As a result, two novel V-added HEAs were designed and the designed HEAs show higher strength than the previously developed non-equiatomic CoCrFeMnNi HEAs, as predicted from the empirical equation.     

A QM/MM study of the catalytic mechanism of aspartate ammonia lyase

Aspartate ammonia lyase (Asp) is one of three types of ammonia lyases specific for aspartate or its derivatives as substrates, which catalyzes the reversible reaction of l-aspartate to yield fumarate and ammonia. In this paper, the catalytic mechanism of Asp has been studied by using combined quantum-mechanical/molecular-mechanical (QM/MM) approach. The calculation results indicate that the overall reaction only contains two elementary steps. The first step is the abstraction of C_β proton of l-aspartate by Ser318, which is calculated to be rate limiting. The second step is the cleavage of C_αN bond of l-aspartate to form fumarate and ammonia. Ser318 functions as the catalytic base, whereas His188 is a dispensable residue, but its protonation state can influence the active site structure and the existing form of leaving amino group, thereby influences the activity of the enzyme, which can well explain the pH dependence of enzymatic activity. Mutation of His188 to Ala only changes the active site structure and slightly elongates the distance of C_β proton of substrate with Ser318, causing the enzyme to remain significant but reduced activity.     

Molecular dynamics simulations of mutated Mycobacterium tuberculosis l-alanine dehydrogenase to illuminate the role of key residues

l-Alanine dehydrogenase from Mycobacterium tuberculosis (l-MtAlaDH) catalyzes the NADH-dependent interconversion of l-alanine and pyruvate, and it is considered to be a potential target for the treatment of tuberculosis. The experiment has verified that amino acid replacement of the conserved active-site residues which have strong stability and no great changes in biological evolutionary process, such as His96 and Asp270, could lead to inactive mutants [Ågren et al., J. Mol. Biol. 377 (2008) 1161–1173]. However, the role of these conserved residues in catalytic reaction still remains unclear. Based on the crystal structures, a series of mutant structures were constructed to investigate the role of the conserved residues in enzymatic reaction by using molecular dynamics simulations. The results show that whatever the conserved residues were mutated, the protein can still convert its conformation from open state to closed state as long as NADH is present in active site. Asp270 maintains the stability of nicotinamide ring and ribose of NADH through hydrogen bond interactions, and His96 is helpful to convert the protein conformation by interactions with Gln271, whereas, they would lead to the structural rearrangement in active site and lose the catalytic activity when they were mutated. Additionally, we deduce that Met301 plays a major role in catalytic reaction due to fixing the nicotinamide ring of NADH to prevent its rotation, and we propose that Met301 would be mutated to the hydrophobic residue with large steric hindrance in side chain to test the activity of the protein in future experiment.     

Getting it right: modeling of pH, solvent and "nearly" everything else in virtual screening of biological targets

“Getting it right” refers to the careful modeling of all elements in the living system, i.e. biological macromolecules, ligands and water molecules. In addition, careful attention should be paid to the protonation state of ionizable functional groups on the ligands and residues at the active site. Computational technology based on the empirical HINT program is described to: (1) calculate free energy scores for ligand binding; (2) include the implicit and explicit effects of water in and around the ligand binding site; and (3) incorporate the effects of global and local pH in molecular models. This last point argues for the simultaneous consideration of a number of molecular models, each with different protonation profiles. Data from recent studies of protein–ligand systems (trypsin, thrombin, neuraminidase, HIV-1 protease and others) are used to illustrate the concepts in the paper. Also discussed are experimental factors related to accurate free energy predictions with this and other computational technologies.     

Identification of potential bivalent inhibitors from natural compounds for acetylcholinesterase through in silico screening using multiple pharmacophores

The symptomatic cure observed in the treatment of Alzheimer's disease (AD) by FDA approved drugs could possibly be due to their specificity against the active site of acetylcholinesterase (AChE) and not by targeting its pathogenicity. The AD pathogenicity involved in AChE protein is mainly due to amyloid beta peptide aggregation, which is triggered specifically by peripheral anionic site (PAS) of AChE. In the present study, a workflow has been developed for the identification and prioritization of potential compounds that could interact not only with the catalytic site but also with the PAS of AChE. To elucidate the essential structural elements of such inhibitors, pharmacophore models were constructed using PHASE, based on a set of fifteen best known AChE inhibitors. All these models on validation were further restricted to the best seven. These were transferred to PHASE database screening platform for screening 89,425 molecules deposited at the “ZINC natural product database”. Novel lead molecules retrieved were subsequently subjected to molecular docking and ADME profiling. A set of 12 compounds were identified with high pharmacophore fit values and good predicted biological activity scores. These compounds not only showed higher affinity for catalytic residues, but also for Trp86 and Trp286, which are important, at PAS of AChE. The knowledge gained from this study, could lead to the discovery of potential AChE inhibitors that are highly specific for AD treatment as they are bivalent lead molecules endowed with dual binding ability for both catalytic site and PAS of AChE.     

A study of hybrid evolutionary algorithms for single machine scheduling problem with sequence-dependent setup times

We present a systematic comparison of hybrid evolutionary algorithms (HEAs), which independently use six combinations of three crossover operators and two population updating strategies, for solving the single machine scheduling problem with sequence-dependent setup times. Experiments show the competitive performance of the combination of the linear order crossover operator and the similarity-and-quality based population updating strategy. Applying the selected HEA to solve 120 public benchmark instances of the single machine scheduling problem with sequence-dependent setup times to minimize the total weighted tardiness widely used in the literature, we achieve highly competitive results compared with the exact algorithm and other state-of-the-art metaheuristic algorithms in the literature. Meanwhile, we apply the selected HEA in its original form to deal with the unweighted 64 public benchmark instances. Our HEA is able to improve the previous best known results for one instance and match the optimal or the best known results for the remaining 63 instances in a reasonable time.     

Molecular modeling study on the effect of residues distant from the nucleotide-binding portion on RNA binding in Staphylococcus aureus Hfq

Hfq is an abundant RNA-binding bacterial protein that was first identified in E. coli as a required host factor for phage Qβ RNA replication. The pleiotrophic phenotype resulting from the deletion of Hfq predicates the importance of this protein. Two RNA-binding sites have been characterized: the proximal site which binds sRNA and mRNA and the distal site which binds poly(A) tails. Previous studies mainly focused on the key residues in the proximal site of the protein. A recent mutation study in E. coli Hfq showed that a distal residue Val43 is important for the protein function. Interestingly, when we analyzed the sequence and structure of Staphylococcus aureus Hfq using the CONSEQ server, the results elicited that more functional residues were located far from the nucleotide-binding portion (NBP). From the analysis seven individual residues Asp9, Leu12, Glu13, Lys16, Gln31, Gly34 and Asp40 were selected to investigate the conformational changes in Hfq–RNA complex due to point mutation effect of those residues using molecular dynamics simulations. Results showed a significant effect on Asn28 which is an already known highly conserved functionally important residue. Mutants D9A, E13A and K16A depicted effects on base stacking along with increase in RNA pore diameter, which is required for the threading of RNA through the pore for the post-translational modification. Further, the result of protein stability analysis by the CUPSAT server showed destabilizing effect in the most mutants. From this study we characterized a series of important residues located far from the NBP and provide some clues that those residues may affect sRNA binding in Hfq.     

Structure-based design of inhibitors of NS3 serine protease of hepatitis C virus

We have designed small focused combinatorial library of hexapeptide inhibitors of NS3 serine protease of the hepatitis C virus (HCV) by structure-based molecular design complemented by combinatorial optimisation of the individual residues. Rational residue substitutions were guided by the structure and properties of the binding pockets of the enzyme's active site. The inhibitors were derived from peptides known to inhibit the NS3 serine protease by using unusual amino acids and alpha-ketocysteine or difluoroaminobutyric acid, which are known to bind to the S1 pocket of the catalytic site. Inhibition constants (Ki) of the designed library of inhibitors were predicted from a QSAR model that correlated experimental Ki of known peptidic inhibitors of NS3 with the enthalpies of enzyme-inhibitor interaction computed via molecular mechanics and the solvent effect contribution to the binding affinity derived from the continuum model of solvation. The library of the optimised inhibitors contains promising drug candidates-water-soluble anionic hexapeptides with predicted Ki* in the picomolar range.     

Linear interaction energy models for β-secretase (BACE) inhibitors: Role of van der Waals,electrostatic, and continuum-solvation terms

Computing the binding affinity of a protein–ligand complex is one of the most fundamental and difficult tasks in computer-aided drug design. Many approaches for computing binding affinities can be classified as linear interaction energy (LIE) models as they rely on some type of linear fit of computed interaction energies between ligand and protein. We have examined the computed interaction energies of a series of β-secretase (BACE) inhibitors in terms of van der Waals, coulombic, and continuum-solvation contributions to ligand binding. We have also systematically examined the effect of different protonation states of the protein and ligands. We find that the binding affinities are relatively insensitive to the protonation state of the protein when neutral ligands are considered. Inclusion of charged ligands leads to large deviations in the coulomb, solvation, and even van der Waals terms. The latter is due to increased repulsive van der Waals interactions in the complex due to the strong coulomb attraction found between oppositely charged functional groups in the protein and ligand. In general, we find that the best models are obtained when the protein is judiciously charged (e.g. Asp32⁻, Arg235⁺) and the potentially charged ligands are treated as neutral.     

Iron(II) triggered conformational changes in Escherichia coli fur upon DNA binding: a study using molecular modeling

In order to identify the Fur dimerization domain, a three-dimensional structure of the ferric uptake regulation protein from Escherichia coli (Fur EC) was determined using homology modeling and energy minimization. The Fur monomer consists of turn- helix -turn motif on the N-terminal domain, followed by another helix-turn-helix-turn motif, and two beta-strands separated by a turn which forms the wing. The C-terminal domain, separated by a long coil from the N-terminal, and consisting of two anti parallel beta strands, and a turn-helix-turn-helix-turn motif. Residues in central domain were found to aid the dimer formation, residues 45-70 as evident in the calculated distances; this region is rich in hydrophobic residues. Most interactions occur between residues Val55, Leu53, Gln52, Glu49 and Tyr56 with closest contacts occurring at residues 49-56. These residues are part of an alpha-helix (alpha(4)) near the N-terminal. Upon raising the Fe(2+) concentration the binding of Fur dimer to DNA was enhanced, this was evident when, the Fur EC dimer was docked onto DNA "iron box" (it was found to bind the AT-rich region) and upon addition of Fe(2+) the helices near the N-terminal bound to the major groove of the DNA. Addition of high Fe(2+) concentration triggered further conformational changes in the Fur dimer as was measured by distances between the two subunits, Fe(2+) mediated the Fur binding to DNA by attaching itself to the DNA. At the same time DNA changed conformation as was evident in the distortion in the backbone and the shrinking of major groove distance from 11.4 to 9.3A. Two major Fe(2+) sites were observed on the C-terminal domain: site 1, the traditional Zn site, the cavity contains the residues Cys92, Cys95, Asp137, Asp141, Arg139, Glu 140, His 145 and His 143 at distances range from 1.3 to 2.2A. Site 2 enclave consists of His71, Ile50, Asn72, Gly97, Asp105 and Ala109 at very close proximity to Fe(2+). The closest contacts between Fur dimer and DNA at the AT-rich region were at residues Ala11, Gly12, Leu13, Pro18 and Arg19 mostly hydrophobic residues near the N-terminal domain. Close contacts repeated at His87, His88 and Arg112, and a third region near the C-terminal at Asn137, Arg 139, Glu140, Asn141, His143, Asn141 and His145. Fur dimer has three major contact regions with DNA, the first on the N-terminal domain, a second smaller region at His87, His88 and Arg112 mediated by Fe(2+) ions, and a third region on the C-terminal domain consisting mainly of hydrophobic contacts and mediated by Fe(2+) ions at high concentration.     

Linear interaction energy models for beta-secretase (BACE) inhibitors: Role of van der Waals, electrostatic, and continuum-solvation terms

Computing the binding affinity of a protein–ligand complex is one of the most fundamental and difficult tasks in computer-aided drug design. Many approaches for computing binding affinities can be classified as linear interaction energy (LIE) models as they rely on some type of linear fit of computed interaction energies between ligand and protein. We have examined the computed interaction energies of a series of β-secretase (BACE) inhibitors in terms of van der Waals, coulombic, and continuum-solvation contributions to ligand binding. We have also systematically examined the effect of different protonation states of the protein and ligands. We find that the binding affinities are relatively insensitive to the protonation state of the protein when neutral ligands are considered. Inclusion of charged ligands leads to large deviations in the coulomb, solvation, and even van der Waals terms. The latter is due to increased repulsive van der Waals interactions in the complex due to the strong coulomb attraction found between oppositely charged functional groups in the protein and ligand. In general, we find that the best models are obtained when the protein is judiciously charged (e.g. Asp32⁻, Arg235⁺) and the potentially charged ligands are treated as neutral.