Interaction of HIV-1 aspartic protease with its inhibitor, by molecular dynamics and ab initio fragment molecular orbital method

For the three complex crystal structures of HIV-1 aspartic protease (an enzyme of AIDS) with its inhibitor in the Protein Data Bank, molecular dynamics of the generalized Born surface area and the ab initio fragment molecular orbital of an ABINIT-MP calculation was performed to obtain the binding free energy, the molecular orbital energy, the interaction energy of residues with an inhibitor and the charge transfer at the active site. The inhibitors are five symmetric cyclic ureas, of which three were modelled, and an asymmetric dipeptide. The interaction energy of the inhibitor at the active sites of aspartic acid is as great as 50 kcal mol⁻¹, coinciding with a tetrahedral transition state. For the inhibitor with a higher affinity, charge was transferred to the inhibitor from the active site. The difference in symmetry of the inhibitor was not evident. Binding free energy corresponds to the experimental value of the binding constant, while molecular orbital energy does not always, which is considered to be an entropy effect.

1. Introduction

A retrovirus, human immunodeficiency virus (HIV), is the etiology of acquired immunodeficiency syndrome (AIDS). HIV proliferates under its own protease. HIV-1 PR is an important target enzyme for the inhibition of HIV proliferation. HIV-1 PR consists of two chains, which constitute a twofold rotational C₂-symmetric homo-dimer. Each chain consists of 99 amino acid residues and has the characteristic disposition Asp–Thr–Gly of aspartic protease at positions 25–27. Fig. 1 shows the structure of HIV-1 PR complexed with a cyclic urea inhibitor XK2 [Protein Data Bank (PDB) ID 1hvr; Lam et al., 1994] viewed perpendicular to the C₂ axis. The aspartic acid of the active site hydrolyzes the peptide bond of the substrate catalytically via a tetrahedral transition state (Doi et al., 2004). To analyze the enzymatic reactions, the binding free energy of an inhibitor to the enzyme is given by the equation from the measured binding constant K_i, where R is the gas constant and T is the absolute temperature. It is also calculated by molecular dynamics taking the water effect into account to compare with the experimental value (Gohlke & Case, 2004). To find the charge transfer between the enzyme and an inhibitor, which seems important in enzymatic reactions, quantum mechanics has to be applied after structural optimization by classical mechanics. It takes a tremendous amount of time to complete quantum mechanical calculations of macromolecules such as proteins. Here, a new method of quantum mechanics for proteins, the fragment molecular orbital method (FMO) `ABINIT- MP', developed by one of the authors (Kitaura, Sawai et al., 1999; Kitaura, Ikeo et al., 1999; Nakano et al., 2000, 2002), is used for both binding energy and charge calculation.

Figure 1 Structure of HIV-1 PR complexed with a cyclic urea inhibitor XK2 (Accelrys, 2001). Light-blue residues are less hydrophobic, red slightly hydrophobic, and white Gly and others.

2. Structural data

Two types of HIV-1 PR complex sample were selected: a complex with a symmetric inhibitor corresponding to the C₂-symmetric homo-dimer of the PR, and a complex with an asymmetric inhibitor. Cyclic urea was selected as the symmetric inhibitor and peptide derivative as the asymmetric inhibitor. Of the X-ray diffraction structural data deposited in the PDB, two complexes with cyclic urea (1hvr and 1ajx; Hultén et al., 1997; denoted AH1) and one with a peptide derivative (1d4h; Andersson et al., 2003; denoted BEH) were selected. From 1hvr, three analogues were modelled: XK1, XK3 and XK4. The structural formulae of all six inhibitors, crystals and modelled, are shown in Fig. 2.

Figure 2 Structural formulae of selected inhibitors to HIV-1 PR.

3. Calculation

3.1. Molecular dynamics, `mm_gbsa'

