In conventional molecular dynamics simulation, generally, only the first structure including indivisible molecular fragments would be accounted for. Herein, it can be seen that the molecular fragments of the enzyme active site are preferentially in their ionized and polarized forms. Over the whole trajectory the main contribution arises from the 5th structure that includes the ionized water molecule and the peptide bond with the C==O bond polarized. It is worth noting that the first structure on the average gives only 10% contribution. The 2nd and 4th structures that give on average contributions of ~30 and 10%, respectively, also include either an ionized water molecule or a fragment with a polarized C==O bond.

In the next step eighteen valence structures of the first reaction model (see section AVB Description of the System) were used allowing additionally for the proton transfer from the protonated Asp-25 onto the peptide bond oxygen. A 100-ps run was performed starting from the thermalized coordinates and velocities. The 100-ps simulation yields the rms deviation values from the crystal structure of 1.78, 1.79, and 2.24 Å, and from the averaged thermalized structure of 1.37, 1.37, and 1.78 Å for the trace, backbone, and heavy solute atoms, respectively. These rms deviations in comparison with the MD/AVB simulations accounting only for eight valence structures, are ~0.3 lower with respect to the crystal structure and ~0.1 higher with respect to the averaged thermalized structure.

During this run the proton transfer did not occur. The fluctuations of the distances are similar to those of the previous simulation. There are, however, different contributions of the valence structures. Due to one more possible ionized form, that of Asp-25, the dominant structures are 8th and 11th, shown below, with the c equal on average 0.29 ± 0.02 and 0.37 ± 0.03, respectively.

8.	A2 + H + OH + S + H + A1
11.	A2 + H + OH + S + H + A1

These dominant structures contain the same elements (i.e., the ionized or polarized molecular fragments) as in the previous calculations with eight valence structures. Additionally, they include an ionized form of A1. There are also some smaller contributions of ~8% of each of the structures: 5th, 9th, and 10th.

5.	A1 + OH + H + A2 + S
9.	OH + A2 + S + H + A1
10.	H2O + A2 + S + H + A1

	INTRODUCTION

TOP ABSTRACT INTRODUCTION AVB DESCRIPTION OF THE... SIMULATION METHODS OPTIMIZATION, THERMALIZATION,... RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

HIV-1 protease (HIV-1 PR) is one of the three enzymes encoded by the viral genome. It exhibits its main role during the course of viral maturation. HIV-1 PR processes the gag and pol polyproteins, which further generate important viral enzymes and structural proteins, including HIV-1 PR itself. Human immunodeficiency virus (HIV) is a causative agent of acquired immunodeficiency syndrome (AIDS) disease. It was observed that virions that lack HIV-1 PR are noninfectious (Kohl et al., 1988). Hence, agents that inhibit this enzyme are used to treat patients infected with the HIV virus and may be crucial in the chemotherapy of AIDS. Therefore, for many years people have shown interest in this enzyme, its inhibitors, and also the catalytic mechanism.

HIV-1 PR catalyses the cleavage of peptide bonds. It acts and crystallizes as a homodimer, with 99 residues in each monomer, one related to the other by a crystallographic twofold axis (Wlodawer et al., 1989). The substrate cleavage site is situated at the interface between identical, adjacent subunits. It involves a pair of a catalytic triad, Asp-Thr-Gly, which identifies HIV-1 PR as an aspartyl protease. The two aspartyl side chains, one from each monomer, are closely apposed and form a symmetric dyad. The binding site has a shape of an extended ravine with an entry controlled by two flaps that form a flexible gate for an approaching ligand (Harte et al., 1990; Swaminathan et al., 1991; Ishima et al., 1999; Scott and Schiffer, 2000).

There have been several proposals for the reaction mechanism (e.g., Hyland et al., 1991a,b; Rodriguez et al., 1993; Chatfield and Brooks, 1995; Lee et al., 1996; Silva et al., 1996; Chatfield et al., 1998; Park et al., 2000). Experimental data (Hyland et al., 1991a,b; Rodriguez et al., 1993) indicate the contribution of the lytic water molecule in the catalytic process. Most studies (e.g., Hyland et al., 1991a; Silva et al., 1996) indicate an acid-base mechanism in which the peptide bond hybridization is changed from sp² to sp³ through a nucleophilic attack on its carbonyl group by a lytic water molecule. To elucidate the reaction mechanism ab initio studies on model fragments (Lee et al., 1996; Park et al., 2000) and also classical molecular dynamics (MD) simulations on the HIV-1 PR system (Chatfield and Brooks, 1995; Chatfield et al., 1998) were performed. Some studies with molecular mechanics (MM) and combined quantum-mechanics/molecular mechanics (QM/MM) potential energy function were also carried out (Chatfield et al., 1998; Liu et al., 1996). For an overview of the methods on computational studies of enzymatic reactions, see for instance Field (2002). However, it still remains unclear what is the detailed mechanism of the reaction catalyzed by HIV-1 PR, and there are several questions regarding the catalytic process that cannot be answered by experimental means. First issue is connected with the protonation states of ionizable amino acids and especially of the catalytic dyad in the enzyme with a bound substrate. The activity of HIV-1 PR varies with pH to yield a "bell-shaped" graph of the catalytic rate constant versus pH with a maximum at approximately pH 5 to 6 (Hyland et al., 1991a; Ido et al., 1991). The customary explanation of this fact is that the active form is monoprotonated and exists only at intermediate pH. Nevertheless, the location of protons in the active site may determine the reaction path. This problem has been addressed before experimentally (Hyland et al., 1991a; Ido et al., 1991; Smith et al., 1996; Yamazaki et al., 1994; Wang et al., 1996) and theoretically (Harte and Beveridge, 1993; Chatfield and Brooks, 1995; Geller et al., 1997; Trylska et al., 1999; Piana and Carloni, 2000; Piana et al., 2001). These studies were performed either for the native enzyme or complexed with various inhibitors. Herein, we present theoretical titration results for the complex of HIV-1 PR with a model substrate. Other still questionable points are, e.g., the position of the water molecule in the cleavage site, the mobility of the active site environment, the stability of protein-substrate interactions, and the protonation state of the reaction intermediate. Some of them are addressed in this study.

