help button home button Biophys. J.
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
Originally published as Biophys J. BioFAST on November 3, 2006.
doi:10.1529/biophysj.106.092528
This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow All Versions of this Article:
biophysj.106.092528v1
92/3/798    most recent
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Aruksakunwong, O.
Right arrow Articles by Hannongbua, S.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Aruksakunwong, O.
Right arrow Articles by Hannongbua, S.
Biophysical Journal 92:798-807 (2007)
© 2007 The Biophysical Society

On the Lower Susceptibility of Oseltamivir to Influenza Neuraminidase Subtype N1 than Those in N2 and N9

Ornjira Aruksakunwong *, Maturos Malaisree {dagger}, Panita Decha {dagger}, Pornthep Sompornpisut {dagger}, Vudhichai Parasuk {dagger}, Somsak Pianwanit {dagger} and Supot Hannongbua {dagger}

* Department of Chemistry, Faculty of Science, Rangsit University, Pathumtani, 12000 Thailand; and {dagger} Department of Chemistry, Faculty of Science, Chulalongkorn University, Bangkok 10330, Thailand

Correspondence: Address reprint requests to Supot Hannongbua, Tel.: 66-22-187602; Fax: 66-22-187603; E-mail: supot.h{at}chula.ac.th.


   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 COMPUTATION METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 ACKNOWLEDGEMENTS
 REFERENCES
 
Aiming to understand, at the molecular level, why oseltamivir (OTV) cannot be used for inhibition of human influenza neuraminidase subtype N1 as effectively as for subtypes N2 and N9, molecular dynamics simulations were carried out for the three complexes, OTV-N1, OTV-N2, and OTV-N9. The three-dimensional OTV-N2 and OTV-N9 initial structures were represented by the x-ray structures, whereas that of OTV-N1, whose x-ray structure is not yet solved, was built up using the aligned sequence of H5N1 isolated from humans in Thailand with the x-ray structure of the N2-substrate as the template. In comparison to the OTV-N2 and OTV-N9 complexes, dramatic changes were observed in the OTV conformation in the OTV-N1 complex in which two of its bulky side chains, N-acethyl (–NHAc) and 1-ethylproxy group (–OCHEt2), were rotated to adjust the size to fit into the N1 catalytic site. This change leads directly to the rearrangements of the OTV's environment, which are i), distances to its neighbors, W-178 and E-227, are shorter whereas those to residues R-224, E-276, and E-292 are longer; ii), hydrogen bonds to the two nearest neighbors, R-224 and E-276, are still conserved in distance and number as well as percentage occupation; iii), the calculated ligand/enzyme binding free energies of –7.20, –13.44, and –13.29 kcal/mol agree with their inhibitory activities in terms of the experimental IC50 of 36.1–53.2 nM, 1.9–2.7 nM, and 9.5–17.7 nM for the OTV-N1, OTV-N2, and OTV-N9 complexes, respectively; and iv), hydrogen-bond breaking and creation between the OTV and neighborhood residues are accordingly in agreement with the ligand solvation/desolvation taking place in the catalytic site.


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
COMPUTATION METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 ACKNOWLEDGEMENTS
 REFERENCES
 
The avian influenza A virus (H5N1) is highly pathogenic and has been shown to cross the species barrier and infect humans. Since its first emergence in Hong Kong in 1997 (1Go), the virulent H5N1 infections have caused many human deaths in various countries in Asia and Europe during the past few years. These events have led to global concern about a worldwide pandemic due to the possibility of antigenic shift allowing human-to-human transmission.

The influenza A virus is a negative stranded RNA virus. Its surface contains two membrane glycoproteins, hemagglutinin (HA) and neuraminidase (NA) (2Go). The HA functions as the receptor-binding and membrane fusion glycoprotein in cell entry, whereas the NA functions as the receptor-destroying enzyme by cleaving the {alpha}-ketosidic linkage between terminal sialic acid and an adjacent sugar residue (3Go). Influenza virus HA and NA are both glycosylated proteins, but the terminal sialic acids are released by the enzymatic activity of NA. The virus is further classified into subtypes based on i), the subtype of the glycoproteins HA and NA present on the virus, and ii), ability of the virus to cause disease. Currently, 15 HA subtypes (H1–H15) and nine NA subtypes (N1–N9) have been identified for the influenza A virus (4Go).

Because of the essential role of NA in influenza replication and its highly conserved active site, the NA inhibitors were rationally designed to block this active site of the influenza virus. For the treatment or prophylaxis of influenza infections, the NA inhibitors are a new class of antiinfluenza drugs that have antiviral activities against a broad range of influenza viruses. Several potent and specific inhibitors of the NA enzyme have been developed through structure-based rational drug design, and two inhibitors, zanamivir (Relenza) and oseltamivir (OTV) (Tamiflu), have been approved for human use (5Go,6Go).