Molecular dynamics calculations were performed using AMBER (Case et al., 2002), using the force field gaff `general Amber force field' and the generalized Born surface area model, mm_gbsa (Srinivasan et al., 1998), of type 2, which takes a much shorter time than the model that arranges water explicitly. Charges were given to an inhibitor by RESP (restrained electrostatic potential) using Antechamber, `an accessory software package for molecular mechanical calculations'. The electrostatic potential input for RESP was obtained by the quantum mechanical program GAUSSIAN (Frisch et al., 1998) at the HF/6-31G* level. For molecular dynamics, Sander (simulated annealing with NMR-derived energy restraints) was run for 50 ps at 300 K, during which time ten snapshots were sampled. For two complex crystals, 1ajx and 1d4h, Sander was minimized only and equilibrated over 5000 steps.

From the snapshots of the complex obtained above, receptors and ligands were extracted to obtain the free energy of the receptor, G_receptor, and that of the ligand, G_ligand. Binding free energy was calculated by the equation below,

The electrostatic contribution to solvation free energy is calculated by generalized Born (GB) methods. The hydrophobic contribution to the solvation free energy is determined with a term dependent on the solvent-accessible surface area (Sitkoff et al., 1994). The surface area is computed using the MOLSURF program, which is based on an idea primarily developed by Connolly (1983). The calculated was compared with the experimental value obtained from the experimental binding constant k_i,

3.2. Fragment molcular orbital, `ABINIT-MP'

Molecular orbital calculations were performed via the fragment molecular orbital method, using the message passing interface (MPI), parallel version `ABINIT-MP'. In this method a protein is divided into fragments by residues, at the default two residues per unit, as shown in Fig. 3. Treating each N fragments as a monomer, two fragments further are paired to form N(N − 1)/2 dimers. The total system is constituted of monomers and dimers. It is not necessary to treat all the system at once, and parallel runs are possible to speed up the calculation with only a small energy error. The interaction between fragments in enzymes can be analyzed. The complex, receptor and inhibitor obtained from the last snapshot of molecular dynamics were locally minimized and used for input of the ABINIT-MP molecular orbital to obtain the binding energy . The basis set used was 6-31G. Interactions between the inhibitor and the protease receptor were obtained from the checkpoint file of the output.

Figure 3 Dividing a protein into fragments, two residues per unit.

4. Results and discussion

4.1. Binding free energy and fragment molecular orbital energy

of the complexes with cyclic urea inhibitors [a complex crystal with XK2 (1hvr) and three modelled complexes, XK1, XK3 and XK4] obtained by molecular dynamics and obtained by ABINIT-MP are shown in Table 1(a). The values of the XK series correspond fairly well to the experimental value - RTlnK_i, while the values do not correspond so well. The results for another complex crystal with a cyclic urea, 1ajx, and an asymmetric dipeptide 1d4h inhibitor, which have the largest and smallest binding constant, respectively, are given in Table 1(b). corresponds to the experimental value, while is reversed. E_receptor of both crystals has almost the same value 77491.0 a.u. The difference between and is considered to be entropy and solvent effect (Raha & Merz, 2004). seems to correspond to the quantum mechanical part of .

where ΔH and ΔS are the enthalpy and entropy change, respectively, LJ is the Lenard–Jones potential, ΔSA_C,N,O,S is the solvent entropy based on the surface area burial C, N, O and S atoms, num(rot bonds) is the number of rotational bods in the ligands, and sub- and superscripts b, solv, f and g denote binding, solvation, heat of formation and gas phase, respectively. The dispersive part of nonpolar interaction (1/R⁶)LJ is calculated using the attractive part of the Lennard–Jones potential. It has been found by isothermal titration calorimetry of the binding that an inhibitor with high affinity is strongly exothermic (favorable enthalpy change, −7.6 kcal mol⁻¹) and has a more balanced distribution of enthalpic and entropic interactions (Velazquez-Campoy et al., 2001). The entropic term should be calculated by normal mode analysis in molecular dynamics.