We have parameterized the approximate valence bond (AVB) method, very fast quantum-mechanical generator of the potential energy surface, for the HIV-1 PR catalyzed reaction, based on the density functional theory (DFT) calculations. This procedure is exactly described in Trylska et al. (2001). It allowed us to obtain the parameters needed to generate the potential energy surface for the active region where the breakage and formation of the bonds occurs during the catalytic process. These parameters are applied in the MD/AVB dynamics.

An empirical version of the valence bond (VB) method was first proposed and developed by Warshel (Warshel and Weber, 1980; Warshel, 1991). The similar AVB approach uses ab initio type calculations and a different parameterization strategy (Grochowski et al., 1996; Trylska et al., 2001). The method was successfully applied for simulations of the whole enzymatic reaction catalyzed by phospholipase A₂ (Bala et al., 1998, 2000).

MD/AVB simulations of the first steps of the reaction were performed in this study for HIV-1 PR complexed with a model substrate with explicit water and mobile ions. The nucleophilic attack is a rate-limiting step and to observe the reaction in a real, shorter time scale than applicable in the dynamics, the potential energy surface for one degree of freedom was deformed with a harmonic potential as in the umbrella sampling technique of Torrie and Valleau (1974). Then the additional restraints were removed, the reaction did not go back, and next proton transfer was observed.

	AVB DESCRIPTION OF THE SYSTEM

TOP ABSTRACT INTRODUCTION AVB DESCRIPTION OF THE... SIMULATION METHODS OPTIMIZATION, THERMALIZATION,... RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

The modeled HIV-1 PR structure was divided into two regions. The quantum one, described with the AVB method, comprises molecular fragments involved directly in the enzymatic reaction, i.e., side chains of catalytic aspartates (25 and 125, residues of one monomer are usually numbered 1-99 and of the other 101-199; Fig. 1), the lytic water, and hydrolyzable peptide bond between both methionines (Met-203 and Met-204). The second, classical region consists of the surrounding protein and solvent environment and will be described in the MD/AVB simulations with a conventional force field.

View larger version (12K): [in this window] [in a new window]   FIGURE 1   AVB region with molecular fragments involved in the reaction: peptide bond between Met-203 and Met-204, catalytic water (207), and two aspartate side chains. Balls represent carbon atoms bonded to other fragments of the protein or the substrate.

To parameterize the AVB method, one needs to perform QM calculations in vacuum for various molecular fragments involved in the reaction. The chosen fragments are presented in Fig. 2 and they may be divided into four groups: 1) the side chains of catalytic aspartates with C carbon modified into the CH₃ group (A1, A1, A2, A2); 2) the catalytic water molecule (H₂O, OH, H); 3) the peptide bond and the transition state with C carbons modified into the CH₃ group (S, S, I1, I2, I3); and 4) the products (P1, P1, P2, P3). The QM calculations were performed using the DFT theory implemented in the Gaussian'94 program (Frish et al., 1995). The B3LYP exchange and correlation functional (Becke, 1993; Lee et al., 1988) was applied with 6-31G(d,p) and 6-31G+(d,p) basis sets. The AVB method and the parameterization procedure has been described in detail in our previous work (Trylska et al., 2001).

View larger version (24K): [in this window] [in a new window]   FIGURE 2   Molecular fragments of the enzyme reaction center used to reproduce the most frequently proposed reaction path. The key protons are labeled with letters: w, water; a, acid. Arrows denote the direction of the reaction. (R is CH3).

The schematic representation of the catalytic mechanism most often appearing in literature and the fragments used to build the valence bond structures in our model are presented in Fig. 2. First stage involves a proton transfer from the lytic water molecule onto Asp-125 (A2 in Fig. 2) and the nucleophilic attack of the hydroxy anion on the peptide bond carbon with the change in the hybridization of that bond. Second, proton transfer between Asp-25 (A1 in Fig. 2) and the carbonyl oxygen of the peptide bond occurs followed by or together with a conformational transition of the CN bond allowing the nitrogen lone pair to accept the hydrogen of Asp-125. Next, the CN bond breakage occurs (in our model through a polarized P1 molecule) and one of the hydrogens is transferred back from one of the hydroxyl groups of the tetrahedral intermediate onto Asp-25.

In the AVB method one constructs a matrix of the electronic Hamiltonian expressed in the basis of valence bond structures that represent different arrangements of covalent and ionic bonds in a molecular system. The ground many-electron state is approximated by a linear combination of the electronic wave functions _I corresponding to various valence structures

(1)

In case of the HIV-1 PR active site the minimal number of valence structures allowing to reproduce the mechanism presented in Fig. 2 is 23 (listed below). The plus signs between the fragments are symbolic and are included just for clarity.