Many crystal structures of an influenza virus complexed with their substrate and inhibitors have been determined (7Go–10Go). Structure analysis of the NA enzyme demonstrated that there are several invariant residues which form part of an active site center. Some of these residues interact directly with the substrate (functional residues) and others scaffold for the functional residues (framework residues) (11Go). These residues include three arginines at positions 118, 292, and 371. The x-ray structures of N2 and N9 complexed with substrate and inhibitor show that these three residues form hydrogen bonds in the same manner with the carboxylate group of the sugar substrate and side chain of the inhibitor (2Go,7Go,9Go). Other importance residues are E-276 and E-119, which were observed to form hydrogen bonds with the substrate (2Go). It was also reported that interactions between zanamivir or OTV and the residues in the NA active center are similar to those with the natural substrate (2Go). Although the active sites of N2 and N9 are highly conserved, the inhibitory activity of OTV leads to the mean IC50 of 1.9–2.7 nM and 9.5–17.7 nM, respectively (12Go).

A key question is why NA subtype N1 has a mean IC50 of 36.1–53.2 nM, significantly higher than those of N2 and N9 (12Go). Unfortunately, the x-ray structure of the OTV-H5N1 complex is not yet solved. The only information helping to understand this fact is the proposal of Moscona (13Go) based on the x-ray structure of OTV and zanamivir complexed with NA N9. The active site of N1 was proposed to change its shape to accommodate the bulky side chain of OTV, whereas such change is not required for zanamivir. This rearrangement involves the rotation of E-276 and its bonding to R-224 (13Go).

To provide an answer to the above question, it is our intention to investigate structural, dynamic, and thermodynamic properties as well as solvation of the three NA subtypes, N1, N2, and N9, complexed with OTV using molecular dynamics (MD) simulations. Here, the initial structure of the OTV-N1 complex was built up using the aligned sequence H5N1 isolated from humans in Thailand and the x-ray structure of the N2-substrate as the template.


   COMPUTATION METHODS
TOP
 ABSTRACT
 INTRODUCTION
 COMPUTATION METHODS
RESULTS AND DISCUSSION
 CONCLUSIONS
 ACKNOWLEDGEMENTS
 REFERENCES
 
Validation of the N1 model
The sequence of the human influenza NA N1, which is the highly pathogenic H5N1 A virus isolated during the 2004 influenza outbreaks in Thailand, was taken from Genbank LOCUS ID AAS89006 (A/Thailand/3(SP-83)/2004). This sequence consists of 449 amino acids. Currently, the crystal structure of the N1 protein from H5N1 is not yet solved. However, the three-dimensional (3D) structure of N1 can be modeled based on the available crystal structures. In this study, the crystal structure of N2 complexed with substrate (Protein Data Bank (PDB) code 2BAT, 2.0-Å resolutions) and N9 complexed with OTV (PDB code 2QWK, 1.8-Å resolutions) were selected as the template. Note that the sequences of N2 and N9 contain 388 and 387 amino acids, respectively. The N1, N2, and N9 sequence alignment was performed using the homology module in the InsightII program (Accelrys, San Diego, CA). The final alignment was shown in Fig. 1. The N1/N2, N1/N9, and N2/N9 identities were 46%, 49%, and 48%, respectively.


Figure 1
View larger version (53K):
[in this window]
[in a new window]

  		FIGURE 1  Sequence alignment of N1 (A/Thailand/3(SP-83)/2004) with N2 (PDB code 2BAT) and N9 (PDB code 2QWK). The N1 protein sequence used in this study is from the highly pathogenic human H5N1 A virus isolated during the 2004 influenza outbreak in Thailand.

 
Three 3D models of N1 were built up using the sequence alignments and the atomic coordinates of the N2-substrate, N9-OTV, and both structures. They are, for simplicity, named N1–(N2), N1–(N9), and N1–(N2, N9), respectively. The Needleman-Wunsch algorithm (14Go) was used for pairwise alignment to identify the correspondence region. The scores were obtained for each pair of amino acid in the sequences, i.e., gaps can be placed in the proper locations and high scores can be yielded for good matches. The coordinates of the structurally conserved regions (SCRs)—in which amino acid sequences in the template and the model are identical—of the template were copied to be those of the model. For the similar SCRs, only the backbone coordinates were replicated and the side-chain atoms were added in a manner that preserves the model protein's residue type. The coordinates for the loop or variable regions were generated by finding peptide segments in other proteins that fit properly into the model's spatial environment.

Superimpositions between the three models, N1–(N2)//N1–(N9), N1–(N2)//N1–(N2,N9), and N1–(N9)//N1–(N2,N9), were performed. The corresponding root mean-square deviations (RMSDs) are 1.46, 1.39, and 1.41 Å, respectively. Inconsistent with the percentage identity between the three sequences mentioned above, the RMSD values indicate that the three 3D models do not differ significantly.

From the characterization of the 1918 Spanish influenza virus (H1N1) NA gene, it was shown that this N1 protein sequence was placed within and near the root of the avian sequence, N2. Those sequences share many characteristics with both mammalian and avian viruses. In that study, the authors suggest that it is possible to align the N1 sequence with N2 subtype and examine the N2 antigenic site for evidence of drift in N1 (15Go). With these reasons, the model built from N2, N1–(N2), was chosen to represent NA subtype N1 and used for MD simulations. The final 3D structure of the selected model superimposed with the N2 template was displayed in Fig. 2.


