Molecular dynamics of HIV-1 protease

Characterizing Protein-Ligand Binding Using Atomistic Simulation and Machine Learning: Application to Drug Resistance in HIV-1 Protease

Over the past several decades, atomistic simulations of biomolecules, whether carried out using molecular dynamics or Monte Carlo techniques, have provided detailed insights into their function. Comparing the results of such simulations for a few closely related systems has guided our understanding of the mechanisms by which changes such as ligand binding or mutation can alter the function. The general problem of detecting and interpreting such mechanisms from simulations of many related systems, however, remains a challenge. This problem is addressed here by applying supervised and unsupervised machine learning techniques to a variety of thermodynamic observables extracted from molecular dynamics simulations of different systems. As an important test case, these methods are applied to understand the evasion by human immunodeficiency virus type-1 (HIV-1) protease of darunavir, a potent inhibitor to which resistance can develop via the simultaneous mutation of multiple amino acids. Complex mutational patterns have been observed among resistant strains, presenting a challenge to developing a mechanistic picture of resistance in the protease. In order to dissect these patterns and gain mechanistic insight into the role of specific mutations, molecular dynamics simulations were carried out on a collection of HIV-1 protease variants, chosen to include highly resistant strains and susceptible controls, in complex with darunavir. Using a machine learning approach that takes advantage of the hierarchical nature in the relationships among the sequence, structure, and function, an integrative analysis of these trajectories reveals key details of the resistance mechanism, including changes in the protein structure, hydrogen bonding, and protein-ligand contacts.

NMR and MD studies combined to elucidate inhibitor and water interactions of HIV-1 protease and their modulations with resistance mutations

Over the last two decades, both the sensitivity of NMR and the time scale of molecular dynamics (MD) simulation have increased tremendously and have advanced the field of protein dynamics. HIV-1 protease has been extensively studied using these two methods, and has presented a framework for cross-evaluation of structural ensembles and internal dynamics by integrating the two methods. Here, we review studies from our laboratories over the last several years, to understand the mechanistic basis of protease drug-resistance mutations and inhibitor responses, using NMR and crystal structure-based parallel MD simulations. Our studies demonstrate that NMR relaxation experiments, together with crystal structures and MD simulations, significantly contributed to the current understanding of structural/dynamic changes due to HIV-1 protease drug resistance mutations.

Learning the Relationship between the Primary Structure of HIV Envelope Glycoproteins and Neutralization Activity of Particular Antibodies by Using Artificial Neural Networks

The dependency between the primary structure of HIV envelope glycoproteins (ENV) and the neutralization data for given antibodies is very complicated and depends on a large number of factors, such as the binding affinity of a given antibody for a given ENV protein, and the intrinsic infection kinetics of the viral strain. This paper presents a first approach to learning these dependencies using an artificial feedforward neural network which is trained to learn from experimental data. The results presented here demonstrate that the trained neural network is able to generalize on new viral strains and to predict reliable values of neutralizing activities of given antibodies against HIV-1.

Human Cytomegalovirus Protease: Why is the Dimer Required for Catalytic Activity?

Human cytomegalovirus (HCMV) is a pathogenic agent responsible for morbidity and mortality in immunocompromised and immunosuppressed individuals. HCMV encodes a serine protease that is essential for the production of infectious virions. In this work, we applied molecular dynamics (MD) simulations on HCMV protease models in order to investigate the experimentally observed (i) catalytic activity of the enzyme homodimer and (ii) induced-fit mechanism upon the binding of substrates and peptidyl inhibitors. Long and stable trajectories were obtained for models of the monomeric and dimeric states, free in solution and bound covalently and noncovalently to a peptidyl-activated carbonyl inhibitor, with very good agreement between theoretical and experimental results. The MD results suggest that HCMV protease indeed operates by an induced-fit mechanism. Also, our analysis indicates that the catalytic activity of the dimer is a result of more favorable interactions between the oxyanion in the covalently bound state and the backbone nitrogen of Arg165, resulting in a reaction that is 7.0 kcal/mol more exergonic and a more significant thermodynamic driving force. The incipient oxyanion in the transition state should also benefit from the stronger interactions with Arg165, reducing in this manner the intrinsic activation barrier for the reaction in the dimeric state.

Effects of PRE and POST therapy drug-pressure selected mutations on HIV-1 protease conformational sampling