1.	A1 + H2O + A2 + S
2.	A1 + OH + H + A2 + S
3.	A1 + OH + A2 + S
4.	A1 + H2O + A2 + S
5.	A1 + OH + H + A2 + S
6.	A1 + OH + A2 + S
7.	H2O + A2 + S + H + A1
8.	A2 + H + OH + S + H + A1
9.	OH + A2 + S + H + A1
10.	H2O + A2 + S + H + A1
11.	A2 + H + OH + S + H + A1
12.	OH + A2 + S + H + A1
13.	A1 + A2 + H + I1
14.	A1 + A2 + I1
15.	A2 + H + H + A1 + I1
16.	A2 + H + A1 + I1
17.	A2 + H + A1 + I2
18.	A2 + A1 + I2
19.	A2 + A1 + I3
20.	A2 + A1 + P1 + P3
21.	A2 + A1 + P1 + P3
22.	A1 + A2 + P2 + P3
23.	A2 + H + A1 + P2 + P3

Even if an acid-base catalysis with the lytic water donating the hydroxy anion for the nucleophilic attack is assumed, there still are various possibilities of the formation of the reaction intermediate and its breakage. The first model (Fig. 2) proposed that after the I1 molecule is formed, H proton is transferred onto one of the oxygens that results in I2. This is followed by the H proton transfer onto the nitrogen and I3 is formed. However, after the formation of I1 some other paths may be possible what is shown in Fig. 3. The I1 molecule may first accept H proton on its nitrogen. In such case the I2' molecule is formed and it can either dissociate into products (in our case through a polarized molecule P2) or accept the H proton on the peptide oxygen O and form the I3 molecular fragment. The I3 molecule dissociates into products in the same manner as presented in Fig. 2.

View larger version (13K): [in this window] [in a new window]   FIGURE 3   Set of molecular fragments used to reproduce other possibilities of the reaction path after the formation of the I1 molecule. The key protons are labeled with letters: w, water; a, acid. Arrows denote the direction of the reaction. (R is CH3).

Considering the path I1  I2'  I3 and then the dissociation into products, also 23 structures describe the whole reaction, as shown above, but instead of structures 17th and 18th, 17'th, and 18'th were used:

17'.	A2 + H + A1 + I2'

18'.  A1 + A2 + I2'

Considering the path I1  I2'  P2 + P3, a total of 22 structures model the reaction with the first sixteen as presented above. The further six structures are listed below:

17'.	A2 + H + A1 + I2'

18'.  A1 + A2 + I2'

19.  A1 + A2 + P2 + P3

20.  A2 + H + A1 + P2 + P3

21.  A1 + A2 + P2 + P3

22.  A2 + H + A1 + P2 + P3

Therefore, depending on the choice of the set of structures one may model various possible reaction mechanisms. Of course, all the models may be accounted for and the whole set of molecular fragments appearing in Figs. 2 and 3 may be used to construct the set of AVB structures and further on to perform MD/AVB simulations. However, such simulations would be far more time consuming.

	SIMULATION METHODS

TOP ABSTRACT INTRODUCTION AVB DESCRIPTION OF THE... SIMULATION METHODS OPTIMIZATION, THERMALIZATION,... RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

The AVB code was incorporated into Gromos'96 package (van Gunsteren et al., 1996; Scott et al., 1999) and the MD/AVB simulations were performed as described earlier with a system divided into two regions: quantum AVB and classical. The energy of all interactions in the AVB region and its electrostatic and polarization interactions with the classical region were calculated within the AVB formalism. The energy of all atomic interactions inside the classical region, as well as the Lennard-Jones potential and bonding interactions across the border between the classical and AVB regions, were calculated with the conventional Gromos'96 force field. The total potential energy surface of the system is, therefore, a sum of both the classical and AVB energy terms. The MD/AVB simulations were performed for the system that consisted of the enzyme, substrate, explicit water, and ions.

The enzyme-substrate complex was prepared using the crystallographic structure of HIV-1 PR complexed with the MVT-101 inhibitor [N-acetyl-Thr-Ile-Nleu-[CH₂NH]-Nleu-Gln-Arg-amide] (Miller et al., 1997). This inhibitor was constructed on the basis of the hydrolyzed X/p9 junction [Thr-Ile-Met-/-Met-Gln-Arg] of the HIV-1 polyprotein. Therefore, to mimic a substrate, each norleucine of the MVT-101 inhibitor was replaced by a methionine, and the reduced [CH₂NH] link was modified into a regular peptide bond. The lytic water was also added to the structure because it was not present in the HIV-1 PR:MVT-101 complex. This is probably due to the fact that transition state analogues have no space for the water molecule next to the catalytic aspartic acid residues.

Preceding MD/AVB simulations the pK_as of all titratable residues in the HIV-1 PR complexed with a model substrate were calculated based on the Poisson-Boltzmann model. The methodology has been described previously (Antosiewicz et al., 1994, 1996a,b; Trylska et al., 1999). Two titration models differing in parameter sets were applied. The full-charge/4 model assigns to the solute interior a dielectric constant of 4 and uses a detailed model of ionization, whereas the single-site/20 model assigns a dielectric constant of 20 and uses a simplified model of ionization. The structure for calculations was prepared as described above. The positions of hydrogens were generated with the HBUILD (Brunger and Karplus, 1988) command of CHARMm22 (CHARMm Version 22., 1992, Molecular Simulations Inc., Waltham, MA). The electrostatic calculations were performed with the UHBD program (Davis et al., 1991; Madura et al., 1995). The temperature was set to 293 K, ionic strength to 150 mM, and solvent dielectric constant to 80. Other technical details concerning the two titration models are presented elsewhere (Trylska et al., 1999).

	OPTIMIZATION, THERMALIZATION, AND EQUILIBRATION PROCEDURES

