Theoretical investigation of the interactions in binding pocket of Reverse Transcriptase.

Interactions in proteins have been studied using several chemical information techniques including quantum chemical methods that are applied to truncated systems composed of the ligand molecule and the surrounding amino acids of the receptor. In this work we adopt an approach to study these interactions accounting for as many as possible explicit solvent molecules and without the need of a fragmented calculation. Furthermore, we embed our quantum chemical calculations within a molecular dynamics framework that enables a fundamentally fast system for quantum molecular dynamic simulations (QCMD). Central to this new system for QCMD is the tight binding QC system, newly developed in our laboratories, and which combined with the MD paradigm results in an ultra accelerated QCMD method for protein-ligand interaction evaluations. We have applied our newly developed method to the Nevirapine (NVP)-Reverse Transcriptase (RT) system. We show how the proposed method leads us to new findings. The advanced QCMD was applied to a system of RT with NVP and it has led to the knowledge of specific groups and atoms that interact with surrounding amino acids of RT and help in drug binding. The information derived from this calculation may be used in designing drugs for NVP resistant virus strains that have binding capability like NVP.

Introduction

HIV-1 Reverse Transcriptase (HIV-1 RT) is an important and extensively studied antiviral target for the chemotherapy of AIDS because of its key role in virus replication. The inhibitors of HIV-1 RT can be divided into two main classes, nucleoside RT inhibitors (NRTIs) and non-nucleoside RT inhibitors (NNRTIs) (Rizzo et al., 2001, Mitsuya et al., 1990, Katz and Skala, 1994, Wang et al., 2001). The NNRTIs analogs such as Nevirapine, TIBO, HEPT, and Efavirenz are non-competitive inhibitors that lock the polymerase active site in an inactive conformation (Archer et al., 2001, De Clercq, 1996, Jonckheere et al., 2000). Although NNRTIs are highly specific and less toxic than nucleoside inhibitors, their therapeutic effectiveness is limited by relatively rapid emergence of drug-resistant HIV-1 strains. NNRTIs include structurally unrelated subclasses of compounds that bind to a common allosteric site, adjacent to the NRTI binding site, by a similar three-dimensional arrangement (Archer et al., 2001, De Clercq, 1994, Hannongbua et al., 2001). Nevirapine, a first generation NNRTI, has been approved for clinical use for the therapy of AIDS. Mutations of amino acid residues in the binding pocket of RT reduce sensitivity significantly, and, therefore, new potent drugs have been widely developed.

However, despite intensive experimental investigations, the detailed origin of enzyme–inhibitor interaction, a question of crucial importance, remains to be unclear. A better understanding is vital for the analysis of activities of mutant or designed proteins, and for the design of inhibitors as pharmaceutical lead compounds. Therefore, theoretical investigation has been an alternative method for studies of the enzyme–inhibitor interaction in detail. However, such investigation of larger molecular system is limited by the computational effort required and the accuracy of the method used. Recently, accurate molecular modeling for larger molecules, such as those in molecular biology, became more feasible due to new developments in computational chemistry (Maseras and Morokuma, 1995, Humbel et al., 1996, Svensson et al., 1996).

HIV-1 Reverse Transcriptase (RT) is an important target for drugs used in the treatment of AIDS. Drugs known as non-nucleoside RT inhibitors (NNRTI) appear to alter the structural and dynamical properties of RT, which in turn inhibit RT’s ability to transcribe. Molecular dynamics (MD), principal component analysis (PCA), and binding free energy simulations are employed to explore the dynamics of RT and its interaction with the bound NNRTI Nevirapine, for both wild-type and mutant (V106A, Y181C, Y188C) RT. These three mutations commonly arise in the presence of Nevirapine and result in resistance to the drug. The mutations cause a loss of van der Waals interactions between the drug and the binding pocket. This suggests that a good inhibitor should efficiently enter and maximally occupy the binding pocket, thereby interacting effectively with the amino acids around the binding pocket (Zhou et al., 2005).

Method

Periodic density functional theory method