Conformational sampling of pre- and post-therapy subtype B HIV-1 protease sequences derived from a pediatric subject infected via maternal transmission with HIV-1 were characterized by double electron-electron resonance spectroscopy. The conformational ensemble of the PRE construct resembles native-like inhibitor bound states. In contrast, the POST construct, which contains accumulated drug-pressure selected mutations, has a predominantly semi-open conformational ensemble, with increased populations of open-like states. The single point mutant L63P, which is contained in PRE and POST, has decreased dynamics, particularly in the flap region, and also displays a closed-like conformation of inhibitor-bound states. These findings support our hypothesis that secondary mutations accumulate in HIV-1 protease to shift conformational sampling to stabilize open-like conformations, while maintaining the predominant semi-open conformation for activity.

Inhibition of HIV-1 protease: the rigidity perspective

MOTIVATION

HIV-1 protease is a key drug target due to its role in the life cycle of the HIV-1 virus. Rigidity analysis using the software First is a computationally inexpensive method for inferring functional information from protein crystal structures. We evaluate the rigidity of 206 high-resolution (2 Å or better) X-ray crystal structures of HIV-1 protease and compare the effects of different inhibitors binding to the enzyme.

RESULTS

Inhibitor binding has little effect on the overall rigidity of the protein homodimer, including the rigidity of the active site. The principal effect of inhibitor binding on rigidity is to constrain the flexibility of the β-hairpin flaps, which move to allow access to the active site of the enzyme. We show that commercially available antiviral drugs which target HIV-1 protease can be divided into two classes, those which significantly affect flap rigidity and those which do not. The non-peptidic inhibitor tipranavir is distinctive in its consistently strong effect on flap rigidity.

CONTACT

jack.heal@warwick.ac.uk; r.roemer@warwick.ac.uk

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Structural basis of drug resistance by genetic variants of HIV type 1 clade c protease from India

Using computer modeling of three-dimensional structures and structural information available on the crystal structures of HIV-1 protease, we investigated the structural effects of mutations, in treatment-naive and treatment-exposed individuals from India and postulated mechanisms of resistance in clade C variants. A large number of models (14) have been generated by computational mutation of the available crystal structures of drug bound proteases. Localized energy minimization was carried out in and around the sites of mutation in order to optimize the geometry of interactions present. Most of the mutations result in structural differences at the flap that favors the semiopen state of the enzyme. Some of the mutations were also found to confer resistance by affecting the geometry of the active site. The E35D mutation affects the flap structure in clade B strains and E35N and E35K mutation, seen in our modeled strains, have a more profound effect. Common polymorphisms at positions 36 and 63 in clade C also affected flap structure. Apart from a few other residues Gln-58, Asn-83, Asn-88, and Gln-92 and their interactions are important for the transition from the closed to the open state. Development of protease inhibitors by structure-based design requires investigation of mechanisms operative for clade C to improve the efficacy of therapy.

Do collective atomic fluctuations account for cooperative effects? Molecular dynamics studies of the U1A-RNA complex

A complete understanding of gene expression relies on a comprehensive understanding of the protein-RNA recognition process. However, the study of protein-RNA recognition is complicated by many factors that contribute to both binding affinity and specificity, including structure, energetics, dynamical motions, and cooperative interactions. Several recent studies have suggested that energetic coupling between residues contributes to formation of the complex between the U1A protein and stem loop 2 of U1 snRNA as a consequence of a cooperative network of interactions. We have performed molecular dynamics simulations on the U1A-RNA complex, including explicit water and counterions, and analyzed the results based on the calculated positional cross-correlations of atomic fluctuations. The results indicate that cross-correlations calculated on a per residue basis agree well with the observed inter-residue cooperativity and predict that the networks identified to date may also be coupled into an extensive hyper-network that reflects the intrinsic rigidity of the RNA recognition motif. In addition, we report a comparison of the MD calculated correlations with the results of a positional covariance analysis based on the sequences of 330 RNA recognition motifs, including U1A. The calculated inter-residue cross-correlations agree very well with the results of the sites exhibiting positional covariance. Collectively, these results strongly support the hypothesis that collective fluctuations contribute to cooperativity and the corresponding observed thermodynamic coupling. Predictions of additional sites in U1A that may be involved in cooperative networks are advanced.

Hydrophobic sliding: a possible mechanism for drug resistance in human immunodeficiency virus type 1 protease

Hydrophobic residues outside the active site of HIV-1 protease frequently mutate in patients undergoing protease inhibitor therapy; however, the mechanism by which these mutations confer drug resistance is not understood. From analysis of molecular dynamics simulations, 19 core hydrophobic residues appear to facilitate the conformational changes that occur in HIV-1 protease. The hydrophobic core residues slide by each other, exchanging one hydrophobic van der Waal contact for another, with little energy penalty, while maintaining many structurally important hydrogen bonds. Such hydrophobic sliding may represent a general mechanism by which proteins undergo conformational changes. Mutation of these residues in HIV-1 protease would alter the packing of the hydrophobic core, affecting the conformational flexibility of the protease. Therefore these residues impact the dynamic balance between processing substrates and binding inhibitors, and thus contribute to drug resistance.