TOP ABSTRACT INTRODUCTION AVB DESCRIPTION OF THE... SIMULATION METHODS OPTIMIZATION, THERMALIZATION,... RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Preceding the molecular dynamics simulations of the modeled complex a special relaxation procedure had to be applied in the active site region due to the addition of a catalytic water molecule (207) and earlier mentioned alterations of the inhibitor. This procedure was performed for a system in vacuum to allow for the flap and substrate surrounding movement to adjust the position of the water molecule 207 and the substrate. Polar hydrogens were added to the protein, substrate and crystallographic water molecules with the HBUILD command of CHARMm22 (Brooks et al., 1983; CHARMm Version 22., 1992 Molecular Simulations Inc.). Second, four minimization cycles (100 steps of the steepest descents method followed by 200 steps of the conjugate gradients minimization method) with the Gromos'96 package were applied: first to the positions of crystallographic water hydrogens; second, to protein hydrogens; next, to the positions of all atoms of crystallographic waters (except 207 and 301); and finally to positions of the active site atoms of Asp-25 and -125, the substrate residues Met-203 and -204, H₂O-207 and 301, and the residues of the flap region 48 to 53 and 148 to 153. The minimizations were based on the Gromos'96 43B1 force field with a dielectric constant of 1. However, for Arg-8, Asp-25, Asp-29, Arg-108, Asp-129, and Asp-125 residues Gromos'96 43A1 force field was used, and charges of Glu-21, Glu-121, Lys-41, Lys-141, Lys-55, Lys-155, Lys-70, and Arg-114 side chains on the surface were set to zero. -Amino-N-butyric acid group that replaces Cys-67 and -95 in each monomer was build on the basis of alanine, and the charge on the additional C carbon was set to zero. The modified inhibitor is neutralized from one side by the acetyl group CH3-CO and from the other side by the CO-NH2 group. The CO charges were taken as in the peptide bond, for the CH3 group was set to zero, and for the NH2 group as in the Gln and Asn residues.

The prepared structure was then subject to thermalization in vacuum. The potential energy surface included the AVB parameterization in the active site region (for parameters see Trylska et al. (2001)) and Gromos'96 force field parameters for the rest of the enzyme. However, this thermalization attempt was not successful even though some alterations of the standard Gromos'96 force field parameters to neutralize the system have been applied as listed above. When the temperature reached 200 K in the canonical (NVT) ensemble the modeled system became unstable and the MD simulations in vacuum could not be continued. It was therefore decided to perform MD/AVB simulations with explicit water and ions.

Following such conclusions the protein was solvated in a truncated octahedral box of simple point change (SPC) water molecules (Smith and van Gunsteren, 1994). The width of the box was set to 7 Å from the border of the solvent accessible surface. Most (94 of 113) water molecules observed in the crystal structure were preserved. It involved the waters buried inside the protein and those situated close to the surface and hydrogen bonded to the amino acids on the surface of the enzyme. In total, the box was filled with 5572 water molecules. The parameters of the force field used for the classical region were changed to those used for the system with explicit solvent (43A1). The dielectric constant of 1 was used.

A 20-ps MD simulation of the solvent at 300-K temperature (with the protein fixed) was performed with a time step of 0.25 fs. Velocities were generated every 5 ps from the Maxwellian distribution at temperature of 100 K. The relaxation time of the thermal bath was set to 0.01 ps. The twin-range method (Berendsen, 1985) was used for the nonbonded interactions. They were evaluated every step with a short range cutoff radius of 12 Å. The pair list was updated every 10 steps. The longer range cutoff radius used for nonbonded interactions evaluated less frequently (every 10th simulation step) and in the reaction-field calculations was set to 14 Å. The truncated octahedron periodic boundary conditions were used.

To neutralize the modeled system 14 water molecules were replaced by 10 Cl and 4 Na⁺ ions. The ions were located at positions of the chosen water oxygens. They were added to neutralize the singly charged residues on the surface of the protein and the groups of three charged amino acids situated close to one another.

Then to adjust both the solvent and ions to the protein crystallographic structure MD simulations of the water and ionic environment were carried out. Ten 0.1-ps runs at 300 K with a time step of 0.25 fs with velocities reassigned each time at temperature of 100 K were performed. Only the water and ion atoms were free to move. Next a 20-ps dynamics of the solvent and ions at 300 K with velocities reassigned every 5 ps was carried out.

