Insights into Saquinavir Resistance in the G48V HIV-1 Protease: Quantum Calculations and Molecular Dynamic Simulations

INTRODUCTION

The spread of acquired immune deficiency syndrome (AIDS) has constantly threatened the world because the disease leads to a significant loss of morbidity and mortality. Unfortunately, chemotherapy for the disease has, in many cases, failed to achieve complete viral suppression (Deeks, 2003). This relies on the fact that the human immunodeficiency virus (HIV) develops resistance to antiretroviral drugs by genetic mutation. The development of novel drug targets and HIV vaccine is promising but the results of those studies remain far from the clinical stage. Understanding the mutations that confer resistance to available drugs is thus an urgent issue in HIV chemotherapy.

The HIV type-1 protease (HIV-1 PR) is an important target for AIDS chemotherapy. This viral protein cleaves the gag and pol nonfunctional polypeptide into functional proteins essential for maturation of infectious HIV particles (Debouck et al., 1987). The protein is a homodimer. Each protein monomer consists of 99 amino acids (Meek et al., 1989). From x-ray data (Fig. 1 B), the substrate/inhibitor binding site is located at the dimer interface (Hong et al., 1996, 2000; Jaskolski et al., 1991; Krohn et al., 1991; Swain et al., 1990; Vondrasek and Wlodawer, 2002). As a member of the aspartyl protease family, HIV-1 PR is composed of the conserved sequences, so-called the binding triads: D25-T26-G27 and D25'-T26'-G27', of which D25 and D25' are known as the active site residues. These two ionizable residues play a major role in the catalytic reaction.

Because of the therapeuticalIy important enzyme, structural and functional studies have been carried out to gain understanding of molecular mechanisms of the proteolytic cleavage process (Hyland et al., 1991; Northrop, 2001; Okimoto et al., 1999; Scott and Schiffer, 2000; Smith et al., 1996). The size and the availability of high-resolution x-ray structures of HIV-1 PR are amenable for molecular dynamics (MD) technique to investigate the relationship of structure, dynamics, and function of the enzyme (Collins et al., 1995; Harte et al., 1992, 1990; Levy and Caflisch, 2003; Piana et al., 2002; Scott and Schiffer, 2000) as well as to serve as a test system for developing computational methodology (Piana et al., 2004, 2001; Wang and Kollman, 2000; York et al., 1993a). The MD approach has provided insightful information on the enzyme-substrate interactions and binding conformations, the protonation states of the active site residues, the role of the flexible flap and the binding waters, and drug resistance. Characterizations of structural intermediates have been made useful for rational drug design (Randolph and DeGoey, 2004; Roberts et al., 1990; Rodriguez-Barrios and Gago, 2004).

Saquinavir (SQV, Fig. 1 A), a peptidomimetic protease inhibitor, is clinically used to treat infected HIV patients. The inhibitor containing a nonhydrolyzable hydroxyethylene isostere was designed based on the transition state structure in the enzyme-substrate complex. Combination of PR and reverse transcriptase inhibitors appears to be a highly effective treatment against HIV (Boucher, 1996). The PR inhibitor blocks the maturation step of the HIV life cycle, which is the crucial stage in the formation of new viral particles. Nevertheless, the current cure with SQV has introduced several resistant variants of HIV-1 PR, some of which can dramatically reduce drug susceptibility (Vondrasek and Wlodawer, 2002). G48V and L90M are considered as the primary mutations commonly occurring in vivo or in vitro (Eberle et al., 1995; Vaillancourt et al., 1999). These "signature" residue mutations can be associated with a dramatic decrease in drug susceptibility. According to K^sub i^ values, the G48V, L90M, and G48V/L90M mutants decrease saquinavir sensitivity by 13.5-, 3-, and 419-fold with respect to that of the wild-type (wt) protease (Ermolieff et al., 1997).

Among common mutations associated with antiretroviral drug resistance, G48V is a unique mutation characteristically generated by SQV. In a view of substituted-type residue, glycine was replaced by a bulkier side-chain residue. The sleric conflict of the mutant should involve in a destabilization of the complex. Although several x-ray structures of the HIV-1 PR provided valuable information on the inhibitor binding, this is not the case for the primary resistance to SQV. The crystal structure of the G48V complex is not yet available. With the aid of the available x-ray data, the molecular modeling techniques offer an opportunity to investigate the structural basis of the mutant enzyme (Prabu-Jeyabalan et al., 2003; Swain et al., 1990).