Dynamical behavior of the vascular endothelial growth factor: Biological implications

The vascular endothelial growth factor (VEGF) seems to be the most important regulator of physiological and pathological angiogenesis, being, for this reason, a favorite target for therapies against angiogenesis-related diseases. VEGF is a homodimer in which the monomers are formed by beta-strands interconnected on the poles by three loops. A recent work showed that an intimate relationship between loops-1 and -3 is required for high affinity binding to the receptors (Kiba et al., J Biol Chem 2003;278:13453-13461). In this work, we report the results of a 10-ns molecular dynamics simulation of VEGF. We analyzed the dynamical behavior of the protein (using a dynamical cross-correlation map) and found that it is governed by a high degree of correlation between the motions of the loops. We also performed a principal component analysis and found an overall motion in which the opposite poles are projected against each other, just like the movement of the wings of a butterfly. From the biological point of view, it is likely that this motion would facilitate receptor binding since VEGF must enter a restricted cavity formed by the two subunits of the receptor.

High performance computing in biology: multimillion atom simulations of nanoscale systems

Computational methods have been used in biology for sequence analysis (bioinformatics), all-atom simulation (molecular dynamics and quantum calculations), and more recently for modeling biological networks (systems biology). Of these three techniques, all-atom simulation is currently the most computationally demanding, in terms of compute load, communication speed, and memory load. Breakthroughs in electrostatic force calculation and dynamic load balancing have enabled molecular dynamics simulations of large biomolecular complexes. Here, we report simulation results for the ribosome, using approximately 2.64 million atoms, the largest all-atom biomolecular simulation published to date. Several other nano-scale systems with different numbers of atoms were studied to measure the performance of the NAMD molecular dynamics simulation program on the Los Alamos National Laboratory Q Machine. We demonstrate that multimillion atom systems represent a 'sweet spot' for the NAMD code on large supercomputers. NAMD displays an unprecedented 85% parallel scaling efficiency for the ribosome system on 1024 CPUs. We also review recent targeted molecular dynamics simulations of the ribosome that prove useful for studying conformational changes of this large biomolecular complex in atomic detail.

A study of collective atomic fluctuations and cooperativity in the U1A-RNA complex based on molecular dynamics simulations

Cooperative interactions play an important role in recognition and binding in macromolecular systems. In this study, we find that cross-correlated atomic fluctuations can be used to identify cooperative networks in a protein-RNA system. The dynamics of the RRM-containing protein U1A-stem loop 2 RNA complex have been calculated theoretically from a 10 ns molecular dynamics (MD) simulation. The simulation was analyzed by calculating the covariance matrix of all atomic fluctuations. These matrix elements are then presented in the form of a two-dimensional grid, which displays fluctuations on a per residue basis. The results indicate the presence of strong, selective cross-correlated fluctuations throughout the RRM in U1A-RNA. The atomic fluctuations correspond well with previous biophysical studies in which a multiplicity of cooperative networks have been reported and indicate that the various networks identified in separate individual experiments are fluctuationally correlated into a hyper-network encompassing most of the RRM. The calculated results also correspond well with independent results from a statistical covariance analysis of 330 aligned RRM sequences. This method has significant implications as a predictive tool regarding cooperativity in the protein-nucleic acid recognition process.

Restrained molecular dynamics simulations of HIV-1 protease: the first step in validating a new target for drug design

To test the anticorrelated relationship that was recently displayed in conventional molecular dynamics (MD) simulations, several different restrained MD simulations on a wild type and on the V82F/I84V drug-resistant mutant of HIV-1 protease were performed. This anticorrelated relationship refers to the observation that compression of the peripheral ear-to-cheek region of HIV protease (i.e., the elbow of the flap to the fulcrum and the cantilever) occurred as the active site flaps were opening, and, conversely, expansion of that ear-to-cheek region occurred as both flaps were closing. An additional examination of this anticorrelated relationship was necessary to determine whether it can be harnessed in a useful manner. Consequently, six different MD experiments were performed that incorporated pairwise distance restraints in that ear-to-cheek region (i.e., the distance between the alpha-carbons of Gly40 and Gln61 was restrained to either 7.7 or 10.5 A, in both monomers). Pushing the backbones of the ear and the cheek regions away from each other slightly did force the flaps that guard the active site to remain closed in both the wild type and the mutant systems-even though there were no ligands in the active sites. Thus, these restrained MD simulations provided evidence that the anticorrelated relationship can be exploited to affect the dynamic behavior of the flaps that guard the active site of HIV-1 protease. These simulations supported our hypothesis of the mechanism governing flap motion, and they are the first step towards validating that peripheral surface as a new target for drug design.