Figure 2
View larger version (41K):
[in this window]
[in a new window]

  		FIGURE 2  3D structure of N1 modeled according to the template of 2BAT and the crystal structure of N2 complexed with the substrate (sialic acid) (a). The NA enzyme is shown in ribbon, where the models N1 and N2 are colored in green and color by element, respectively. The substrate is presented in yellow. The two-dimensional structures of the substrate (sialic acid) (b) and OTV (c) are also displayed.

 
Protein preparation for MD simulations
The x-ray structure of N9-OTV (OTV-N9) used for MD simulation was taken from the PDB. To prepare the OTV-N1 and OTV-N2 complexes, superimpositions of the following pairs were performed, N1 (model) with OTV-N9 and N2 (x-ray structure excluded substrate) with OTV-N9. The N9 enzyme coordinates were, then, removed and the OTV-N1 and OTV-N2 were, then, obtained. All missing hydrogens of the protein were added using the LEaP module in the AMBER 7 software package (16Go).

The starting structure and force field parameters for the inhibitor were obtained as follows. Hydrogens were added to the x-ray coordinates of OTV (2QWK) 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 lengths involving hydrogens. Then, the restrained electrostatic potential fitting procedure was employed to calculate partial atomic charges of the inhibitor (17Go). Force field parameters of the inhibitor were assigned based on the atom types of the force field model (18Go). Gaussian 98 (19Go) 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 (20Go).

To incorporate the solvent and counterions under consideration, each system was solvated with 10 Å of the TIP3P water model (21Go) and neutralized by the counterions using the LEaP module. The crystallographic water molecules were also kept in the 78 x 86 x 80 simulation box. This leads to total atoms of 43998, 43905, and 41572 for the OTV-N1, OTV-N2, and OTV-N9, respectively.

Molecular dynamics simulations
Three MD simulations were performed for the OTV-N1, OTV-N2, and OTV-N9 complexes, where their initial structures were described above. Energy minimization and MD simulations were carried out using the SANDER module of AMBER 7 (16Go) with the Cornell force field (18Go). The periodic boundary condition with the NPT ensemble was employed. A Berendsen coupling time of 0.2 ps was used to maintain the temperature and standard pressure of the system (22Go). The SHAKE algorithm (23Go) was applied to constrain all bonds involving hydrogens. The simulation time step of 2 fs was used. All MD simulations were run with a 12-Å residue-based cutoff for nonbonded interactions and the particle mesh Ewald method was used for an adequate treatment of long-range electrostatic interactions (24Go). The simulation steps consist of thermalization, equilibration, and production phases. Initially, the temperature of the system was gradually heated from 0 K to 298 K during the first 60 ps. Then, the system was maintained at 298 K until MD reached 1 ns of the simulation. Finally, the production phase started from 1 ns to 2 ns of the simulation. The convergence of energies, temperature, pressure, and global RMSD was used to verify the stability of the system.

The MD trajectory was collected every 0.2 ps. 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 (25Go). Graphic visualization and presentation of protein structures were done using RasMol, Swiss-Pdb Viewer (26Go), WebLab Viewer (Accelrys), and MolScript (27Go).

The binding free energy ({Delta}Gbinding) calculations
In general, the free energy of the inhibitor binding, {Delta}Gbinding, can be obtained from the difference between the free energy of the receptor-ligand complex (Gcpx) and the unbound receptor (Grec) and ligand (Glig) as follows:

Formula 1(1)

The MM/PBSA approach calculates {Delta}Gbinding on the basis of a thermodynamic cycle. Therefore, Eq. 1. can be approximated as

Formula 2(2)
where {Delta}EMM is related to the enthalpic changes in the gas phase upon binding and is obtained from molecular mechanics van der Waals and electrostatic energies, T{Delta}S involves the entropy effect, and {Delta}Gsol is the free energy of solvation. Then {Delta}Gsol is composed of the electrostatic and nonpolar contributions (28Go), and therefore can be expressed as

Formula 3(3)
where {Delta}GPB was calculated using a continuum solvent model with a Poisson-Boltzmann solution (29Go), and {Delta}GSA is estimated from the solvent-accessible surface area (SASA) (30Go). The {Delta}Gbinding was obtained using the MM/PBSA module in the program AMBER 7, which interfaces the program DelPhi 4 (31Go). To calculate electrostatic free energy of solvation, the grid resolution of 0.33 Å with the boundary conditions of Debye-Huckel potentials was employed. Atomic charges were taken from the Cornell force field (18Go). The water and protein dielectric values were set to 80 and 1, respectively. The SASA was calculated using a 1.4-Å probe radius. Atomic radii were taken from the PARSE parameter set (30Go). The nonpolar free energy of solvation was obtained by 0.00542 x SASA + 0.92 kcal/mol (30Go). In the study, the contribution of the entropy (T{Delta}S) was not included. An estimation of the entropy effect from normal mode analysis requires high computational demands. The entropic difference, among the three systems, should be very small because all system models, NA subtypes N1, N2, and N9 complexed with OTV, are almost similar. In addition, interpretation will focus only on the relative values of the binding free energy.


   RESULTS AND DISCUSSION