Next, the structure underwent the standard thermalization procedure with the potential energy surface including the AVB parameters for the active site and Gromos'96 force field for the rest of the system. The MD/AVB runs were performed for the first valence bond structure: A1 + H₂O + A2 + S (see AVB Description of the System). To allow for some flexibility of the system the solvent atoms were positionally restrained with a harmonic force constant of 2.5 × 10⁴ [kcal/mol × Å²]. A 5-ps run in the NVT ensemble at 10-K temperature was performed (step I). The time step in this and all further MD simulations was set to 0.5 fs. Then a 30-ps run (steps II-VII) with velocities generated every 5 ps and temperature increased by 50 K from 50 to 300 K, also every 5 ps was carried out. The solute atoms were free to move, however, some additional constraints were applied to some distances in the active site as listed for each step (I-VIII) in Table 1. The constraints were applied to conserve the hydrogen bonds in the active center and to keep the position of the lytic water between the aspartic acids during the thermalization. The distances Gly-27:N-Asp-25:OD1, Gly-27:N-Asp-125OD1, Gly-127:N-Asp-25:OD1, Gly-127:N-Asp-125:OD1 were restrained using the harmonic potential energy functions with a changeable force constant at a reference distance of 3 Å. The above hydrogen bonds and strong dipoles of the Thr-26-Gly-27 and Thr-126-Gly-127 pointing toward the negative charge of aspartyl dyad are said to be responsible for maintaining the position and planarity of aspartates (Piana and Carloni, 2000). The other active site atoms listed in Table 1 were restrained using one-sided harmonic potential energy functions, i.e., the restraining function was applied only when the actual distance between the restrained atoms was larger than the reference distance.

The last step (VIII, 120 ps) was to gradually free the solvent atoms and to allow for the relaxation of the whole system. Two 10-ps simulations with the harmonic force constant, applied for restraining the solvent atoms, decreased by 2 to 1.25 × 10⁴ [kcal/mol × Å²] and by 4 to 0.625 × 10⁴ [kcal/mol × Å²], respectively, and a 100-ps simulation at 300 K with all atoms free to move and the constraints decreased (see step VIII in Table 1) were performed. The relaxation time of the thermal bath was set to 0.1 ps.

The root mean square (rms) deviation of the averaged thermalized structure (30 coordinates of the last 15 ps of thermalization taken at equal intervals) with the crystallographic structure taken as a reference was calculated. The rms values were equal to 1.78, 1.80, and 2.14 Å for the C carbons (408 atoms), backbone (1616), and the heavy solute (3146) atoms, respectively.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION AVB DESCRIPTION OF THE... SIMULATION METHODS OPTIMIZATION, THERMALIZATION,... RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Protonation states in the enzyme-substrate complex

Both titration models presented in the Simulation methods section result in aspartic and glutamic acids negatively charged and lysines and arginines positively charged. On the average histidine pK_as are shifted up by 0.3 units. For brevity the results in Table 2 are presented only for the aspartic dyad. Catalytic water (207) in the active site and water 301 corresponding to the water linking the MVT-101 inhibitor with the flaps (Miller et al., 1997) are either included explicitly or as a test case these waters are removed and treated as a part of bulk solvent.

The full-charge/4 calculations show that in the presence of the model substrate Asp-25 possesses a hydrogen. The pK_as shift in comparison with the same data for the HIV-1 PR:MVT-101 complex (Trylska et al., 1999) that were 10.8 and 0.6 for Asp-25 and -125, respectively, with all crystal waters included. The data suggest that the probability of protonation of Asp-25 in the complex with the substrate is higher than in the complex with the MVT-101 inhibitor.

The single-site/20 calculations also shift the pK_as but only a few units. However, they still indicate that both aspartates are deprotonated at pH 7. This difference may result from the fact that the single-site/20 model is not sensitive to small changes of conformation and leads to underestimates of charge-charge interactions. This may be of great importance for the dyad situated inside the protein and being inaccessible to bulk solvent. The full-charge/4 calculations are consistent with the acid-base model of the reaction mechanism in which one of the aspartyl groups has to be protonated.

The results of the full-charge/4 model are in agreement with experimental data of Hyland and coworkers (Hyland et al., 1991a), which indicate that substrates bind only to a form with one of the catalytic aspartyl residues protonated. It is worth noting that some pK_as computed with the full-charge/4 model show exaggerated shifts that are clearly unrealistic because the protein is expected to denature at a much lower pH than 20 (Table 2). Such values are interpreted as indications that the group remains fixed in a single ionization state at any pH for which the complex is actually stable in solution.

MD/AVB simulations

The first step of the reaction, i.e., the dissociation of the lytic water molecule and the nucleophilic attack of the resulting hydroxy anion on the peptide bond carbon involved only eight valence structures in the AVB region (1-6 and 13, 14 of the first proposed model, see AVB Description of the System). A 100-ps MD run of the whole system in the NVT ensemble at 300 K with no constraints applied was performed. Fig. 4 shows the changes of the rms deviations of the system from the modified crystal structure and the averaged thermalized structure. The data indicate that there are no significant changes in the geometry of the system during this simulation. The proton transfer did not occur during this 100-ps dynamics but obviously the time scale was too short. However, the distance fluctuations in the active site are significant and in a larger time scale it may be possible that the proton is transferred onto Asp-125, and a hydroxy anion is formed. At some steps of the dynamics HW2 proton transfer onto OD2 of Asp-125 was observed but the proton did not stay at OD2 for a long time and the nucleophilic attack did not follow (for atom names appearing in this and further sections, see Fig. 1). One should, nevertheless, emphasize that in our MD/AVB dynamics proton was allowed to sample more broad configurational space regions than in classical dynamics.

View larger version (23K): [in this window] [in a new window]   FIGURE 4   Changes of the rms value during the 100 ps of the MD/AVB simulation. Filled and opaque symbols show the rms deviation from the modified crystal and averaged thermalized (30 coordinates of the last 15 ps of thermalization taken at equal intervals) structures, respectively. The rms value was calculated every 2.5 ps (last 30 frames were taken at equal intervals and average coordinates were generated). Circles refer to the main chain atoms (1616) and squares to heavy solute atoms (3146).