Functional Plasticity in the Substrate Binding Site of β-Secretase

The aspartic protease beta-secretase (BACE) cleaves the amyloid precursor protein into a 42 residue beta-peptide, which is the principal biochemical marker of Alzheimer's disease. Multiple explicit-water molecular dynamics simulations of the apo and inhibitor bound structures of BACE indicate that both open- and closed-flap conformations are accessible at room temperature and should be taken into account for inhibitor design. Correlated motion is observed within each of the two lobes of BACE, as well as for the interfacial region. A self-inhibited conformation with the side chain of Tyr71 occupying the S(1) pocket is present in some of the unbound simulations. The reversible loss of the side chain hydrogen bond between the catalytic Asp32 and Ser35, due to the concomitant reorientation of the Ser35 hydroxyl group and a water molecule conserved in pepsin-like enzymes, provides further evidence for the suggestion that Ser35 assists in proton acceptance and release by Asp32 during catalysis.

The role of loop ZA and Pro371 in the function of yeast Gcn5p bromodomain revealed through molecular dynamics and experiment

Biological experiments were combined with molecular dynamics simulations to understand the importance of amino acidic residues present in the bromodomain of the yeast histone acetyltransferase Gcn5p. It was found that residue Pro371 plays an important role in the molecular recognition of the acetylated histone H4 tail by Gcn5p bromodomain. Crystallographic analysis of the complex showed that this residue does not directly interact with the histone substrate. It has been demonstrated that a double mutation Pro371Thr and Met372Ala in the Gcn5p bromodomain impairs chromatin remodeling activity. It is demonstrated here that, in this double mutant and in the fully deleted bromodomain strain, there is lower growth under amino acid deprivation conditions. By in vitro surface plasmon resonance (Biacore) experiments it is shown that the binding affinity of the double mutation to acetyl lysine 16 histone H4 peptide decreases. Molecular dynamics simulations were used to explain this loss in acetyl lysine-Gcn5p bromodomain affinity, in the double mutant. By comparing nanosecond molecular dynamics trajectories of the native as well as the single and doubly mutated bromodomain, it is concluded that the presence of Pro371 is important to the functionality of the Gcn5p bromodomain. In the simulation a point mutation involving this highly conserved residue induced an increase in the flexibility of the ZA loop, which in turn modulated the exposure of the binding pocket to the acetyl lysine. The combined double mutations (Pro371Thr-Met372Ala) not only markedly perturb the motion of the ZA loop but also destabilize the entire structure of the bromodomain.

Insights into saquinavir resistance in the G48V HIV-1 protease: quantum calculations and molecular dynamic simulations

The spread of acquired immune deficiency syndrome has increasingly become a great concern owing largely to the failure of chemotherapies. The G48V is considered the key signature residue mutation of HIV-1 protease developing with saquinavir therapy. Molecular dynamics simulations of the wild-type and the G48V HIV-1 protease complexed with saquinavir were carried out to explore structure and interactions of the drug resistance. The molecular dynamics results combined with the quantum-based and molecular mechanics Poisson-Boltzmann surface area calculations indicated a monoprotonation took place on D25, one of the triad active site residues. The inhibitor binding of the triad residues and its interaction energy in the mutant were similar to those in the wild-type. The overall structure of both complexes is almost identical. However, the steric conflict of the substituted valine results in the conformational change of the P2 subsite and the disruption of hydrogen bonding between the -NH of the P2 subsite and the backbone -CO of the mutated residue. The magnitude of interaction energy changes was comparable to the experimental K(i) data. The designing for a new drug should consider a reduction of steric repulsion on P2 to enhance the activity toward this mutant strain.

Investigation of glucose binding sites on insulin

Possible insulin binding sites for D-glucose have been investigated theoretically by docking and molecular dynamics (MD) simulations. Two different docking programs for small molecules were used; Multiple Copy Simultaneous Search (MCSS) and Solvation Energy for Exhaustive Docking (SEED) programs. The configurations resulting from the MCSS search were evaluated with a scoring function developed to estimate the binding free energy. SEED calculations were performed using various values for the dielectric constant of the solute. It is found that scores emphasizing non-polar interactions gave a preferential binding site in agreement with that inferred from recent fluorescence and NMR NOESY experiments. The calculated binding affinity of -1.4 to -3.5 kcal/mol is within the measured range of -2.0 +/- 0.5 kcal/mol. The validity of the binding site is suggested by the dynamical stability of the bound glucose when examined with MD simulations with explicit solvent. Alternative binding sites were found in the simulations and their relative stabilities were estimated. The motions of the bound glucose during molecular dynamics simulations are correlated with the motions of the insulin side chains that are in contact with it and with larger scale insulin motions. These results raise the question of whether glucose binding to insulin could play a role in its activity. The results establish the complementarity of molecular dynamics simulations and normal mode analyses with the search for binding sites proposed with small molecule docking programs.