TOP
 ABSTRACT
 INTRODUCTION
 COMPUTATION METHODS
 RESULTS AND DISCUSSION
CONCLUSIONS
 ACKNOWLEDGEMENTS
 REFERENCES
 
Overall enzyme and inhibitor structures
RMSDs with respect to the initial configuration for the overall structure of NA and OTV were separately calculated and plotted in Fig. 3. Substantial changes were observed for the OTV-N1 complex, where the RMSD of the enzyme increases slowly during the first ~700 ps, whereas that of the inhibitor raises rapidly at ~1 ns. No significant difference was found for the RMSDs of the other two systems, OTV-N2 and OTV-N9.


Figure 3
View larger version (22K):
[in this window]
[in a new window]

  		FIGURE 3  RMSD plots of protein backbone atoms (a) and all atoms of OTV (b) with respect to simulation time.

 
Changes of the oseltamivir structure
Detailed evaluations were focused on the changes of the molecular structure of OTV in the OTV-N1 complex (Fig. 3 b). Structure changes were monitored in terms of torsional angles, tor1–tor4, defined, respectively, by the following sets of four atoms, C9-C2-C1-O1, C4-C5-N1-C6, C5-C8-O4-C10, and C8-O4-C10-C12 in Fig. 4 a. The plots were for all three systems, evaluated after equilibration, at 1.0 ns. The results are compared in Fig. 4, b–e. Note that all torsional angle plots show sharp peaks, indicating that the side chains were rotated only in a narrow range around their preferential configurations.


Figure 4
View larger version (20K):
[in this window]
[in a new window]

  		FIGURE 4  Definition of OTV's torsional angles, tor1–tor4 (a) and their distributions, (b–e), in the OTV-N1, OTV-N2, and OTV-N9 complexes, sampling after equilibration, at 1.0 ns.

 
Considering the –COO– side chain of OTV (see Fig. 4 a), tor1 for the three systems, the OTV in N1, N2, and N9 appears at almost the same angle, 2° (Fig. 4 b). This means that the –COO– group was observed to lie almost in the ring plane.

Dramatic changes were detected in tor2–tor4 of OTV in the NA N1 compared with the initial conformation of OTV, whereas such changes were not found for the N2 and N9 complexes. The tor2 angle in OTV-N2 and OTV-N9 appears between –100° and –120° (Fig. 4 c), indicating that the C5 -N1 bonds (see the –NHAc group in Fig. 4 a) of the OTV in the two complexes were located above the ligand's plane and slightly tilted, by 10–30°, to the –OCHEt2 group. The tor3 angle of 92° for both systems (Fig. 4 d) implies that the C8-O4 bond (see Fig. 4 a) lies perpendicular to the ligand plane. In addition, the detected values of tor4 for both systems at –162° (Fig. 4 e) leads to the conclusion that the end of this side chain points backward to the OTV molecular plane.

In comparing the OTV-N1 with the OTV-N2 and OTV-N9 systems, a different situation was observed for the conformation of OTV in the OTV-N1 complex in which the tor2–tor4 angles were turned by 40–60° (from –100° or –120° to –160°), 55° (from 92° to 147°), and 56° (from –162° to –98°), respectively. Changes of the tor2 and tor3 angles indicate the rotation of the C5-N1 bond in the direction approaching the –OCHEt2 group and rotation of the C8-O4 bond in the direction toward the –COO– side chain (see Fig. 4 a), respectively. Interestingly, the tor4 angle of OTV-N1 was turned backward, from –162° to –98° (Fig. 4 e). This means that the end of the –OCHEt2 side chain, defined by the O4-C10 bond, for the three systems points with the same angle to the molecular plane.

In conclusion, changes of the OTV structure in the NA N1 pocket can directly affect its interaction with the enzyme as well as the structure of its environment, solvation, and nearest enzyme residues (see details in the next sections), whereas this is not the case for OTV in N2 and N9.

Changes of the neuraminidase cavity
As mentioned earlier, changes of the OTV structure (rotation of tor2 – tor4 as shown in Fig. 4) can directly affect the interaction, position, and orientation of the neighboring residues in the NA pocket. To visualize and quantify such effects, distances between the atoms of OTV and its neighboring residues, R-118, R-152, W-178, R-224, E-227, E-276, and R-292, were measured for the three systems and plotted in Fig. 5.


Figure 5
View larger version (30K):
[in this window]
[in a new window]

  		FIGURE 5  Distance to NA residues located in the neighborhood of the OTV in the OTV-N1, OTV-N2, and OTV-N9 systems.

 