The missing hydrogens in the structural data have led to studies of the ionization state of the active site residues D25/D25' (Smith et al., 1996; Wang et al., 1996; Wlodawer and Vondrasek, 1998; Yamazaki et al., 1994). This subject is important for drug design in a way to optimize the interactions of the inhibitor with the enzyme. Different protonation models were found depending upon the local environment of the enzyme-inhibitor complex. The single protonation at one of the two acidic residues has been most commonly observed with the binding of the hydroxyl-ethylene-based inhibitors (Baldwin et al., 1995; Chen and Tropsha, 1995; Hyland et al., 1991; Smith et al., 1996; Wang et al., 1996; Wlodawer and Vondrasek, 1998). From NMR experiments, the neutral D25/D25' side chain (diprotonation) was determined in the presence of inhibitor diol groups (Yamazaki et al., 1994), whereas the dianionic form (unprotonation) was observed in the free enzyme (Smith et al., 1996; Wang and Kollman, 2000).

In this study, we employed a computational approach to access information regarding molecular structure and dynamics of the G48V HIV-1 protease conferring to saquinavir resistance. The MD simulations were carried out for the wt and the G48V HIV-1 protease complexed with saquinavir in explicit aqueous solution. The study of the protonation state of the HlV-1 PR complexed with SQV has been carried out before exploring the structure and dynamic data of the signature resistance. The density functional theory (DFT), ONIOM, and molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) methods have been performed to identify the protonation model of the active site residues. The MM/PBSA approach offers an efficient computation for calculating the binding free energy of biomolecules (Kollman el al., 2000; Srinivasan et al., 1998; Wang and Kollman, 2000). The method has been extensively used to study protein-ligand complexes. The quantum-based approach, DFT and ONIOM, has been useful in providing accurate energy information of the interested region. In particular, the hybrid quantum mechanical/molecular mechanical (QM/ MM) method, ONIOM (our own N-layered integrated molecular orbital and molecular mechanics), developed by Morokuma, has been extended from small molecules to biological applications (Friesner and Beachy, 1998; Morokuma, 2002; Prabhakar et al., 2004; Torrent et al., 2002). Its efficiency has been improved over the years. Simply, the concept of the ONIOM approach is partitioning a large molecular system into onion skin-like layers, and applying the quantum mechanics and molecular mechanics methods to the defined different parts (Morokuma, 2002). In the partitioned system, the high-level quantum computations engage the essential part of the central activity, whereas the lower-level energy calculations take into account the contribution of the remaining region. The comparison of MD results of the two systems provides insightful details of how the G48V mutant is associated with saquinavir resistance. The study provided fundamental principles on the molecular mechanism of inhibitor binding and resistance, which will be useful for designing an anti-HIV inhibitor to combat AIDS.

METHODS

Before starting the MD simulations, the problem of the protonated state of the active site residues was addressed, since such information cannot be directly obtained from the x-ray data. Our approach consists of molecular orbital energy calculations and solvent continuum free-energy calculations. However, structural data of the complexes, particularly the structure of the G48V-SQV, are not available, and the quantum-based computation for the whole enzyme-inhibitor complex is not feasible. A strategy is developed. First, MD simulations of the four protonation states for the wt and for the G48V complexes were performed to obtain the protonation models, which were subsequently subjected to calculate interaction and binding energies. Details of the methodology are described as follows.

Initial structure