Improved prediction of HIV-1 protease-inhibitor binding energies by molecular dynamics simulations

BACKGROUND

The accurate prediction of enzyme-substrate interaction energies is one of the major challenges in computational biology. This study describes the improvement of protein-ligand binding energy prediction by incorporating protein flexibility through the use of molecular dynamics (MD) simulations.

RESULTS

Docking experiments were undertaken using the program AutoDock for twenty-five HIV-1 protease-inhibitor complexes determined by x-ray crystallography. Protein-rigid docking without any dynamics produced a low correlation of 0.38 between the experimental and calculated binding energies. Correlations improved significantly for all time scales of MD simulations of the receptor-ligand complex. The highest correlation coefficient of 0.87 between the experimental and calculated energies was obtained after 0.1 picoseconds of dynamics simulation.

CONCLUSION

Our results indicate that relaxation of protein complexes by MD simulation is useful and necessary to obtain binding energies that are representative of the experimentally determined values.

Drug resistance in HIV-1 protease: Flexibility-assisted mechanism of compensatory mutations

The emergence of drug-resistant variants is a serious side effect associated with acquired immune deficiency syndrome therapies based on inhibition of human immunodeficiency virus type 1 protease (HIV-1 PR). In these variants, compensatory mutations, usually located far from the active site, are able to affect the enzymatic activity via molecular mechanisms that have been related to differences in the conformational flexibility, although the detailed mechanistic aspects have not been clarified so far. Here, we perform multinanosecond molecular dynamics simulations on L63P HIV-1 PR, corresponding to the wild type, and one of its most frequently occurring compensatory mutations, M46I, complexed with the substrate and an enzymatic intermediate. The quality of the calculations is established by comparison with the available nuclear magnetic resonance data. Our calculations indicate that the dynamical fluctuations of the mutated enzyme differ from those in the wild type. These differences in the dynamic properties of the adducts with the substrate and with the gem-diol intermediate might be directly related to variations in the enzymatic activity and therefore offer an explanation of the observed changes in catalytic rate between wild type and mutated enzyme. We anticipate that this "flexibility-assisted" mechanism might be effective in the vast majority of compensatory mutations, which do not change the electrostatic properties of the enzyme.

Relation between sequence and structure of HIV-1 protease inhibitor complexes: a model system for the analysis of protein flexibility

The flexibility of different regions of HIV-1 protease was examined by using a database consisting of 73 X-ray structures that differ in terms of sequence, ligands or both. The root-mean-square differences of the backbone for the set of structures were shown to have the same variation with residue number as those obtained from molecular dynamics simulations, normal mode analyses and X-ray B-factors. This supports the idea that observed structural changes provide a measure of the inherent flexibility of the protein, although specific interactions between the protease and the ligand play a secondary role. The results suggest that the potential energy surface of the HIV-1 protease is characterized by many local minima with small energetic differences, some of which are sampled by the different X-ray structures of the HIV-1 protease complexes. Interdomain correlated motions were calculated from the structural fluctuations and the results were also in agreement with molecular dynamics simulations and normal mode analyses. Implications of the results for the drug-resistance engendered by mutations are discussed briefly.

Molecular dynamics simulations of human rhinovirus and an antiviral compound

The human rhinovirus 14 (HRV14) protomer, with or without the antiviral compound WIN 52084s, was simulated using molecular dynamics and rotational symmetry boundary conditions to model the effect of the entire icosahedral capsid. The protein asymmetrical unit, comprising four capsid proteins (VP1, VP2, VP3, and VP4) and two calcium ions, was solvated both on the exterior and the interior to fill the inside of the capsid. The stability of the simulations of this large system (~800 residues and 6,650 water molecules) is comparable to more conventional globular protein simulations. The influence of the antiviral compound on compressibility and positional fluctuations is reported. The compressibility, estimated from the density fluctuations in the region of the binding pocket, was found to be greater with WIN 52084s bound than without the drug, substantiating previous computations on reduced viral systems. An increase in compressibility correlates with an entropically more favorable system. In contrast to the increase in density fluctuations and compressibility, the positional fluctuations decreased dramatically for the external loops of VP1 and the N-terminus of VP3 when WIN 52084s is bound. Most of these VP1 and VP3 loops are found near the fivefold axis, a region whose mobility was not considered in reduced systems, but can be observed with this simulation of the full viral protomer. Altered loop flexibility is consistent with changes in proteolytic sensitivity observed experimentally. Moreover, decreased flexibility in these intraprotomeric loops is noteworthy since the externalization of VP4, part of VP1, and RNA during the uncoating process is thought to involve areas near the fivefold axis. Both the decrease in positional fluctuations at the fivefold axis and the increase in compressibility near the WIN pocket are discussed in relationship to the antiviral activity of stabilizing the virus against uncoating.