Effects of the rotation of the OTV side chain can be clearly seen from this plot. Compared with OTV-N2, which was used as the template in building up the initial OTV-N1 structure by homology modeling, changes of the OTV conformation do not affect the positions of the R-118 and R-152 residues (Fig. 5, b and c). Rearrangements of the N1's active site around the –NHAc and the –OCHEt2 groups were significantly observed to accommodate the OTV's bulky side chains. Rotation of the –NHAc side chain (tor2 in Fig. 4 a) in the direction to the –OCHEt2 group shifts the W-178 (Fig. 5 d) and E-227 (Fig. 5 f) closer to the OTV by more than 2 Å. Consequently, the R-224, E-276, and R-292 residues were repelled by the bulky –OCHEt2 group to move away by more than 1 Å (Fig. 5, e, g, and h). The above findings may be one of the important reasons OTV is less susceptible to human influenza NA subtype N1 than those in N2 and N9. The importance of R-292 in OTV-N2 and OTV-N9 is signified by the existence of the mutant strain R-292K (32Go).

An interaction between the R-224 and E-276 residues in the neighborhood of the OTV molecule was found indirectly from the x-ray structure of the NA N9 complexed with OTV (33Go). With the distances shown in Table 1 and the conformation shown in Fig. 6 a, this suggested clearly that two hydrogen bonds were formed. Therefore, this type of binding was supposed to play a role either in stabilizing the enzyme structure or in the catalytic process.


View this table:
[in this window]
[in a new window]

  		TABLE 1  Percentage occupations (%) and the corresponding distances (R, in Å) of direct hydrogen bonding

 

Figure 6
View larger version (37K):
[in this window]
[in a new window]

  		FIGURE 6  Hydrogen bonds between the two residues R-224 and E-276 of the NAs extracted from (a) the x-ray structure of N9 (PDB code 2QWK) and the MD simulations for the (b) OTV-N1, (c) OTV-N2, and (d) OTV-N9 systems.

 
To monitor this event for the NA N9 in solution as well as those of the N1 and N2 subtypes, hydrogen bonding between the two residues was analyzed using the CARNAL module in AMBER 7 with the following criteria: i), the distance between proton donor (D) and acceptor (A) atoms was ≤3.5 Å, and ii), the D-H....A angle was ≥120°. The results are summarized in Table 1. The corresponding configurations for the three systems, as well as that of the x-ray data, are depicted in Fig. 6.

Effects of the rotation of the side chains of OTV (tor2–tor4 in Fig. 4) in the OTV-N1 complex can be seen in terms of the R-224–E-276 binding shown in Table 1 and Fig. 6. The position and orientation of the R-224 and E-276 residues in the crystallographic and liquid phases are noticeably different. In the solution, the numbers of detected hydrogen bonds are three for OTV-N9 and one for OTV-N1 and OTV-N2, whereas the x-ray structure has two hydrogen bonds. This observation indicates destabilization in the vicinity of the NA N1 and N2 pockets in comparison with that of N9, both in solid and liquid phases.

Interestingly, the positions of R-224 and E-276 in the OTV-N1 complex were observed to shift by more than 1 Å away from the OTV molecule (Fig. 5, e and g) in comparison with OTV-N2. However, the relative orientation between the two residues remains unchanged, i.e., distances as well as numbers and percentage of hydrogen-bond occupation of the two systems are very similar (Table 1). This confirms the previous reports that proposed that the active site of the NA is highly conserved (11Go).

Taking into account the structural data summarized above, the OTV structure does not fit well into the NA N1 cavity, leading to the rotation of its side chain in the OTV-N1 complex. This directly affects its binding with the neighboring residues in the pocket of the NA. Consequently, the OTV's side-chain rotation was found to significantly influence the position and orientation of the five nearest neighbors, W-178, R-224, E-227, E-276, and R-292. The movements of E-276 and R-224 residues appear, in agreement to some degree with the proposal by Moscona (13Go), to be the source of OTV resistance for the NA subtype N1.

The molecular mechanism of NA resistance to OTV was proposed by Moscona (13Go). The NA active site must go through a rearrangement to change its shape to accommodate the bulky side chain of OTV, whereas such a change is not required for zanamivir. This rearrangement involves the rotation of the E-276 residue and its bonding to R-224. The mutations R-292K, N-294S, and H-274Y inhibit this rotation, resulting in resistance to OTV. As Moscona did not state how the molecular structures of both complexes were obtained, we assume this proposed mechanism to be derived from information from the x-ray structures of OTV and zanamivir complexed with NA N9. In our work, such rearrangements in the NA N1, N2, and N9 pockets as well as dynamical properties and solvent effects were additionally quantified.

Solvation structure and oseltamivir-enzyme hydrogen bonding
To monitor the ligand's solvation, especially its changes due to the rotation of the side chains of OTV (tor2–tor4 in Fig. 4) and radial distribution function (RDF, gxy(r))—the probability of finding a particle of type x in a sphere of radius r around the particle of type y—from acceptor atoms of OTV (O1, O2, O3, O4, N1, and N2 as labeled in Fig. 4 a) to an oxygen atom of water, were evaluated for the three systems and plotted in Fig. 7. In addition, ligand-enzyme hydrogen bonding was analyzed using the same criteria as applied for those of the R-224 and E-276 residues (Table 1 and Fig. 6). The results are summarized in Table 2.


Figure 7
View larger version (16K):
[in this window]
[in a new window]

  		FIGURE 7  RDFs centered at donor atoms of OTV (labeled in Fig. 4 a) to oxygen atoms of water molecules for the three simulated systems, OTV-N1, OTV-N2, and OTV-N9, evaluated after equilibration, at 1.0 ns.

 