The x-ray structure of the wt HIV-1 PR complexed with Ro 31-8959, saquinavir, (Protein Data Bank code 1HXB; 2.3 [Angstrom] resolution) was used as a starting model. All missing hydrogens of the protein were added using the LEaP module in the AMBER 7 software package (case et al., 2002). The protonation state of the ionizable residues, the C- and the N-termini, except for D25/25', was assigned based on the predicted pKa values at pH 7. The pKas of ionizable residues were calculated based on the Poisson-Boltzmann free-energy calculations (see the pKa prediction). The results concluded that all Lys, Arg, Glu, Asp, and the terminal groups are charged, whereas His was in the neutral form. Protonation states of D25/25' were explicitly assigned into four different ionizable states, including unprotonation, monoprotonation (each site of D25 and D25'), and diprotonation (protonated at both aspartyl residues). For the wt study, the simulated systems were labeled as wt-unpro, wt-mono25, wt-mono25', and wt-dipro, respectively.

The starting structure and force-field parameters for the inhibitor were obtained as follows. Hydrogens were added to the x-ray coordinates of SQV (1HXB) by taking into account the hybridization of the covalent bonds. Geometric optimization was subsequently performed at the Hartree-Fock level with 6-31G** basis functions to adjust the bond-length involving hydrogens. Then, the RESP fitting procedure was employed to calculate partial atomic charges of the inhibitor (Cornell et al., 1993). Force-field parameters of the inhibitor were assigned based on the atom types of the Cornell et al. (1995) force-field model. Gaussian 98 (Frisch et al., 2002) was used to optimize the molecular structure, generate electrostatic potentials, and calculate ab initio energies. Partial charge generation and assignment of the force field were performed using the Antechamber suite (Wang et al., 2001).

The preparation of the initial structure for the simulation of the G48V mutant-SQV complex was similar to that of the wt complex. The comparative model of the mutant was constructed based on 1HXB because the three-dimensional structure of the G48V-SQV complex is not available. It should be noted that the x-ray structure of the double mutant, G48V/L90M-SQV complex (1FB7) could be considered as an alternative template. However, the x-ray coordinates of the second monomer of the double mutant are not available. Thus, 1HXB is considered to be more appropriate as a template. The simulated systems of the mutant complex consist of four protonated states, which were defined similar to those of the wt complex, i.e., unprotonation (mt-unpro), monoprotonations (mt-mono25 and mt-mono25'), and diprotonation (mt-dipro).

The next step was to incorporate the solvent and counterions into the models previously prepared. The crystallographic waters were also included in the simulations. Each model was solvated with the TIP3P waters (Jorgensen et al., 1983) and neutralized by the counterions using the LEaP module. The total number of the TIP3P waters in the periodic box for all systems was in a range of 9100-9900 molecules.

Molecular dynamics simulations

Energy minimization and MD simulations were carried out using the SANDER module of AMBER 7 (Case et al., 2002) with the Cornell force field (Cornell et al., 1995). The whole systems were subjected to energy minimization within a range of 200-5000 steepest descent steps to avoid bad contacts. It should be noted that position-restrained minimizations of some particular regions were carried out for systems that clashed during minimization because of incidental overlay of atoms. This procedure was repeated until there was no sign of bad contacts. The resulting protein structure was compared with the before-minimized structure. Root mean-square displacement (RMSD) for nonhydrogen atoms of the compared protein structures showed that there were no RMSDs exceeding 0.1 [Angstrom] in all systems.

The MD simulation was performed employing the periodic boundary condition with the NPT ensemble. A Berendsen coupling time of 0.2 ps was used to maintain the temperature and pressure of the systems (Berendsen ct al., 1984). The SHAKE algorithm (Ryckaert et al., 1977) was employed to constrain all bonds involving hydrogens. The simulation time step of 2 fs was used. All MD simulations were run with a 12 [Angstrom] residue-based cutoff for nonbonded interactions and the particle-mesh Ewald method was used for an adequate treatment of long-range electrostatic interactions (York et al., 1993a).

The simulation consists of thermalization, equilibration, and production phases. Initially, the temperature of the system was gradually heated from 0 to 298 K during the first 60 ps. Then, the systems were maintained at 298 K until MD reached 400 ps of the simulation. Finally, the production phase started from 400 ps to 1 ns of the simulation. The convergence of energies, temperature, pressure, and global RMSD was used to verify the stability of the systems. The MD trajectory was collected every 0.1 ps. The 600 ps trajectory of the production phase was used to calculate the average structure. All MD simulations were carried for 1 ns. Analysis of all MD trajectories i.e., RMSD, distances, torsion angles, etc. was carried out using the CARNAL and Ptraj modules of AMBER 7. The geometry and stereochemistry of the protein structure were validated using PROCHECK (Laskowski et al., 1996). In summary, a total of eight systems for the MD simulations were carried out.

Graphic visualization and presentation of protein structures were done using RasMol, Swiss-Pdb Viewer (Guex and Peitsch, 1997), WebLab Viewer (Accelrys, San Diego, CA), and MolScript (Kraulis, 1991).

The pKa prediction