Theoretical studies on the molecular basis of HIV-1RT/NNRTIs interactions

Abstract

Molecular dynamics simulations (MD) of the human immunodeficiency virus type 1 reverse transcriptase (HIV-1 RT) complexed with the four non-nucleoside reverse transcriptase inhibitors (NNRTIs): efavirenz (EFV), emivirine (EMV), etravirine (ETV) and nevirapine (NVP), were performed to examine the structures, binding free energies and the importance of water molecules in the binding site. The binding free energy, calculated using molecular mechanics Poisson-Boltzmann surface area (MM-PBSA), was found to decrease in the following order: EFV ∼ ETV > EMV > NVP. The decrease in stability of the HIV-1 RT/NNRTI complexes is in good agreement with the experimentally derived half maximal inhibitory concentration (IC₅₀) values. The interaction energy of the protein-inhibitor complexes was found to be essentially associated with the cluster of seven hydrophobic residues, L100, V106, Y181, Y188, F227, W229 and P236, and two basic residues, K101 and K103. Moreover, these residues are considered to be the most frequently detected mutated amino acids during treatment by various NNRTIs and therefore, those most likely to have been selected in the population for resistance.

Keywords::

• Molecular dynamics simulation
• reverse transcriptase
• HIV-1
• NNRTIs

Introduction

Since human immunodeficiency virus type 1 (HIV-1), a causative agent of acquired immune deficiency syndrome (AIDS), was first discovered in 1981 [1], HIV-1 reverse transcriptase (RT) has been the subject of intensive studies. HIV-1 RT is part of the HIV capsid and has an essential role in the replication of the AIDS virus. This enzyme transcribes the single-stranded RNA of the AIDS retrovirus to a double-stranded DNA which can then be integrated into the host genome. The RT is a heterodimer composed of two subunits, the 66 kDa (p66, 560 residues) and the 51 kDa (p51, 440 residues) domains. Both subunits have a polymerase domain, but it is only functionally active in the p66 subunit, whereas the p51 subunit is simply a cleaved version of p66, and lacks the ribonuclease H (RNaseH) domain [2]. Due to its key role in the HIV life cycle, RT is an important target for antiviral agents in the treatment of AIDS.

Anti-HIV-1 RT drugs can be divided into two main classes: nucleoside reverse transcriptase inhibitors (NRTIs) and non-nucleoside reverse transcriptase inhibitors (NNRTIs). NRTIs are competitive to the nucleotide substrates and act as chain terminators when incorporated into viral DNA by HIV-1 RT. However, they instigate serious side effects, especially damage to mitochondria [3–5]. NNRTIs are non-competitive inhibitors which are highly specific for HIV-1 RT at a common allosteric site approximately 10 Å apart from the polymerase active site. In addition, they are less likely to cause adverse side effects by disruption of normal DNA polymerase activity [4–5]. Consequentially, interest is currently focused upon the NNRTIs group, aiming to understand and to provide a detailed insight into the HIV-1/NNRTIs interactions.

Residues that have been reported to play an important role in binding to the NNRTI consist of approximately 15 amino acid residues from the p66 subunit, L100, K101, K103, V106, T107, V108, V179, Y181, Y188, V189, G190, W227, W229, L234 and Y318 plus the E138 residue from p51. It appears that these binding pocket residues are flexible to some extent depending on the molecular shape and size, the specific chemical structure, and the binding mode of the individual NNRTIs [6–8]. The NNRTIs considered in this study are efavirenz (EFV), emivirine (EMV), etravirine (ETV) and nevirapine (NVP). Figure 1 shows their structures in the corresponding binding pockets of HIV-1 RT, their chemical compositions and numbering of selected atoms. NVP is classified as a member of the first generation of NNRTIs. EMV has been found to inhibit HIV strains that have developed resistance against NVP, albeit at a higher concentration than that required to inhibit the wild-type virus, but it failed in clinical trials [9]. EFV and ETV are classified as second generation inhibitors. EFV is the most prescribed NNRTI used for patients for first-line antiretroviral therapy because of its potent activity on wild-type HIV-1 [10]. ETV shows potent efficacy and retains potency in patients infected with NNRTI-resistant HIV-1 variants [7,11]. Experimentally derived half maximal inhibitory concentration (IC50) values of the wild-type HIV-1 RT treated with NNRTIs were reported in the following order: EFV > ETV > EMV > NVP [7,12–14]. However, more NNRTIs are needed to overcome the drug resistance mutations in the enzyme. Evidently, there exists no final treatment solution against the AIDS. NNRTI-resistant mutations cluster between position 98–108, 179–190, and 22–238 in the p66 subunit. K103N is probably the predominant mutation observed in patients receiving NNRTIs. Y188L has also been linked to high-level cross-resistance, mainly induced by NVP. In addition there are ten other residues with selected mutations, L100I, K101P, V106M/A, V108I, Y181C/I/V, G190A/S/C/E/Q, P225H, M230L, P236L and Y318F, have been observed [15].