View this table:
[in this window]
[in a new window]

  		TABLE 2  Percentage occupation of direct hydrogen bonding between OTV (atom types were labeled as in Fig. 4 a) and specific residues of NAs N1, N2, and N9 during MD simulations

 
In terms of the solvent accessible area, space around the inhibitor can be clearly classified into two zones. No water was detected to approach the O3, N1, and N2 atoms of OTV due to the presence of R-152, E-119, D-151, W-178, and E-227. Those residues were found to form hydrogen bonds with the O3, N1, and N2 atoms (Table 2), preventing them from being accessed by water molecules. The radius of the water-free region is 3–4 Å for OTV-N2 and OTV-N9 and 5–6 Å for the OTV-N1 system (Fig. 7, c, e, and f).

The situation is different for the regions around the O1, O2, and O4 atoms of OTV where almost all RDFs show a sharp first peak centered at 2.8 Å (Fig. 7, a, b, and d) with the corresponding coordination number, integrated up to the first minimum. A clear conclusion is that those atoms were solvated by forming hydrogen bonding with one water molecule. Characteristics of the plots are fully consistent with the ligand-enzyme hydrogen-bond formation shown in Table 2. For example the O1 atom was found to form four hydrogen bonds with R-292 and R-371 (Nos. 7, 8, 9, and 11 in Table 2) for the OTV-N2 and OTV-N9 but only one hydrogen bond (No. 11 in Table 2) for OTV-N1 system. That means that R-292 in OTV-N1 was shifted apart from the O1 atom, generating empty space that can accumulate a water molecule to form a strong hydrogen bond with O1. This leads to a much higher and sharper first peak of the RDF centered at the O1 atom of OTV than those of O2 and O4.

For the O2 atom of OTV, the numbers of detected hydrogen bonds are two (with R-118 and R-371, Nos. 1 and 12 in Table 2), four (two with R-118 and two with R-371, Nos. 1, 2, 10, and 12 in Table 2), and two (with R-371, Nos. 10 and 12 in Table 2) for OTV-N1, OTV-N2, and OTV-N9, respectively. Only in the OTV-N2 system was the hydrogen bond to R-118 not detected. These are the reasons the first peak of the N2 RDF is lower than those of N1 and N9. Correspondence between the two sources of data, ligand-water RDFs, and ligand-enzyme hydrogen bonds, was clearly observed for the O4 atom, where only one hydrogen bond with R-152 (Table 2, No. 4) was found only for the OTV-N1 system, i.e., water molecules were prevented from approaching the O4 atom. This leads directly to the flattening of its RDF up to 4 Å (Fig. 7 d). This is not the case for the other two systems, OTV-N2 and OTV-N9, where the O4 atom was evidently solvated by a water molecule.

Enzyme/inhibitor interactions
Interactions between OTV and enzymes were analyzed in terms of binding free energies computed using MM/PBSA calculations. The results in Table 3 support the above conclusions regarding lower susceptibility of the inhibitor to the NA subtype N1 in comparison with N2 and N9. The calculated total binding free energy for OTV-N1 is –7.20 kcal/mol, whereas it is –13.44 kcal/mol and –13.29 kcal/mol for OTV-N2 and OTV-N9, respectively. These values agree well with the experimental free energies of –10.2 to –10.0, –12.0 to –11.7, and –11.0 to –10.6 kcal/mol (converted from the IC50 values of 36.1–53.2 nM, 1.9–2.7 nM, and 9.5–17.7 nM) for OTV-N1, OTV-N2, and OTV-N9, respectively.


View this table:
[in this window]
[in a new window]

  		TABLE 3  Comparison of components of calculated binding free energies* of the three simulated systems: OTV-N1, OTV-N2, and OTV-N9

 
Referred to the calculation details where the entropic term, T{Delta}S, was neglected based on the assumption that the difference would be small among the three simulated systems due to their similarity. This was, somehow, supported by the thermodynamics study of HIV-1 protease complexed with saquinavir, ritonavir, nelfinavir, and indinavir, which show that the binding entropy for the four drugs are –14.0, –11.2, –14.2, and –15.7 kcal/mol, respectively (34Go). The differences are relatively small although sizes of the four drugs are relatively difference. Therefore, the order of the binding free energy shown in Table 3 for the simulated systems should not be changed due to the exclusion of the entropic term. In addition, interpretation will focus only on the relative values of the binding free energy.

These binding interactions correspond to the biological activity of the drug toward the various strains of viruses as can be seen from Table 3. Very strong bindings are observed for all complexes in the gas phase. However, these interactions are offset by the complex-solvent interaction. Negative electrostatic and van der Waals interactions for these complexes imply that the binding is mainly of hydrogen-bonding type, in agreement with our previous conclusion.


   CONCLUSIONS
TOP
 ABSTRACT
 INTRODUCTION
 COMPUTATION METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
ACKNOWLEDGEMENTS
 REFERENCES
 