Periodic DFT calculations were performed by using the DMol3 software (DMol3, version4.0, Accelrys, 1999). In order to reduce the computational cost, all core electrons were represented by effective core pseudopotentials. In this study, double numerical sets with polarization base were employed. Vosk–Wilk–Nusair (VWN) (Vosko et al., 1980) local correlation functional was used to optimize geometries. Generalized gradient approximation (GGA) in terms of Perdew–Wang exchange and correlation functional (Perdew et al., 1992) was used to evaluate energies. The orbital population and bond populations were obtained using truncated models by the DFT method that employs Cambridge serial total energy package (CASTEP). In CASTEP calculations, we used GGA in terms of Perdew–Burke–Ernzerhof functional (PBE) (Perdew et al., 1996). The crystal lattice size used for DFT calculations on truncated Nevirapine-Tyr188 acid complex was 20 Å × 20 Å × 20 Å.

Tight-binding quantum chemical molecular dynamics method

Ultra-accelerated quantum chemical molecular dynamics (UA-QCMD) simulator

Our novel UA-QCMD simulator consists of two parts. The first part is the tight-binding QCMD (TB-QCMD) simulator. In this simulator, an electronic structure calculation is performed by solving the Schrödinger equation (; , , ε, and refer to the Hamiltonian matrix, eigenvectors, eigenvalues, and overlap integral matrix, respectively) with the diagonalization condition (= ; refers to the unit matrix). In order to determine the off-diagonal elements of , and H, the corrected distance-dependent Wolfsberg-Helmholtz formula (Calzaferri et al., 1989) (Eq. (1)) was used.

In order to solve the Schrödinger equation in this simulator, parameters for Hamiltonian matrix are used, which will be explained later. For electronic structure calculations using our TB-QCMD simulator, the total energy of a system is obtained by using the following equation,where the first, second, and third terms on the right-hand side refer to the molecular orbital (MO) energy, columbic energy, and exchange-repulsion energy, respectively. The first term on the right-hand side of Eq. (2) is rewritten as follows,where the first and second term on the right-hand side refer to the mono-atomic contribution to the binding energy and the diatomic contribution to the binding energy, respectively (n is the number of electrons occupied in kth molecular orbital). A binding energy calculated from the second term of Eq. (3) is used for the determination of the DAB parameter in the Morse-type 2-body interatomic potential function described as Eq. (4),where EAB, DAB, βAB, RAB, and RAB refer to the interatomic potential energy between atoms A and B, binding energy between atoms A and B, factor for potential curve, interatomic distance between atoms A and B, and equilibrium interatomic distance between atoms A and B, respectively. The second part of our UA-QCMD simulator is the classical MD simulator. In our MD simulator, a Verlet algorithm is employed to integrate equation of motion.

Quantum chemical calculations were carried out using Colors program, which is based on our original tight-binding approximation (Kubo et al., 2004, Elanany et al., 2003, Luo et al., 2003, Jung et al., 2003). In the “Colors” program, high efficiency of computation is realized by adopting parameters in the Hamiltonian. In our program, the parameters are determined on the basis of DFT calculation results. This can realize both high calculation speed and high accuracy (Elanany et al., 2002, Elanany et al., 2003, Sasata et al., 2003, Jung et al., 2003). The crystal lattice size used for ‘Colors’ calculations on truncated Nevirapine-Tyr188 complex was 20 Å × 20 Å × 20 Å.

We used a combination of MD simulation with the quantum chemical method that allows the computation of electronic properties of the system with conventional MD to simulate the dynamic behavior of the system when local properties such as charges and the free field vary. The combination of MD with QC enables the identification of the interacting ligand atoms and groups of the enzyme and characterization of the main terms leading to the interaction. The simulation was carried out by Ryudo program for 150,000 steps with an integration time of 0.1 fs. Prior to this simulation, the truncated system was subjected to a single point quantum chemical calculation by Colors program. The Morse potential between different pairs of atoms was used for MD simulation by New-Ryudo program. This combination of Colors and New-Ryudo was repeated for three times to get the stability and accuracy of results.

Results and discussion

We carried out quantum chemistry-based molecular dynamics simulation using “Colors” and “New-Ryudo” program to investigate the interactions of anti-HIV drug Nevirapine and Reverse Transcriptase enzyme.

Fig. 1 shows Nevirapine (shown as space fill model) and protein part falling under 10 Å radii from Nevirapine with water molecules. This model was used for our MD simulation. This model contains 11,520 atoms, which includes 2699 water molecules. The water molecules were added in a volume of 50 Å × 50 Å × 50 Å. This system was simulated for 15 ps.