The minimal, maximal, and mean values with the standard deviations of the above distances during the 100-ps simulation are shown in Table 3.

If the nucleophilic attack is to occur, the change in the hybridization of the peptide bond carbon (C) and the formation of the bond between OW and C atoms is a must. The C···OW distance in the 100-ps simulation varies between 2.6 and 3.2 Å. The [CA2 N C O] improper angle minimal and maximal values are 2° and 34° with an average value of only 16 ± 5°. These values, however, are far too small for the formation of any of the tetrahedral intermediate.

As mentioned before, the crystal structure with the MVT-101 inhibitor includes a water molecule (301) that links the inhibitor with the flaps via two pairs of hydrogen bonds. It is assumed that this water molecule helps bind the substrate in a conformation conductive to cleavage, and is important for maintaining the structure of such complex as well as the transition state but does not participate in the reaction directly (Chatfield and Brooks, 1995). During the 100-ps dynamics its hydrogen bonds with the enzyme atoms are conserved. The minimal, maximal, and mean values with standard deviations of the distances between atoms hydrogen bonded to this water molecule are shown in Table 4.

The changes of the Asp-125:OD1-Asp-25:OD1 "inner" oxygens distance were also traced in the simulation. This distance in all known crystal structures, native and complexed with different inhibitors, varies between 2.6 and 3 Å. For example in the 3HVP native structure it equals 3.0 Å (Wlodawer et al., 1989), in the complex with MVT-101 inhibitor (Miller et al., 1997) it reaches 2.7 Å. In the studied case, the minimal distance between the "inner" oxygens is 2.93 Å, and the maximal one 4.26 Å with the average value of 3.47 ± 0.18 Å over the 100-ps trajectory.

The short distance between the "inner" oxygens indicates that if these atoms do not share a proton (as in the studied case), they should form a hydrogen bond with a lytic water molecule 207. Such hydrogen bond was found between Asp-125:OD1 and HW1 of the lytic water with a mean distance of 3.47 ± 0.21 Å. However, this water molecule experiences large mobility also forming and breaking the hydrogen bonds with the other atom of aspartyl side chain (Asp-25:OD1).

Also the planarity of the active site aspartates was reported. Most probably, the almost coplanar conformation of the dyad is crucial for the enzymatic function and for the binding of a substrate or inhibitor. The minimal, maximal, and mean values of the [Asp-125:OD2 Asp-125:OD1 Asp-25:OD1 Asp-25:OD2] dihedral angle are equal to 81, 30, and 21 ± 17°, respectively. Thus, on the average, the dyad keeps its planarity in the simulations even though large fluctuations are observed. The Asp-25:C-Asp-125:C average distance in this 100-ps simulation was equal to 5.5 ± 0.2 Å.

A hydrogen bond was observed between Asp-125:OD2 and Met-204:N of the cleaved peptide bond. The average distance between these atoms in the 100-ps simulation equals 3.06 ± 0.25 Å with the values ranging from 2.30 to 4.85 Å. The CA2···CA3 carbon minimal, maximal, and mean distances of the peptide bond carbons (see Fig. 1) were found to be 3.4, 4.0, and 3.8 ± 0.1 Å, respectively.

The MD/AVB simulation was performed with eight valence structures (1-6, 13, and 14 of the first model presented in the section AVB description of the system). Their contributions, i.e., c, were analyzed and are presented below ("" indicates negligible contributions in the 100-ps simulation; for the fragment names see Fig. 1):

Steered MD/AVB simulations of the first steps of the reaction

The reactions occurring in peptides and proteins happen most often on a time scale much larger than 1 ns, i.e., the typical scale for present MD simulations. Taking into account that the nucleophilic attack by water is usually the slow step of the cleavage process, the 100 ps of MD/AVB simulation was, of course, not sufficient for the reaction to occur spontaneously. Therefore, for an improved sampling of the configurational space a time-saving technique was applied, which is based on a deformation of the potential energy surface in such a manner that the system is steered to cross the energy barrier. Hence, in the HIV-1 PR case, to gain the nucleophilic attack of the hydroxy anion in a shorter time scale than in reality, external forces were applied to the oxygen nucleus of the lytic water molecule (OW) and to the peptide carbon nucleus (C) (see Fig. 1). The distance between these nuclei served as a reaction coordinate (r). The force is a consequence of an external harmonic potential U = 1/2k[r  (r₀  vt)]² added to the Hamiltonian of the system (k being the stiffness of the harmonic spring; r₀, equilibrium distance; v, velocity with which the harmonic spring is pulled; t, time step). As a result the C and OW atoms were driven toward each other to an r₀  vt distance value decreasing in time. If the restraining potential changes slowly in comparison with the atomic motion, within the simulation period may be approximated as constant. One should be, however, aware of the fact that such force driven simulations do not address the question of the reaction trajectory and time scale. We are able to test the mechanism based on the assumption that in reality those two atoms come close to each other.

Several simulations were performed that differed in the number of valence structures, the starting geometry, the stiffness of the harmonic spring and the cutoff used for the long range interactions. We show below exemplar results of a 72-ps-steered MD/AVB simulation with eight selected valence structures. The structures are numbered as in section AVB Description of the System and colored as in Fig. 5:

1.	A1 + H2O + A2 + S  black line
2.	A1 + OH + H + A2 + S  red line
3.	A1 + OH + A2 + S  green line
4.	A1 + H2O + A2 + S  orange line
5.	A1 + OH + H + A2 + S  yellow line
6.	A1 + OH + A2 + S  blue line
13.	A1 + A2 + H + I1  turquoise line
14.	A1 + A2 + I1  violet line

View larger version (70K): [in this window] [in a new window]   FIGURE 5   Contributions of the probability amplitudes (c) during the force driven dynamics of the first step of the reaction (for colors see text).

The starting coordinates were chosen after 50 ps of molecular dynamics from the equilibrated ensemble. The initial distance between the C and OW atoms was equal to r₀ = r = 3.15 Å and it was changed up to 1.35 Å (length somewhat smaller than the equilibrium for the COH bond) with k = 2000 [kJ/mol × Å²] and v = 0.025 Å/ps. Hence, the equilibrium length of the harmonic spring was shortened slowly, 0.0000125 Å every time step. The nonbonded interactions were calculated with a twin-range method with the cutoff set to 25 and 30 Å. The relaxation time of the thermal bath was set to 0.1 ps. In this case the first proton transfer, HW2 onto Asp-125:OD2, occurred around the 52 ps when the C···OW distance was equal roughly 2.2 Å (Fig. 6). The distance between the donor, OW, and the acceptor, OD2, reached then its shortest value of ~2.4 Å, which grew after the transfer to 3.2 Å.

View larger version (54K): [in this window] [in a new window]   FIGURE 6   Proton transfer from the catalytic water molecule to Asp-125 in the force driven MD/AVB simulations with eight valence bond structures. Dark blue, green, and red lines correspond to OW···HW2, OD2···HW2, and OW···OD2 distances, respectively. Black line shows the equilibrium length of the harmonic spring, and the turquoise line shows the actual C···OW distance at each MD/AVB step.

The course of this step of the reaction and the corresponding changes in the electronic state of the system are reflected in the contributions of the structure coefficients (Fig. 5). It is worth noting that at the beginning the contributions of structures including water molecule in a covalent form (1st and 4th structure) are quite low, whereas those with a dissociated water molecule are much higher. This indicates that the surrounding environment (especially catalytic aspartates and the substrate) facilitate the dissociation of the lytic water molecule.

Before the proton transfer the contributions of the 2nd and especially the 5th structure are dominant in the AVB Hamiltonian. Both structures contain a dissociated water molecule. The 5th structure with all ionized and polarized forms of molecular fragments, just before the proton transfer, grows to its maximal value (its c reaches 0.7). Next, due to the proton transfer, the nucleophilic attack, and the growth of the OD2···OW distance, the coefficients of the 2nd structure become negligible, of the 5th structure drop to 0.35 and of the 13th valence structure gain significant contribution. Also, the structures 6th and 14th yield higher contributions.

To check if the obtained transition state geometry is metastable and the reaction does not go back immediately after the steered MD/AVB simulations of the nucleophilic attack a further run was carried out with no additional force and no constraints applied. The geometry with the C ... OW distance close to equilibrium, 1.5 Å was chosen. As a continuation of the previous steered MD/AVB, a 50-ps run was performed with no constraints and using the same combination of eight valence structures. During such simulation the system conserves the geometry of the transition state. The mean C-OW distance was 1.51 ± 0.05 Å, which is very close to the equilibrium distance of this bond in the RCOOHNHR (R==CH₃) molecular fragment obtained with DFT calculations (see Trylska et al. (2001)). The OD2···HW2, OW···HW2, OW···OD2 mean distances were 1.06 ± 0.03 Å, 2.3 ± 0.1 Å, and 3.2 ± 0.1 Å, respectively, showing that after the transfer the HW2 proton was still bonded to the Asp-125:OD2 atom. The main contribution belongs to the 5th and 13th structure, their c mean values were 0.43 ± 0.02 and 0.45 ± 0.03, respectively. It is worth noting that the contributions of the first four structures, describing the state of the system before this reaction step were negligible.

To check if the proton transfer and the following nucleophilic attack did not occur accidentally, another similar force driven simulation, 66 ps, was performed from a different starting geometry (taken after 42.3 ps of the simulation from the equilibrated ensemble). The C···OW starting distance was equal to 2.95 Å. The twin-range cutoff was set to 30 and 35 Å. The character of Figs. 5 and 6 and the geometric features of the transition state at equilibrium distance of COW (~1.5 Å) were similar.

Also, as a test case, calculations with various force constants were performed (2000, 1500, and 1000 [kJ/mol × Å²]) for two starting geometries described earlier. To save the computer time, the twin-range cutoff was decreased to 12 and 14 Å. In all cases the first reaction step occurred and the characteristic features of Figs. 5 and 6 did not change.

After the nucleophilic attack, the I1 molecular fragment is formed with an sp³ hybridization on the peptide carbon. The CA2···CA3 (see Fig. 1) distance during and after the nucleophilic attack is shortened, in the 72-ps simulation, from an average value of 3.8 ± 0.1 Å to an average of 3.2 ± 0.1 Å. The [CA2 N C O] dihedral angle reaches ~45°, and the planar angles around the C and N atoms gain the sp³-like values.