Influence of transketolase substrates on its conformation

Dynamics stimulation of the holotransketolase molecule revealed that the enzyme's conformation in crystal was different from that in solution. It was shown also that dissolved holotransketolase can bind aldose (the acceptor substrate) even in the absence of ketose (the donor substrate). The holotransketolase conformation did not change upon aldose binding unlike in the case of ketose binding/cleavage. Therefore the conformation of a catalytic complex of holotransketolase with an intermediate-i.e., a glycolaldehyde residue formed upon binding and subsequent cleavage of ketose-differed, at least in solution, from the conformation of both the free and aldose-complexed holotransketolase. Some structural peculiarities of the holotransketolase with the intermediate were established by means of molecular dynamics stimulation.

Conformational flexibility of the catalytic Asp dyad in HIV-1 protease: An ab initio study on the free enzyme

The enzyme protease from the human immunodeficiency virus type 1 (HIV-1 PR) is one of the main targets for therapeutic intervention in AIDS. Computer modeling is useful for probing the binding of novel ligands, yet empirical force field-based methods have encountered problems in adequately describing interactions of the catalytic aspartyl pair. In this work we use ab initio dynamic methods to study the molecular interactions and the conformational flexibility of the Asp dyad in the free enzyme. Calculations are performed on model complexes that include, besides the Asp dyad, the conserved Thr26 and Gly27 residues and water molecules present in the active site channel. Our calculations provide proton location and binding mode of the active-site water molecule, which turn out to be different from those of the eukariotic isoenzyme. Furthermore, the calculations reproduce well the structural features of the aspartyl dyad in the protein. Finally, they allow the identification of both dipole/charge interactions and a low-barrier hydrogen bond as important stabilizing factors for the peculiar conformation of the active site. These findings are consistent with site-directed mutagenesis experiments on the 27, 27; positions (Bagossi et al., Protein Eng 1996;9:997-1003). The electric field of the protein frame (included in some of the calculations) does not affect significantly the chemical bonding at the cleavage site. Proteins 2000;39:26-36.

Three-dimensional structure of an Fab-peptide complex: structural basis of HIV-1 protease inhibition by a monoclonal antibody

F11.2.32, a monoclonal antibody raised against HIV-1 protease (Kd = 5 nM), which inhibits proteolytic activity of the enzyme (K(inh) = 35(+/-3)nM), has been studied by crystallographic methods. The three-dimensional structure of the complex between the Fab fragment and a synthetic peptide, spanning residues 36 to 46 of the protease, has been determined at 2.2 A resolution, and that of the Fab in the free state has been determined at 2.6 A resolution. The refined model of the complex reveals ten well-ordered residues of the peptide (P36 to P45) bound in a hydrophobic cavity at the centre of the antigen-binding site. The peptide adopts a beta hairpin-like structure in which residues P38 to P42 form a type II beta-turn conformation. An intermolecular antiparallel beta-sheet is formed between the peptide and the CDR3-H loop of the antibody; additional polar interactions occur between main-chain atoms of the peptide and hydroxyl groups from tyrosine residues protruding from CDR1-L and CDR3-H. Three water molecules, located at the antigen-antibody interface, mediate polar interactions between the peptide and the most buried hypervariable loops, CDR3-L and CDR1-H. A comparison between the free and complexed Fab fragments shows that significant conformational changes occur in the long hypervariable regions, CDR1-L and CDR3-H, upon binding the peptide. The conformation of the bound peptide, which shows no overall structural similarity to the corresponding segment in HIV-1 protease, suggests that F11.2.32 might inhibit proteolysis by distorting the native structure of the enzyme.

Domain motions in dihydrofolate reductase: a molecular dynamics study