Figure 1.  Geometries of EFV (A), EMV (B), ETV (C) and NVP (D), in their binding pockets of HIV-1 RT. Chemical structures of NNRTIs and the definition of atoms are also given.

This work aims to investigate the amino acid residues that are responsible for the HIV-1 RT–NNRTIs interactions by means of molecular dynamics and free energy calculation methods. Here, we constructed four model systems consisting of the HIV-1 RT with 1000 residue-long complexed with the four NNRTIs, (EFV, EMV, ETV and NVP). The size of the systems for the MD simulations was considerably large compared to what has been studied previously using the same protein with the same approach [9,13,16–18]. The calculation was focused upon examination of the structure, binding free energy and water accessibilities in the binding site of each of the four HIV-1 RT/NNRTI complexes. In addition, the distribution and binding of water molecules in the cavity of the enzyme-inhibitor complex, which is known to play an important role in the drug binding, was investigated and is discussed with reference to those found in the crystal structures. This information will be helpful for the rational design of new anti-HIV drugs with improved resistance profiles for anti-HIV therapy, as well as prediction of new resistance mutations.

Methods

Initial structures of HIV-1 RT and their complexes

The initial crystallographic structures of HIV-1 reverse transcriptase complexed with EFV, EMV, ETV and NVP were obtained from the Protein Data Bank (RCSB PDB, http://www.rcsb.org), entry codes 1FK9 [19], 1RT1 [20], 1SV5 [7] and 1VRT [19], respectively. The wild type system was reconstructed from the structure of the mutation K103N (1SV5) by replacing the glutamine by a lysine using the LEaP module in the AMBER 9 software package (UCSF, San Francisco, CA) [21]. The missing residues of 1FK9, 1RT1 and 1VRT were reconstructed with the help of another X-ray structure (PDB entry code 1SV5) and the subsequent low temperature annealing and energy minimisation procedures of the added residues whilst constraining the known X-ray residue positions. All the missing hydrogen atoms of the protein were added using the LEaP module in the AMBER 9 software package. The ionisation states of the amino acids with charged side chains were assigned using the PROPKA program (Boston, MA) [22].

Force field parameters of inhibitors

The starting structures and force field parameters for the four inhibitors were obtained in the following way: Ligand geometries were optimised using Gaussian03 [23] with the HF/6-31G* method. Then, single-point calculations were carried out to compute the electrostatic potentials around each compound using the same basis set and level of theory. The electrostatic potential was calculated by RESP [24]. Partial charge generation and assignment of the force field [25] were performed using the antechamber suite [26].

Molecular dynamics simulations

Energy minimisation and MD simulations were performed using the SANDER module of the AMBER 9 software package. An all-atom representation of the system was used employing the ff03 force field [27,28] to assign parameters for the standard amino acids. To incorporate the solvent and counter ions under consideration, each system was solvated using the TIP3P water model [29]. Neutralisation was performed by counter ions using the LEaP module. The total numbers of atoms were 142662, 141659, 141220 and 143341 for the RT/EFV, RT/EMV, RT/ETV and RT/NVP systems, respectively. The periodic boundary condition with the NPT ensemble was applied. The simulation steps consisted of thermalisation, equilibration and production phases. Initially, the temperature of the system was gradually increased from 0 K to 298 K during the first 50 ps. Then, the system was maintained at 298 K until the MD simulations reached 1.5 ns. Finally, the production phases were held up from 1.5 ns to 3.0 ns. A Berendsen coupling time of 0.2 ps was used to maintain the temperature and standard pressure of the system [30]. The SHAKE algorithm [31] was applied to constrain all the bonds involving hydrogen atoms, and a simulation time step of 2 fs was used. All MD simulations were run with a 10 Å residue-based cut-off for non-bonded interactions and the particle-mesh Ewald (PME) method was used as an adequate treatment of long-range electrostatic interactions [32]. The convergence of energies, temperature, pressure and root mean square displacements (RMSD) was monitored to verify the stability of the systems. After reaching an almost stable RMSD value, the production phase selected was adjusted from 1.5 ns to 3.0 ns for the four systems.

The binding free energy calculations

The binding free energies (ΔG_binding) of the four systems were estimated by the MM-PBSA approach [16,33]. In this method, the ΔG_binding value can be calculated according to: 1) the molecular mechanics interaction energy of the protein-inhibitor complex in vacuum, composed of the electrostatic interaction and the van der Waals interaction, 2) the contribution of solute entropy to binding and 3) the solvation free energy, composed of the polar and the nonpolar contributions. An extensive determination of the entropic effect from the normal mode analysis requires high computational demands and is therefore costly. However, the entropic contribution of the inhibitors may not be negligible for calculating the binding free energy, therefore, the binding free energies (ΔG_binding) incorporating the contribution of the entropy (TΔS) of the four systems were calculated. Owing to the system size (already 1,000 residues for HIV-1 RT), the entropic contributions were not taken into account for the whole system. Therefore, one hundred MD snapshots from 1500 to 3000 ps of the full-length RT complexes with a truncated HIV-1 RT protein composed of residues 1–350 were extracted for calculating ΔG_binding including the TΔS by using the NMode program [34]. These residues were selected because their location covers the binding pocket of HIV-1 RT. The solvation energy was calculated using the dielectric constant of 1 for the protein interior and 80 for the water molecule. The electrostatic solvent energy was calculated using the PBSA program [35], whilst the nonpolar contribution of the solvent was determined using the MolSurf program [36].