Table 1
Binding free energy and fragment molecular orbital energy (kcal mol⁻¹) of HIV-1 PR with its inhibitor compared with experimental values: (a) a complex crystal with cyclic urea (XK2) and three (XK1, XK3, XK4) modelled from XK2 (1hvr); (b) two complex crystals with cyclic urea (AH1; 1ajx) and a dipeptide BEH (1d4h)

Inhibitor	Ki (nM)	- RTlnKi
(a)
XK1	4.7	−11.43	−32.37	−39.34
XK2	0.31	−13.05	−57.03	−39.85
XK3	2.14	−11.94	−37.28	−36.96
XK4	0.27	−13.13	−46.56	−49.51
(b)
AH1	12.2	−10.86	−46.22	−70.01
BEH	0.10	−13.73	−54.38	−47.00

4.2. Interaction of inhibitors with HIV-1 protease residues

As stated in the Introduction, in HIV-1 PR, Asp25 and Asp124 at the active sites hydrolyze a peptide bond of the substrate catalytically via not a triangular but a tetrahedral transition state bearing a negative charge like other PRs (Ser PR for example; Branden & Tooz, 1999), as shown in Fig. 4 (Doi et al., 2004). ABINIT-MP is able to calculate not only the total binding energy but also the interaction between residues of the receptor and the inhibitor, and the charge transfer from the receptor to the inhibitor at the active site, , where the complex consists of a receptor and an inhibitor (Nakano & Kato, 2004). The electric charge was calculated using Mulliken's method. The interaction energies of the inhibitors and the protease are shown in Fig. 5. The interactions at the active sites, Asp25 and Asp124, are as great as 50 kcal mol⁻¹, corresponding to the tetrahedral transition state. The interactions are not necessarily balanced at both sites in symmetrical cyclic urea inhibitors (XK series and AH1), but conversely the lower the binding constant, the more balanced the interaction. Hydrogen bonds are formed between the inhibitors, AH1 and BEH, and the active sites, Asp25 and Asp124, as shown in Fig. 6, corresponding to the transition state. In the symmetric cyclic urea inhibitor AH1, the two hydrogen-bond lengths between the hydroxy group of AH1 and sites Asp25 and Asp124 are almost equal (1.961 and 1.870 Å, respectively), while in the asymmetric dipeptide BEH, the length difference is greater (1.928 and 2.559 Å). The charge transferred from the receptor to the inhibitor, , at the active sites of the two complex crystals is shown in Table 2. At the active site Asp25, the charge transfer from receptor to inhibitor of high-affinity BEH is larger than that of low-affinity AH1, while at the other site, Asp124, the direction of the charge transfer is reversed and the value is low.

Table 2
Charge transfer from receptor to inhibitor, , at the active site (in a.u.)

Active site	Ligand	qcomplex	qreceptor
Asp25	AH1	−0.7354	−0.7558	0.0204
	BEH	−0.7445	−0.8009	0.0564
Asp124	AH1	−0.77786	−0.77473	−0.0031
	BEH	−0.82909	−0.81347	−0.0156

Figure 4 Substrate tetrahedral transition state of HIV-1 PR with negative charge.

Figure 5 Interacton of HIV-1 PR residues and inhibitors XK1–BEH. Energy in kcal mol^−1; numbering of residues is through the dimer; each has 99 residues.

Figure 6 Hydrogen bonds between inhibitor and active sites. PR protein is drawn in Cα, inhibitor, active sites Asp25 and Asp124 in wire frame. C atoms are colored green, O atoms red, N atoms blue and H atoms white. The inhibitor is (a) a cyclic urea (AH1) and (b) a dipeptide (BEH).

Footnotes

‡Present address: Koyano Research Laboratory, Midorigaoka 2-2-9, Meguro-ku, Tokyo, 152-0034, Japan.