Molecular dynamics simulations have been carried out on the enzyme dihydrofolate reductase from Lactobacillus casei complexed with methotrexate, NADPH and 264 crystallographic water molecules. Analysis of correlations in atomic fluctuations reveal the presence of highly correlated motion (correlation coefficient > 0.6) in the region between residues 30 to 35 and 85 to 90 leading to the identification of two domains, an "adenosine-binding domain" and a "large domain", which rotate by 3 to 4 degrees with respect to each other. The strongest correlation (> 0.6) within the large domain involves a coupling between the motions of the "teen-loop", and the spatially contiguous loops linking beta 6-beta 7 and beta 7-beta 8. Moreover, there is a significant correlation (approximately 0.5) between the adenosine fragment of NADPH and the pteridine and p-aminobenzoyl fragments of methotrexate, which are separated by approximately 17 A, and is lost on removal of "rigid-body" motion from the original trajectory. This provides support for the idea that the relative motion of the two domains is a means by which the occupation of the binding site for the adenosine end of the coenzyme can affect methotrexate binding and vice versa. Quasiharmonic vibrational analysis of the trajectory reveals that the overall dynamics of the system are governed by domain motions whose contributions are dominant at low frequencies. In addition, different low-frequency modes are responsible for separately coupling the adenosine-binding site and parts of methotrexate.

A pseudo-particle approach for studying protein-ligand models truncated to their active sites

A molecular dynamics method has been developed to describe the structural and dynamic properties of protein-ligand complexes that are truncated to their active sites. The active site is comprised of the ligand and discontinuous, positionally unrestrained peptide chains. This truncated active-site complex is surrounded by big unspecific pseudo-particles representing the complete protein and the solvent. Thus, knowledge of the folding of the outer parts of the protein is not required, and the method can be applied to protein models, derived from homology modeling. The method has been tested using ligand complexes of adenylate kinase, retinol binding protein, HIV-1 protease, and human leucocyte antigen. Comparisons with their crystal structures and with results from time-demanding simulations of the whole complexes in explicit water solvent show that the ligand binding properties are conserved. Most of the hydrogen bonds between the ligand and the active-site residues are reproduced and, furthermore, the simulation time is reduced.

Resistance of HIV type 1 to proteinase inhibitor Ro 31-8959

During replication of human immunodeficiency virus type 1 (HIV-1), proteolytic cleavage of Gag and Gag-Pol precursor proteins into different functional protein subunits is catalyzed by the viral proteinase, and this enzyme is the target of the antiviral proteinase inhibitor, Ro 31-8959. We investigated in vitro which HIV mutants with reduced sensitivity to Ro 31-8959 emerged during proteinase inhibition treatment; from three different HIV-1 strains, comparable progeny virus resistant to proteinase inhibitor were found, whereas the same experimental protocol detected no resistant HIV-2 mutants. Molecular analysis of the mutations underlying resistance revealed a multistep mechanism in which an amino acid exchange was common to all resistant isolates, and in all experiments preceded further exchanges at position 90 (leucine to methionine) and/or at position 54 (isoleucine to valine). For wild-type strains the 90% inhibitory concentrations of Ro 31-8959 were close to 20 nM, whereas HIV-1 mutants with all 3 amino acid exchanges had more than 50-fold increased 90% inhibitory concentrations (above 1000 nM). The primary event (Gly-48 to valine) occurs at the hinge of the flaps of the proteinase, thus hampering entry of the inhibitor to the active center and suggesting steric hindrance. Detailed knowledge of this stereotypic process could open inhibitor design, thus preventing conceivable escape of resistant virus on proteinase inhibitor action.

Methodological advances in molecular dynamics simulations of biological systems

New advances in the techniques used to simulate specific statistical ensembles provide molecular dynamics algorithms that permit rigorous connections to be made between thermodynamic observables and calculated quantities in simulations of biological molecules. The complete inclusion of electrostatic forces in simulations also improves the comparison between calculations of simple structural measures in crystals and X-ray crystallographic results. These advances coupled with the ongoing development of more accurate inter/intramolecular forcefields with enhanced accuracy provide guidelines and benchmarks for comparison as we move to study more complicated biological problems.

Molecular dynamics simulations of trp apo- and holorepressors: domain structure and ligand-protein interaction

Molecular dynamics simulations of the apo- and holo-forms of the trp-repressor protein were performed under extensively solvated conditions in order to elucidate their dynamic structures and ligand-protein interactions. The root mean square fluctuations calculated from the trajectories agreed with those calculated from X-ray temperature factors. Distance, distance fluctuation, and dynamic cross-correlation maps were drawn to provide information on the dynamic structures and communications among the domains. A three-domain format has been proposed for the crystal structure (Zhang et al., Nature 327:591-597, 1987); namely, helices A-C and F of both subunits make up a central core, and D and E of each subunit forms a DNA binding head. The results of the simulations were mostly consistent with the three-domain format. However, helix F was more flexible and freer than other parts of the central core. The turn DE, the helix-turn-helix DNA binding motif, was free from interactions and correlations with other domains in both forms of the repressor. A comparison of the simulations of the aporepressor and holorepressor showed that tryptophan binding made the DNA-binding helix D more flexible but helix F less flexible. Several amino acid residues in contact with the bound tryptophan were identified as making concerted motions with it. Interaction energies between the corepressor and the amino acid residues of the protein were analyzed; the results were mostly consistent with the mutational experiments.