Results and discussion

Hydrogen bonds between inhibitors and the enzyme

To analyse the ligand-enzyme interaction, the occurrence of hydrogen bonds between the inhibitors and the binding pocket residues together with the percentage of occupations were determined based on the following criteria: 1) a proton donor-acceptor distance of ≤3.5 Å, and 2) a donor-H-acceptor bond angle of ≥120°. The analysis was carried out on the trajectories after equilibration. The results are summarised in Table 1.

Table 1.  Percentage occupation for detected hydrogen bonds between amino acid residues in the binding pocket and the inhibitors in the four simulated systems, RT-EFV, RT-EMV, RT-ETV and RT-NVP (see definition of atoms for NNRTIs in Figure 1).

NNRTI	Type		% Occupation
EFV	N1-H… O=C	K 101	100
	N1…H-N	K 101	13
	O1…H-N	K 101	73
EMV	N1-H… O=C	K 101	100
	O1…H-NZ	K 103	10
	O1…H-NZ	K 103	4
ETV	N1-H… O=C	K 101	5
	N3-H…OE1	E 698	24
NVP	-		-

As can be seen from the table, the hydrogen bond patterns observed for EFV and EMV are slightly different. For EFV, there are two strong and one weak hydrogen bonds (with occupations of 100%, 73% and 13%) respectively, these were detected with residue K101 only and no hydrogen bond interaction with K103 could be found. For EMV, the ligand binds strongly through a single hydrogen bond to the binding pocket via K101 with a 100% hydrogen bond occupation.

In the case of ETV, a less weak hydrogen bond is observed between the nitrogen of the inhibitor and the carboxylate group of E698 (E138 in the p51 subunit) with a 24% occupation. Such an interaction with E698 is unusual for HIV-1 RT/NNRTIs complexes, and is found only in special cases [37].

For the NVP, no hydrogen bond formation within the binding pocket of RT was found. This observation is in agreement with the X-ray structure of the HIV-1 RT-NVP complex [20], and the suggestions of Kohlstaedt et al. [2] and Das et al. [7], who stated that NVP generally contains aromatic rings and forms π-π interactions with the five aromatic amino acid residues (Y181, Y188, F227, W229 and Y318) in the HIV-1 RT binding pocket without any hydrogen bonds.

Taking into account all the hydrogen bond data given above, EFV and EMV bind much more tightly into the catalytic site of the HIV-1 RT than ETV, whilst for NVP no hydrogen bond interactions were found. These hydrogen bond patterns seem to correspond with the binding affinities and subsequently, with the experimentally derived IC₅₀ values of the drugs. However, the HIV-1 RT binding pocket of NNRTIs is a hydrophobic cavity and, consequently other types of interactions also need to be considered too.

Which RT residues are important for binding?