Rotation of the –NHAc and –OCHEt2 side chains of the OTV, which leads consequently to the rearrangements of the N1's cavity, is a primary source of the lower susceptibility of the OTV to influenza NA subtype N1 than to N2 and N9. The observed data suggest that the size of these two side chains is reduced in a way that maintains their interactions with the neighborhood residues of the N1 protein.


   ACKNOWLEDGEMENTS
TOP
 ABSTRACT
 INTRODUCTION
 COMPUTATION METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 ACKNOWLEDGEMENTS
REFERENCES
 
The authors thank the Computer Center for Advanced Research, Faculty of Science, Chulalongkorn University, and the Computational Chemistry Unit Cell, Department of Chemistry, Faculty of Science, Chulalongkorn University for use of their computing facilities throughout this work. Moreover, we thank Prof. Dr. Yong Poovorawan at the Faculty of Medicine, Chulalongkorn University for many suggestions and comments.

This work was supported by Thailand Research Fund Senior Scholar Grant No. RTA468008.

Submitted on July 5, 2006; accepted for publication October 6, 2006.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 COMPUTATION METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 ACKNOWLEDGEMENTS
 REFERENCES
 
1. Bender, C., H. Hall, J. Huang, A. Klimov, N. Cox, A. Hay, V. Gregory, K. Cameron, W. Lim, and K. Subbarao. 1999. Characterization of the surface proteins of influenza A (H5N1) viruses isolated from humans in 1997–1998. Virology. 254:115–123.[CrossRef][Medline]

2. McKimm-Breschkin, J. L. 2000. Resistance of influenza viruses to neuraminidase inhibitors—a review. Antiviral Res. 47:1–17.[CrossRef][Medline]

3. Ha, Y., D. J. Stevens, J. J. Skehel, and D. C. Wiley. 2002. H5 avian and H9 swine influenza virus haemagglutinin structures: possible origin of influenza subtypes. EMBO J. 21:865–875.[CrossRef][Medline]

4. Puthavathana, P., P. Auewarakul, P. C. Charoenying, K. Sangsiriwut, P. Pooruk, K. Boonnak, R. Khanyok, P. Thawachsupa, R. Kijphati, and P. Sawanpanyalert. 2005. Molecular characterization of the complete genome of human influenza H5N1 virus isolates from Thailand. J. Gen. Virol. 86:423–433.[Abstract/Free Full Text]

5. Gubareva, L. V., L. Kaiser, and F. G. Hayden. 2000. Influenza virus neuraminidase inhibitors. Lancet. 355:827–835.[CrossRef][Medline]

6. Wang, M. Z., C. Y. Tai, and D. B. Mendel. 2002. Mechanism by which mutations at His274 alter sensitivity of influenza A virus N1 neuraminidase to oseltamivir carboxylate and zanamivir. Antimicrob. Agents Chemother. 46:3809–3816.[Abstract/Free Full Text]

7. Smith, B. J., P. M. Colman, M. V. Itzstein, B. Danylec, and J. N. Varghese. 2001. Analysis of inhibitor binding in influenza virus neuraminidase. Protein Sci. 10:689–696.[Abstract/Free Full Text]

8. McKimm-Breschkin, J. L., A. Sahasrabudhe, T. J. Blick, M. McDonald, P. M. Colman, G. J. Hart, R. C. Bethell, and J. N. Varghese. 1998. Mutations in a conserved residue in the influenza virus neuraminidase active site decreases sensitivity to Neu5Ac2en-derived inhibitors. J. Virol. 72:2456–2462.[Abstract/Free Full Text]

9. Varghese, J. N., J. L. McKimm-Breschkin, J. B. Caldwell, A. A. Kortt, and P. M. Colman. 1992. The structure of the complex between influenza virus neuraminidase and sialic acid, the viral receptor. Proteins. 14:327–332.[CrossRef][Medline]

10. Varghese, J. N., V. C. Epa, and P. M. Colman. 1995. Three-dimensional structure of the complex of 4-guanidino-Neu5Ac2en and influenza virus neuraminidase. Protein Sci. 4:1081–1087.[Abstract]

11. Colman, P. M., P. A. Hoyne, and M. C. Lawrence. 1993. Sequence and structure alignment of paramyxovirus hemagglutinin-neuraminidase with influenza virus neuraminidase. J. Virol. 67:2972–2980.[Abstract/Free Full Text]

12. Govorkova, E. A., I. A. Leneva, O. G. Goloubeva, K. Bush, and R. G. Webster. 2001. Comparison of efficacies of RWJ-270201, zanamivir, and oseltamivir against H5N1, H9N2, and other avian influenza viruses. Antimicrob. Agents Chemother. 45:2723–2732.[Abstract/Free Full Text]

13. Moscona, A. 2005. Neuraminidase inhibitors for influenza. N. Engl. J. Med. 353:1363–1373.[Free Full Text]

14. Durbin, R., S. Eddy, A. Krogh, and G. Mitchison. 2004. Needleman-Wunsch Algorithm. 1–3. Cambridge University.