In literature concerning the reaction mechanism a rotation or conformational transition of the CN bond was postulated allowing the nitrogen lone pair to form a hydrogen bond with the protonated Asp-125:OD2 atom. The distances between the N···OD2 and N···HW2 atoms during the 72-ps simulation are shown in Fig. 7. It may be clearly seen that the hydrogen bond is formed around the 54th-ps of the MD/AVB simulation after the proton transfer. Fig. 8 shows the characteristic changes in the geometry of the active site atoms during this simulation.

View larger version (79K): [in this window] [in a new window]   FIGURE 7   Formation of the hydrogen bond between the peptide nitrogen and the OD2 atom of Asp-125. Black and gray lines correspond to the N···OD2 and N···HW2 distances, respectively.

View larger version (18K): [in this window] [in a new window]   FIGURE 8   Four frames from the steered MD/AVB simulations of the HW2 proton transfer and the nucleophilic attack: (1a) one of the starting geometries; (1b) position of the HW2 proton, midway between the OD2 and OW atoms; (1c) process of bond formation between the HW2 and OD2 atoms and between the hydroxy anion and the peptide carbon; (1d) the geometry of the transition state just after the nucleophilic attack.

The Asp-25:OD2-Met-203:O minimal and maximal distances were 3.0 and 5.4 Å, respectively, and the mean distance was 4.3 ± 0.3 Å. The Met-203:O and Asp-25:H distance was in the range of 2.1 to 5.5 Å and an average value of 3.5 ± 0.4 Å. The hydrogen bond is therefore not stable and is created and broken during this simulation.

The above calculations were performed with eight valence structures corresponding to the first steps of the reaction up to the formation of the tetrahedral intermediate. To analyze further steps leading to the weakening of the peptide bond, more valence bond structures are needed. We used 21 structures, 19 (nr 1-19) of the set presented for the first reaction model and two more (17' and 18') to check the mechanism in detail (see section AVB Description of the System).

Another similar force-driven simulation, 66 ps, was performed with the starting geometry taken from the previous simulation with the C···OW starting distance of 2.95 Å and a cutoff of 12 and 14 Å. No other constraints were applied. The changes of the distances for the proton transfer and the nucleophilic attack are shown in Fig. 9 for the time steps of interest. The contributions of the structure coefficients (c) are presented in Fig. 10. The inclusion of additional valence structures allows for a deeper insight into the possible mechanisms of the reaction. In the beginning of the process the 11th structure gives main contribution. It consists of all possible ionic fragments that appear before the reaction. Just before the proton transfer, its c reaches the maximal value of around 0.6. The contribution of the 5th structure, which was dominant in the previous simulation with eight valence structures, is much smaller. However, the 11th structure differs only in the description of the Asp-25 fragment allowing for dissociation of the acid protein H.

View larger version (50K): [in this window] [in a new window]   FIGURE 9   Fragment of the dynamics of the HW2 proton transfer process in the steered MD/AVB simulations with 21 valence bond structures. Dark blue, green, and red lines present the OW···HW2, OD2···HW2, and OW···OD2 distances, respectively. The black line shows the equilibrium length of the harmonic spring and the turquoise line the actual C···OW distance at each MD step.

View larger version (67K): [in this window] [in a new window]   FIGURE 10   Contributions of the valence bond structures (c) in the fragment of the force-driven MD/AVB simulation. The legend shows the numbers (see text) of the significant structures and their colors.

Other structures present are 2nd, 5th, 8th, and 10th, which also contain at least one of the ionic molecular fragments. After the proton transfer their contributions drop and the c values of the structures 13th, 17'th, and especially 15th increase. The 15th structure includes the just formed I1 intermediate molecular fragment and two ionized aspartates. Also the 13th structure includes the I1 molecular fragment and one of the aspartates ionized. Describing the system with 21 valence structures, after the formation of the I1 intermediate another two proton transfers are possible (see Figs. 2 and 3), one onto Met-203:O and the other onto Met-204:N of the peptide bond. It is unknown which molecular structure is formed first, I2 or I2', or if both are. The high contribution of the 17'th structure suggests that I2' molecular fragment is formed. The contributions of structures that include the I2 molecular fragment are negligible. Comparing this simulation with the previous one with only eight valence structures included, one may notice that also the structure that contains all, either ionized or polarized molecular fragments, dominates. If other molecular fragments are included (especially I2'), Asp-125 is not necessarily forming a strong covalent bond with the HW2 proton. The contributions of valence structures containing the protonated aspartate (A2) are not noticeable in the simulation. This aspartic acid only accepts the HW2 atom for a short time and passes it onto another proton acceptor in this case the N atom.

MD/AVB simulations of proton transfer from Asp-125 onto peptide nitrogen

After the formation of the I1 molecular fragment there may be a few possible paths that lead to the breakage of the C-N bond. Some of them are accounted for in this study (see Figs. 2 and 3). First, a test was performed to check which molecular intermediate is formed after I1 and before I3. This is a continuation of a test performed in the steered MD/AVB simulation with 21 structures "switched on." There are three possibilities: either a proton transfer from Asp-25 onto Met-203:O may occur forming the I2 molecular fragment or a proton transfer from Asp-125 onto Met-204:N forming I2', or both may happen simultaneously. The inclusion of molecular fragments: A1, A1, H, A2, A2, H, I1, I2, I2', and I3 (see Figs. 2 and 3) results in nine valence structures listed below. The most contributing ones in this simulation are in bold face.

13.	A1 + A2 + H + I1
14.	A1 + A2 + I1

15.  A2 + H + H + A1 + I1

16.  A2 + H + A1