The interaction energies between inhibitors and individual amino acid residues can be calculated using the decomposition energy module of AMBER 9, in order to identify the residues that are important for the binding affinities. The plot of the decomposition energies (DC) of those amino acids of RT, which are located within 5 Å of the binding pocket, is shown in Figure 2. Their interaction energies vary from −4 to −1 kcal/mol for all the systems. It is interesting to note that the nine amino acid residues: L100, K101, K103, V106, Y181, Y188, F227, W229 and L234, show the largest contributions to the HIV-1 RT/NNRTIs interaction energies. These observations support the previous clinical data [15,38], which found that the mutations in viruses from resistant patients, and thus mutations likely to have been selected for resistance, mainly occur on these residues.

Figure 2.  Per residue interaction energies of the HIV-1 RT to EFV, EMV, ETV and NVP with those residues which contribute most to the inhibitor surroundings.

As can be seen in Figure 2, the interaction energies of RT for six of the above nine amino acid residues namely, L100, V106, V179, V189, L234 and P236, with the four inhibitors are not significantly different. The amino acids of HIV-1 RT with the highest interaction energies for EFV are K101 and K103. This seems to be in contrast to the previous section which focused on hydrogen bonding, where no hydrogen bond with K103 was observed. However, interactions calculated in this section are total interactions where not only hydrogen bonds but also other types of interactions such as electrostatic and van der Waals interactions are included. In the RT/NVP system, the interaction energies with the residues Y181 and Y188 were the highest compared with the other three inhibitors. The simulated interaction energies of K101, K103, Y181 and Y188 are comparable to those found experimentally in which the K103N mutation affects a high resistance to EFV by reducing the rate of inhibitor entry [39], whilst the mutation of the Y181C and Y188L amino acid residues cause high resistance levels against NVP by loss of the favourable aromatic ring interactions [7]. Therefore, the data leads us to conclude that the importance of these residues on the interaction energies is a primary consequence of the mutation of NNRTIs.

Inhibitors/HIV-1 RT binding free energies

The total binding free energies (ΔG_binding) averaged over the trajectories between 1.5 and 3 ns were calculated for the RT/EFV, RT/EMV, RT/ETV and RT/NVP complexes using MM-PBSA. The free energies and decomposed free energies are shown in Table 2.

Table 2.  Calculated binding free energy and its components (kcal/mol) as well as the experimental IC₅₀ (in µM) of the NNRTIs complex with RT.

	RT-EFV	R-EMV	R-ETV	R-NVP
ΔEel	−11.0 ± 2.6	−7.7 ± 2.1	−17.5 ± 5.7	−5.7 ± 1.9
ΔEvdw	−42.4 ± 1.8	−49.7 ± 2.7	−50.3 ± 2.9	−42.1 ± 2.0
ΔEsolute	−53.4 ± 2.5	−57.4 ± 3.3	−67.8 ± 5.0	−47.8 ± 3.2
ΔG^nonpolar	−3.6 ± 0.1	−4.1 ± 0.2	−6.5 ± 0.2	−3.6 ± 0.3
ΔG^el	25.0 ± 2.1	32.7 ± 3.0	43.3 ± 3.4	25.5 ± 2.6
ΔGsol	21.4 ± 2.1	28.6 ± 3.0	36.8 ± 3.3	21.9 ± 2.6
ΔGbinding	−31.8 ± 2.4	−28.8 ± 2.8	−31.0 ± 4.5	−25.9 ± 3.7
	0.001^7		0.002^7	0.085^7
I5	0.001^1		0.0014^1	0.76^1
		0.004^1		0.034^1
	0.004^1		0.029^1	0.039^1

The van der Waals interactions, part of the ΔE_solute, appear to be the major contribution to the HIV-1 RT/inhibitor binding (∼80%) for all complexes in accordance with the fact that the binding pocket of the HIV-1 RT is considerably hydrophobic. Amongst the four complexes, the order of ΔE_vdw interactions in the gas phase (solute) are RT/ETV ∼ RT/EMV > RT/EFV ∼ RT/NVP, whilst those for ΔE_ele interactions are RT/ETV > RT/EFV > RT/EMV > RT/NVP. However, the binding free energies (ΔG_binding) of RT with EFV and ETV are of about the same magnitude, but larger than those of the other two inhibitors, EMV and NVP. The solvation free energy (ΔG_sol) is the most important term for total binding free energy and for EMV and ETV the gas phase interaction energies are larger than those of EFV and NVP, which is likely to be a consequence of the overall free solvent energies of these molecules (including the binding pocket and complexes). This is in agreement with the solvation data where more water molecules were found in the neighbourhood of ETV and EMV than that for EFV. The absolute values of the predicted binding free energies are consistent with the reported experimentally derived IC₅₀ values, as shown in Table 2, but it is important to note that the experimentally derived values depend on the method used [7,12–14].