15. Reid, A. H., T. G. Fanning, T. A. Janczewski, and J. K. Taubenberger. 2000. Characterization of the 1918 "Spanish" influenza virus neuraminidase gene. Proc. Natl. Acad. Sci. USA. 97:6785–6790.[Abstract/Free Full Text]

16. Case, D. A., J. W. Pearlman, T. E. Caldwell, J. Cheatham III, W. S. Wang, C. L. Ross, T. A. Simmerling, K. M. Darden, R. V. Merz, A. L. Stanton, J. J. Cheng, M. Vincent, V. Crowley, H. Tsui, R. J. Gohlke, Y. Radmer, J. Duan, I. Pitera, G. L. Massova, U. C. Seibel, P. K. Singh, and P. A. Kollman. 2002. AMBER 7. University of California, San Francisco, CA.

17. Cornell, W. D., P. Cieplak, C. I. Bayly, and P. A. Kollmann. 1993. Application of RESP charges to calculate conformational energies, hydrogen bond energies, and free energies of solvation. J. Am. Chem. Soc. 115:9620–9631.[CrossRef]

18. Cornell, W. D., P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell, and P. A. Kollman. 1995. A second generation forcefield for the simulation of proteins, nucleic-acids, and organic-molecules. J. Am. Chem. Soc. 117:5179–5197.[CrossRef]

19. Frisch, M. J., G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, V. G. Zakrzewski, J. A. Montgomery, J. Stratmann, J. C. Burant, S. Dapprich, J. M. Millam, and others. 2002. Gaussian 98. Gaussian, Inc., Pittsburgh, PA.

20. Wang, J. M., W. Wang, and P. A. Kollman. 2001. Antechamber: an accessory software package for molecular mechanical calculations. J. Am. Chem. Soc. 222:U403. (Abstr.)

21. Jorgensen, W. L., J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein. 1983. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79:926–935.[CrossRef]

22. Berendsen, H. J. C., J. P. M. Postma, W. F. V. Gunsteren, A. DiNola, and J. R. Haak. 1984. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81:3684–3690.[CrossRef]

23. Ryckaert, J. P., G. Ciccotti, and H. J. C. Berendsen. 1997. Numerical integration of the Cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 23:327–341.[CrossRef]

24. York, D. M., T. A. Darden, and L. G. Pedersen. 1993. The effect of long-range electrostatic interactions in simulations of macromolecular crystals: a comparison of the Ewald and truncated list methods. J. Chem. Phys. 99:8345–8348.[CrossRef]

25. Laskowski, R. A., J. A. Rullmannn, M. W. MacArthur, R. Kaptein, and J. M. Thornton. 1996. AQUA and PROCHECK-NMR: programs for checking the quality of protein structures solved by NMR. J. Biomol. NMR. 8:477–486.[Medline]

26. Guex, N., and M. C. Peitsch. 1997. SWISS-MODEL and the Swiss-PdbViewer: an environment for comparative protein modeling. Electrophoresis. 18:2714–2723.[CrossRef][Medline]

27. Kraulis, P. J. 1991. MOLSCRIPT: a program to produce both detailed and schematic plots of protein structures. J. Appl. Crystallogr. 24:946–950.[CrossRef]

28. Srinivasan, J., T. E. Cheatham, P. Cieplak, P. A. Kollman, and D. A. Case. 1998. Continuum solvent studies of the stability of DNA, RNA, and phosphoramidate-DNA helices. J. Am. Chem. Soc. 120:9401–9409.[CrossRef]

29. Gilson, M. K., K. Sharp, and B. Honig. 1987. Calculating electrostatic interactions in biomolecules: method and error assessment. J. Comput. Chem. 9:327–335.[CrossRef]

30. Sitkoff, D., K. A. Sharp, and B. Honig. 1994. Accurate calculation of hydration free energies using macroscopic solvent models. J. Phys. Chem. 98:1978–1988.[CrossRef]

31. Rocchia, W., E. Alexov, and B. Honig. 2001. Extending the applicability of the nonlinear Poisson-Boltzmann equation: multiple dielectric constants and multivalent ions. J. Phys. Chem. B. 105:6507–6514.[CrossRef]

32. Roche Pharmaceuticals. 2005. Tamiflu (oseltamivir phosphate) capsules and for oral suspension. Nutley, NJ.

33. Varghese, J. N., P. W. Smith, S. L. Sollis, T. J. Blick, A. Sahasrabudhe, J. L. McKimm-Breschkin, and P. M. Colman. 1998. Drug design against a shifting target: a structural basis for resistance to inhibitors in a variant of influenza virus neuraminidase. Structure. 6:735–746.[Medline]

34. Todd, M. J., L. Luque, A. Velazquez-Campoy, and E. Freire. 2000. Thermodynamic basis of resistance to HIV-1 protease inhibition: calorimetric analysis of the V82F/I84V active site resistant mutant. Biochemistry. 39:11876–11883.[CrossRef][Medline]





This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow All Versions of this Article:
biophysj.106.092528v1
92/3/798    most recent
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Aruksakunwong, O.
Right arrow Articles by Hannongbua, S.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Aruksakunwong, O.
Right arrow Articles by Hannongbua, S.

HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
Copyright © 2007 by the Biophysical Society.