diffraction structural biology Journal of Synchrotron Radiation

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 À1 , 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.


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 2 -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 2 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 ÁG ¼ ÀRT ln K i 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., 2000Nakano et al., , 2002, is used for both binding energy and charge calculation.

Structural data
Two types of HIV-1 PR complex sample were selected: a complex with a symmetric inhibitor corresponding to the C 2 -symmetric homo-

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. 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. 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 ÁG 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 ÁG was compared with the experimental value ÁG exp obtained from the experimental binding constant k i ,

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 ÁE. The basis set used was 6-31G. Interactions between the inhibitor and the protease receptor were obtained from the checkpoint file of the output.

Results and discussion
4.1. Binding free energy and fragment molecular orbital energy ÁG 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 ÁE obtained by ABINIT-MP are shown in Table 1  Structural formulae of selected inhibitors to HIV-1 PR.

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

Figure 4
Substrate tetrahedral transition state of HIV-1 PR with negative charge. which have the largest and smallest binding constant, respectively, are given in Table 1 (b). ÁG corresponds to the experimental value, while ÁE is reversed. E receptor of both crystals has almost the same value 77491.0 a.u. The difference between ÁG and ÁE is considered to be entropy and solvent effect (Raha & Merz, 2004). ÁE seems to correspond to the quantum mechanical part of ÁH f .
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 6 Þ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 À1 ) 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.

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 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. between residues of the receptor and the inhibitor, and the charge transfer from the receptor to the inhibitor at the active site, Áq ¼ q complex À q receptor , 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 À1 , 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 hydrogenbond 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, Áq, 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.