For the structural differences of the ligands, the entropic contribution of the inhibitors may not be negligible for calculating the binding free energy. Because of the system size we were not able to compute the entropic energy of the entire system because the protein has 1000 residues and is huge. Therefore, one hundred MD snapshots from 1500 to 3000 ps of the full-length RT complexes with a truncated HIV-1 RT protein composing of residues 1–350 were extracted for calculating ΔG_binding including the TΔS. These residues were selected because they constitute the binding pocket for HIV-1 RT and NNRTI inhibitors. We believe that although the inclusion of all residues would be important for determining the binding free energy accurately, the most important contribution should be considered using our 350-residue model. Thus, the relative stabilities of the four NNRTIs inhibitors could be addressed correctly and these results are shown in Table 3. From the calculations, the order of ΔG_binding are RT/EFV > RT/ETV > RT/EMV > RT/NVP. Interestingly, this trend of free energies (where the entropy term is included) is the same as those obtained for a 1000-residue model but without inclusion of the entropic terms. Despite apparent differences in their ligands these inhibitors share the common six-membered heterocyclic ring as the core structure.In their bound state, the conformation of the conserved ring exhibits similar orientations in such a way that the N1 nitrogen of the ring forms a hydrogen bond to K101 (Figure 1) with the exception of NVP. Such a conformational arrangement of the bound NNRTIs is known as a butterfly-like shape and has π-π interactions with the aromatic amino acids (Y181, Y188, F227, W229 and Y318). Therefore, the entropic contribution was not large enough to change the order of the values of ΔG_binding.

Table 3.  Calculated binding free energy (kcal/mol) of the NNRTIs complex with HIV-1 RT for residues position 1-350, residues around the binding pocket.

	RT-EFV	RT-EMV	RT-ETV	RT-NVP
-TΔS	19.7 ± 11.1	20.2 ± 10.3	23.0 ± 10.3	16.9 ± 10.7
ΔGbinding (excluded TΔS)	−57.4 ± 2.5	−54.9 ± 2.9	−59.7 ± 2.9	−47.0 ± 2.5
ΔGbinding (included TΔS)	−37.7 ± 11.4	−34.7 ± 10.7	−36.7 ± 10.7	−30.1 ± 11

Water molecules in the cavity of HIV-1 RT