Molecular dynamics studies of the human CD4 protein

A dynamical model for an N-terminal fragment of the human CD4 protein has been determined by computer simulation. The protein has been studied both in vacuo and in solution. Data from both simulations agree moderately well with each other and with the crystal structure. All elements of secondary structure were retained during simulation. Point mutation and sequence replacement studies have shown that a loop in CD4, residues 40-52 is involved in binding with gp120, the human immunodeficiency virus surface glycoprotein. Our results show that the gp120-binding loop and a few regions which bind to monoclonal antibodies and class II MHC molecules are the most highly motile areas of the protein. These results are consistent with the suggestion that CD4 binds to target molecules by using induced-fit contacts.

Nucleotide-binding properties of adenylate kinase from Escherichia coli: a molecular dynamics study in aqueous and vacuum environments

The complex of adenylate kinase with its transition-state inhibitor has been studied by molecular dynamics simulations in water and in vacuum environments with the GROMOS force field over a period of 300 ps. The adenylate kinase, a member of the nucleotide-binding protein family, was exemplarily chosen for the inspection of the nucleotide-binding properties in the active site. The ligand binding and the domain movements have been studied in detail over the simulation period and compared with the crystal structure. Secondary structure transitions and domain closures defined those parts of the structure which are involved in an induced-fit movement of the enzyme. The presence of more stable hydrogen bonds on the substrate side leads to the assumption that substrate binding is more specific than cosubstrate binding. Reliable results were achieved only if water was explicitly included in the stimulation.

Synthesis of novel inhibitors of the HIV-1 protease: difunctional enols of simple N-protected amino acids

A series of enol HIV-1 protease inhibitors which show competitive inhibition and the structure-activity relationship study which led to the design of these compounds are reported. By systematically modifying simple amino acids, Boc-Phe enol and Boc-Tyr enol derivatives yield nanomolar Kiapp values (Kiapp = 0.485 microM and Kiapp = 0.425 microM, respectively). These enols are of low molecular weight (< 500 g/mol) and of non-peptidic nature. The enols are synthesized in a one step chemical synthesis and modifications to increase their potency could easily be performed. Boc-Phe enol and Boc-Tyr enol showed low inhibitory effect on pepsin, Kiapps of 23 and 149 microM, respectively, and Boc-Phe enol showed a Kiapp of 20 microM for cathepsin D. Neither of these two compounds inhibited renin (< 10% inhibition at 200 microM).

Probing structure-function relationships in human immunodeficiency virus type 1 protease via molecular dynamics simulation

This chapter has focused on the application of molecular dynamics computer simulations and related molecular modeling techniques to the study of HIV protease structure and structure-function relationships. The abundance of crystallographic data provides ample experimental quantities (average structures, temperature factors, and hydrogen bond topography) to validate the computational techniques employed. Furthermore, these studies provide insight into the structure and functional energetics of HIV-1 protease that would be difficult or impossible to study experimentally. This chapter covers studies that investigate correlated motion between and within subunits of the protease, mutants of the protease that disrupt the tertiary structure and dimer formation, and studies of HIV-1 protease-inhibitor complexes that rationalize both the protonation state of the active site and the observed binding strength of these complexes. These studies demonstrate that MD is capable of contributing to our understanding of structure-function relationships and may aid in the design of potential therapeutics.

Studies on the specificity of HIV protease: an application of Markov chain theory

A sequence-coupled (Markov chain) model is proposed to predict the cleavage sites in proteins by proteases with extended specificity subsites. In addition to the probability of an amino acid occurring at each of these subsites as observed from a training set of oligopeptides known cleavable by HIV protease, the conditional probabilities as reflected by the neighbor-coupled effect along the subsite sequence are also taken into account. These conditional probabilities are derived from an expanded training set consisting of sufficiently large peptide sequences generated by the Monte Carlo sampling process. Very high accuracy was obtained in predicting protein cleavage sites by both HIV-1 and HIV-2 proteases. The new method provides a rapid and accurate means for analyzing the specificity of HIV protease, and hence can be used to help find effective inhibitors of HIV protease as potential drugs against AIDS. The principle of this method can also be used to study the specificity of any multisubsite enzyme.