Figure 1

Nevirapine (shown as space fill model) and protein part falling under 10 Å radii from Nevirapine with water molecules.

Table 1 shows the comparison of orbital population of NVP-TYR truncated system (Fig. 2) from DFT and colors calculation.

Table 1

Comparison of orbital population of Nevirapine and Tyrosine 188 system from DFT calculation and colors.

	Castep		Colors
	%s	%p	%s	%p
Hydrogen	100	0	100	0
Carbon	28.36	71.64	30.06	69.94
Nitrogen	27.21	72.79	24.93	75.07
Oxygen	26.98	73.02	23.21	76.79

Figure 2

Truncated model of NVP–Tyrosine complex on which colors and castep calculations were done for comparison of orbital population.

We found that there is a good agreement between %s and %p-orbital populations obtained with DFT and colors calculations. The bond population, bond energy and the distance shown in Table 2 indicate that there is hydrogen-bonding interaction between Nitrogen atom (N1) from the pyridine ring of Nevirapine and hydrogen (H1) from Tyrosine 188 residue from RT enzyme (as shown in Fig. 3). The bond population, bond energy and the distance shown in Table 3 indicate that there is hydrogen-bonding interaction between hydrogen atom (H2) from the cyclopropane ring of Nevirapine and Nitrogen (N2) from Tyrosine 181 residue from RT enzyme (as shown in Fig. 4). The bond population, bond energy and the distance shown in Table 4 indicate that there is hydrogen-bonding Interaction between Oxygen (O3) from Nevirapine and Hydrogen (H3) from Valine106 (as shown in Fig. 5).

Table 2

Bond population, bond energy and distance between Nitrogen atom from the pyridine ring of Nevirapine and hydrogen from Tyrosine 188 residue from RT enzyme.

Bond	Distance (Å)	Energy (kcal/mol)	BP
N1–H1	2.96	−4.6	0.02

Figure 3

Interaction between Nitrogen atom from the pyridine ring of Nevirapine and hydrogen from Tyrosine 188 residue from RT enzyme.

Table 3

Bond population, bond energy and distance between hydrogen atom from the cyclopropane ring of Nevirapine and Nitrogen from Tyrosine 181 residue from RT enzyme.

Bond	Distance (Å)	Energy (kcal/mol)	BP
N2–H2	2.99	−6.8	0.042

Figure 4

Interaction between hydrogen atom from the cyclopropane ring of Nevirapine and Nitrogen from Tyrosine 181 residue from RT enzyme.

Table 4

Bond population, bond energy and distance between Oxygen (O3) from Nevirapine and Hydrogen (H3) from Valine106.

Bond	Distance (Å)	Energy (kcal/mol)	BP
H3–O3	2.9	−4.9	0.02

Figure 5

Interactions between Oxygen from Nevirapine and Hydrogen from Valine106.

Conclusion

Interactions in protein have been studied using several chemical information techniques including quantum chemical methods that are applied to truncated systems composed of the ligand molecule and the surrounding amino acids of the receptor. In this work, we embed our quantum chemical calculations within a molecular dynamics framework that enables a fundamentally fast system for quantum molecular dynamic simulations.

In the present work we propose the use of ultra accelerated quantum chemical molecular dynamics (QCMD) to study the interactions between RT enzyme and anti-HIV drug molecule for the first time. Central to our methodology is the system for QC based on the tight binding methodology, encoded in the COLORS system developed in our laboratories. The proposed methodology enables treatment of a large molecular system including explicit solvent molecules, an elusive aspect to more orthodox quantum chemical calculations based on first principle calculations. The ultra accelerated QCMD was applied to a system of RT with Nevirapine, and it has led to new insights, not unveiled so far, of the atoms and groups of NVP drug and RT enzyme involved in hydrogen bonding interactions in the binding pocket of the enzyme.

We confirmed that tight binding quantum chemical molecular dynamics simulation program developed in our laboratory is an effective tool to investigate bio-molecular systems as demonstrated by this study on RT enzyme. Using our in-house programs, we found those atoms and groups that are actively involved in the interactions between NVP and RT. This information is useful in designing potent analogs of Nevirapine, which may act as substitute for NVP in case of NVP-resistant strains of Human Immunodeficiency Virus.