The distribution of water molecules in the cavity of the enzyme-inhibitor complex is known to play an essential role in drug binding. To analyse this information, the radial distribution functions (RDFs, g_xy(r)), that is the probability of finding the particle of type y in the spherical radius r around the particle of type x, were evaluated. Here, the RDFs for donor and acceptor atoms of the four inhibitors (EFV: N¹, O¹-O², EMV: N¹-N², O¹-O³, ETV: N¹-N⁵, O¹ and NVP: N¹-N⁴, O¹ (see labels in Figure 3), to the oxygen atom of the water molecule were calculated. The results, as well as the running coordination numbers integrated up to the first minimum (marked by an arrow), of the corresponding RDF are summarised in Figure 3, whilst the chemical structures of each inhibitor together with the definition of the atoms are summarised in Figure 4.

Figure 3.  Radial distribution functions, g(r), centred on the inhibitor atoms to the oxygen atoms of modelled water of the four complexes, RT/EFV, RT/EMV, RT/ETV and RT/NVP. The chemical structures of each inhibitor together with the numbering of the atoms and running integration number up to the first minimum (marked as arrow) are also given.

Figure 4.  Snapshot of the MD simulation showing the surrounding water molecules of the first hydration shell of the four NNRTIs in the HIV-1 RT binding pocket.

To estimate the total number of water molecules in the pocket of the HIV-1 RT, any oxygen atom of water lying within the spherical radius of the first minimum of the donor atoms of NNRTIs were counted. The four inhibitors were increasingly solvated in the order of ETV > EMV > NVP ∼ EFV (Figure 3), and this was the case when considered in terms of both the coordination number and distances to solvent molecules. To visualise the above-mentioned hydration, the snapshots which accumulate water molecules lying under the first peak ∼ 4 Å of the RDFs are displayed in Figure 4.

For EFV, the RDF shows a sharp and narrow first peak at ∼2.9 Å, with a corresponding coordination number (CN) of 1, indicating that O¹ is solvated by a single water molecule. In terms of hydrogen bonding (O-O distance 3.0 Å), this accounts for 100% occupation in good agreement with the X-ray structure (PDB: 1FK9). The solvation of the N¹ atom can be accessed by 0.5 water molecules with a distance of ∼4.5 Å (Figure 4), suggesting that this is the same water found at the O¹ site which is swinging back and forth between the nearby binding site. A small broad peak pronounced at 4.5 Å, with a CN of 0.5 water molecules was found for O² indicating a weak hydrogen bond of around 30% occupation between O² of EFV and the water molecule.

For the RDFs of EMV, two water molecules were found in the binding pocket in agreement with the X-ray structure (PDB: 1RT1). The first sharp peak of O¹ and the broad peak of N¹ can be clearly seen, which is pronounced at 3 Å and 4 Å with the CNs of 1 and 0.5 water molecules, respectively. For O², N² and O³, the first peaks were observed at 3, 4.5 and 4.5 Å with CNs of 1, 0.5 and 0.5, respectively. The detailed analysis of the simulated trajectories shows that the various coordination numbers result from one water molecule which is located in the vicinity of the side chain of the inhibitor.

For ETV, solvation of N³ can be accessed by 2.5 water molecules (the configuration is shown in Figure 4) at a distance of ∼3 Å (Figure 3). However, the minimum of the first sharp peak is above zero which demonstrates solvent exchange is occurring between the first and another water molecule with N³. A small broad peak pronounced at 4 Å with a CNs of 0.5, results from the same water molecule found at N³. A small sharp peak, pronounced at 3 Å with a CN of 0.5 water molecules, was found for N⁵ which indicates a strong hydrogen bond between the N⁵ of ETV and the water molecule. However, this hydrogen bond was not seen during the whole simulation time of the MD simulation. No water molecules were detected around the O¹, N², N⁴ and N⁶ atoms.

For NPV, RDFs of N¹ and N³ donor atoms show the first sharp peaks at 3.0 Å, with a CN of 1, which indicates that N¹ and N³ are solvated by one water molecule (this is the same as those found in X-ray; 1VRT). The first broad peak of N² RDF observed at 4.0 Å with a CN of 1 indicates that N² is solvated by the same water molecule as N³. In contrast, N⁴ and O¹ were observed to be free from solvation.

The RDFs calculation is in good agreement with the X-ray structure of EFV, EMV and NVP, whereas the RDF of ETV could not be compared because its X-ray structure does not include water coordinates. Taking into account all the data and the discussion given above, the EFV atoms are much less accessible by water molecules than the other three inhibitors, considered in terms of both coordination number and distance to solvent molecules. This implies that there is less space available between the EFV atoms and the enzyme residues, a conclusion which was strongly supported by the hydrogen bond formation (Table 1) and enzyme/inhibitor binding energies (Table 2), i.e. amongst the four inhibitors, the main contribution of EFV was exhibited in terms of number as well as the percentage occupation of hydrogen bonds, as shown in Table 1.

Conclusion

MD simulations of the HIV-1 RT complexed with EFV, EMV, ETV and NVP provide information on hydrogen bonding, important residues of the HIV-1RT/inhibitor interaction energies, binding free energies and hydration structures of the complexes. The simulation results report good evidence concerned with questions related to the basic mutation and free energy data. With regard to the interaction energies, the nine amino acid residues, L100, K101, K103, V106, Y181, Y188, F227, W229 and L234, show the largest contributions to the enzyme/inhibitor interaction energies. The obtained results support clinical data which revealed that these same nine residues are the most frequent mutated amino acids recovered in viral isolates from NNRTI resistant patients, suggesting selection for these mutations for NNRTI resistance. The data leads us to conclude that the residues which are of high importance for the interaction energy are the primary sources for selected mutations of NNRTIs. In addition, the binding affinities calculated by MM-PBSA, ΔG_binding values were observed in the following order: EFV ∼ ETV > EMV > NVP. This agrees well with the experimental IC₅₀ values. In addition, the distribution of water molecules is comparable to that found in the crystal structures and is in good agreement with enzyme/inhibitor binding energies in which the EFV atoms can be much less accessed by water molecules than by the other three inhibitors.

Acknowledgements

The authors would like to thank the Computational Chemistry Unit Cell, Faculty of Science, Chulalongkorn University, and the Institute of Theoretical Chemistry, University of Vienna, for providing research facilities, software packages and computing time.

Declaration of interest

This work was supported by the Thailand Research Fund and the National Center for Genetic Engineering and Biotechnology (BIOTEC).