**4.2. Computational details**

Six systems were modelled and simulated for oseltamivir bound H5N1, and H1N1pdm wild type and H274Y and N294S mutants. In simSMD1, steered molecular dynamics (SMD) simulations were used to remove oseltamivir from its stable binding site in H5N1 neuraminidase. In simFEQ1-10, equilibration simulations used a starting point generated from simSMD1 in which Oseltamivir has undergone an axial rotation such that it is partially displaced from its binding site. In total, 680 ns of simulations were carried out on system sizes of approximately 35,000 atoms.

The "Ensemble" column lists the variables held constant during the simulations; N, p, T, and V correspond to the number of atoms, pressure, temperature, and volume respectively. Under "Type", EQ denotes equilibration, and SCV denotes constant velocity SMD simulation with the speed of 0.5 Å/ns. SimSMD1 is a steered MD simulation with the starting structure from equilibrated simEQ1 ("preflip" drug position). In simFEQ1 to simFEQ10, the starting structure reflected a "flipped" position of oseltamivir from simSMD1 after 7.5ns of simulation, whereas several of its main stabilizing hydrogen bonds to the protein had already been ruptured. Crystallographically resolved water molecules and a structurally relevant calcium ion near the native binding site for SA were retained and modelled in all simulated systems. The protein complexes were then solvated in a TIP3P water box 23 and ionized by NaCl (0.152M) to mimic physiological conditions.

All simulations were performed using NAMD 2.7 and the CHARMM31 force field (Mackerrel et al., 1998). The ionized systems were minimized for 10,000 integration steps and equilibrated for 20 ns with 1 fs time steps. Following this, a 20 ns unconstrained equilibration production run was performed for subsequent trajectory analysis, with frames stored after each picosecond (every 1000 time steps). Constant temperature (T = 300 K) was enforced using Langevin dynamics with a damping coefficient of 1 ps−1. Constant pressure (p = 1 atm) was enforced using the Nose´-Hoover Langevin piston method with a decay period of 100 fs and a damping time constant of 50 fs. Van der Waals interaction cut-off distances were set at 12Å, (smooth switching function beginning at 10 Å) and long-range electrostatic forces were computed using the particle-mesh Ewald (PME) summation method with a grid size of less than 1Å, along with the pencil decomposition protocol where applicable. SMD simulations fixed the center of mass of neuraminidase -carbons and applied a force to the center of mass of oseltamivir, along a vector connecting the two centers of masses. In simSMD1, a constant velocity protocol was employed, with stretching velocity of 0.5Å/ns. For the SMD spring constant we chose kSMD = 3kBT/Å2 which corresponds to an RMSD value of (pkBT/kSMD)1/2 < 0.6 Å.

### **4.3. Results**

174 Bioinformatics

binding pocket is shown.

**4.2. Computational details** 

sizes of approximately 35,000 atoms.

those mutations on the kinetics of drug binding. Since both H274Y and N294S are nonactive site mutations, prior studies which focused on endpoint interactions between drug and proteins have been unable to provide a full understanding of how these mutations directly impact drug binding. In fact, given that the drug binding kinetics of H5N1 mutants are significantly diminished, it is possible that these mutations alter the binding process, and not necessarily just the specific endpoint interactions. It is well known that electrostatic surface potential of a protein can be an important driving force directing the diffusion of ligands into a protein's active site. The resulting electrostatic maps shown in Figure 7A for H5N1 and Figure 7B for H1N1pdm reveal a highly negatively charged column of residues that forms a pathway 10 Å in length between the SA binding site and the edge of the binding cavity mouth. Electrostatic calculations also reveal that oseltamivir has a highly positive electrostatic surface potential, illustrated in Figure 7C. The question arises whether the negatively charged surface column plays a role in the binding and unbinding of oseltamivir, given a possible mutual attraction between oseltamivir and this column. To answer this question, we employed SMD simulations (described in the Methods section) to pull oseltamivir out of the SA binding site

**Figure 7.** Electrostatic surface potential of the SA binding pocket of H1N1pdm and oseltamivir. Shown in A) and B) are closeup views of the SA binding pocket with drug bound H1N1pdm and avian H5N1 neuraminidase, respectively. The region of the binding pocket where the drug bound possesses a highly negative potential (colored red), whereas the opening of the pocket is surrounded by a highly positive potential ring (colored blue). In C), a detailed surface electrostatic potential for oseltamivir. Shown are the "front" side facing the annulus of the binding pocket, and "back" side facing the interior of the

Six systems were modelled and simulated for oseltamivir bound H5N1, and H1N1pdm wild type and H274Y and N294S mutants. In simSMD1, steered molecular dynamics (SMD) simulations were used to remove oseltamivir from its stable binding site in H5N1 neuraminidase. In simFEQ1-10, equilibration simulations used a starting point generated from simSMD1 in which Oseltamivir has undergone an axial rotation such that it is partially displaced from its binding site. In total, 680 ns of simulations were carried out on system

The "Ensemble" column lists the variables held constant during the simulations; N, p, T, and V correspond to the number of atoms, pressure, temperature, and volume respectively.

and probe possible unbinding pathways. Such simulations have

To simulate the binding of the drug would be the most natural approach, but the computations would be rather expensive. It is possible to the put drug in front of the receptor active site and the attraction forces between the drug and the receptor will drive it to the binding site. This process usually requires microsecond simulations and there is a high chance that the drug will fluctuate around the water box. Unbinding simulations however are feasible and may also reveal features that are characteristic for binding. In simSMD1, a pulling force was applied to rupture all the stabilizing hydrogen bonds between H5N1 and oseltamivir, in order to draw the drug away from the SA binding site. The results of simSMD1 show that the behaviour of oseltamivir under effect of a pulling force can be divided into three distinct stages: 1) from 0 to 8 ns, a buildup of forces during which hydrogen bonds between oseltamivir with E119, D151 and R152 are destabilized; 2) at 8 ns where the remaining stable hydrogen bonds with R292 and R371 rupture; 3) after 8 ns during which the drug is pulled out of the binding pocket. Figure 8 shows the force dependent rupture of hydrogen bonds in simSMD1 by plotting both hydrogen bond stability and (in inset) the force vs. time curve. To our surprise, it was observed that

Incorporating Molecular Dynamics Simulations

into Rational Drug Design: A Case Study on Influenza a Neuraminidases 177

therefore performed additional equilibrium simulations (simFEQ1-10) with oseltamivir already in this flipped state. From these simulations, we were able to observe two distinct outcomes: 1) oseltamivir is able to escape the binding pocket by interacting favourably with the charged binding funnel and 2) the drug returns, not unexpectedly, to its stably bound "preflip" state. Each simulation for FEQ1-10 was run long enough to observe either outcome, with the exception of FEQ5. SimFEQ5 was a special case in which the drug, after following the binding funnel to escape the protein, actually rebound to the SA active site through the same binding funnel after we extended the simulation to follow its movement. A summary of observed outcomes from these simulations is shown in Table 2. In simFEQ1- 5, oseltamivir successfully escaped the SA binding site, whereas in simFEQ6-10, oseltamivir returned to its stably bound "preflip" state. Oseltamivir was observed to diffuse out of the SA active site after strong interaction with the electrostatically charged binding funnel (described above) in five out of ten equilibrium simulations we performed (simFEQ1-5).

**Figure 9.** Forced unbinding of oseltamivir from H5N1 neuraminidase. Shown here are the relative positions of oseltamivir on the electrostatic surface of avian H5N1 neuraminidase during simSMD1.At 0 ns (A), oseltamivir was stably bound within the active site, as seen in simEQ1. Application of force ruptured the stabilizing hydrogen bonds between H5N1 and oseltamivir (see Figure 9), drawing the drug away from its stable binding site within 10 ns, as shown in B. Over the next 2.5 ns of pulling, oseltamivir followed the charged binding funnel (shown in C) until it was completely free of the protein

In four of the cases (simFEQ1-3, 5) oseltamivir was observed to diffuse along the full length of the binding funnel before disassociating with neuraminidase. In simFEQ5 we observed not only a diffusion of oseltamivir through the charged binding funnel, but also the re-entry.

Snapshots from each of these events in FEQ5 are shown in Figure 10. Analysis of interactions of the newly rebound oseltamivir with active site residues from 50 to 100 ns revealed that the drug was stabilized by hydrogen bonds with Y406, R292, D151, E119, and R118, even though the pentyl group had not yet moved to its requisite hydrophobic pocket

binding pocket after 15 ns, as shown in D).

(I222-R224- A246-E276).

**Figure 8.** Distances between hydrogen bond acceptor-donor pairs between oseltamivir and active site amino acids vs. simulation time in simSMD1. Most hydrogen bonds were quickly broken by the pulling force except for those with R292 and R371, which fully ruptured only after 8 ns, corresponding to the peak of the curve of the applied force vs. simulation time (shown in inset). Position of oseltamivir and all its possible hydrogen bonding pairs with residues located along this pathway, it was seen that nonspecific electrostatic attractions formed the predominant interactions between drug and protein.

Oseltamivir followed a lateral escape path through strong interaction with the negatively charged column of residues indentified via electrostatic mapping above, despite application of force to pull the drug straight out of the binding pocket. This divergent path taken by the drug is shown in Figure 9A-D. A second notable observation in simSMD1 was that just before rupture of hydrogen bonds between oselvamivir's carboxyl functional group and residues R292 and R371 of the binding pocket, the drug underwent a rather significant rotation due to the destabilization of hydrogen bonds with E119, D151 and R152. It is likely that this rotation or "flip" is crucial for placing oseltamivir in a position which permits it to leave the SA binding site. Despite application of force to pull oseltamivir straight out of its binding pocket, the drug followed a somewhat lateral path through the electrostatically charged funnel. The "preflip" or stably bound state for oseltamivir in H5N1 is shown in Figure 10A1, 10A2, 10A3, and the flipped transition state shown in Figure 10B1 and 10B2, 10B3.

It was observed in simSMD1 that following the transition to the flipped state, very little force was then required to subsequently draw oseltamivir out of its binding pocket in neuramidase. This result suggests that one may be able to probe the unbinding pathway via diffusion, if the starting state consists of oseltamivir already in its flipped orientation. We therefore performed additional equilibrium simulations (simFEQ1-10) with oseltamivir already in this flipped state. From these simulations, we were able to observe two distinct outcomes: 1) oseltamivir is able to escape the binding pocket by interacting favourably with the charged binding funnel and 2) the drug returns, not unexpectedly, to its stably bound "preflip" state. Each simulation for FEQ1-10 was run long enough to observe either outcome, with the exception of FEQ5. SimFEQ5 was a special case in which the drug, after following the binding funnel to escape the protein, actually rebound to the SA active site through the same binding funnel after we extended the simulation to follow its movement. A summary of observed outcomes from these simulations is shown in Table 2. In simFEQ1- 5, oseltamivir successfully escaped the SA binding site, whereas in simFEQ6-10, oseltamivir returned to its stably bound "preflip" state. Oseltamivir was observed to diffuse out of the SA active site after strong interaction with the electrostatically charged binding funnel (described above) in five out of ten equilibrium simulations we performed (simFEQ1-5).

176 Bioinformatics

10B3.

**Figure 8.** Distances between hydrogen bond acceptor-donor pairs between oseltamivir and active site amino acids vs. simulation time in simSMD1. Most hydrogen bonds were quickly broken by the pulling force except for those with R292 and R371, which fully ruptured only after 8 ns, corresponding to the peak of the curve of the applied force vs. simulation time (shown in inset). Position of oseltamivir and all its possible hydrogen bonding pairs with residues located along this pathway, it was seen that nonspecific electrostatic attractions formed the predominant interactions between drug and protein.

Oseltamivir followed a lateral escape path through strong interaction with the negatively charged column of residues indentified via electrostatic mapping above, despite application of force to pull the drug straight out of the binding pocket. This divergent path taken by the drug is shown in Figure 9A-D. A second notable observation in simSMD1 was that just before rupture of hydrogen bonds between oselvamivir's carboxyl functional group and residues R292 and R371 of the binding pocket, the drug underwent a rather significant rotation due to the destabilization of hydrogen bonds with E119, D151 and R152. It is likely that this rotation or "flip" is crucial for placing oseltamivir in a position which permits it to leave the SA binding site. Despite application of force to pull oseltamivir straight out of its binding pocket, the drug followed a somewhat lateral path through the electrostatically charged funnel. The "preflip" or stably bound state for oseltamivir in H5N1 is shown in Figure 10A1, 10A2, 10A3, and the flipped transition state shown in Figure 10B1 and 10B2,

It was observed in simSMD1 that following the transition to the flipped state, very little force was then required to subsequently draw oseltamivir out of its binding pocket in neuramidase. This result suggests that one may be able to probe the unbinding pathway via diffusion, if the starting state consists of oseltamivir already in its flipped orientation. We

**Figure 9.** Forced unbinding of oseltamivir from H5N1 neuraminidase. Shown here are the relative positions of oseltamivir on the electrostatic surface of avian H5N1 neuraminidase during simSMD1.At 0 ns (A), oseltamivir was stably bound within the active site, as seen in simEQ1. Application of force ruptured the stabilizing hydrogen bonds between H5N1 and oseltamivir (see Figure 9), drawing the drug away from its stable binding site within 10 ns, as shown in B. Over the next 2.5 ns of pulling, oseltamivir followed the charged binding funnel (shown in C) until it was completely free of the protein binding pocket after 15 ns, as shown in D).

In four of the cases (simFEQ1-3, 5) oseltamivir was observed to diffuse along the full length of the binding funnel before disassociating with neuraminidase. In simFEQ5 we observed not only a diffusion of oseltamivir through the charged binding funnel, but also the re-entry.

Snapshots from each of these events in FEQ5 are shown in Figure 10. Analysis of interactions of the newly rebound oseltamivir with active site residues from 50 to 100 ns revealed that the drug was stabilized by hydrogen bonds with Y406, R292, D151, E119, and R118, even though the pentyl group had not yet moved to its requisite hydrophobic pocket (I222-R224- A246-E276).

Incorporating Molecular Dynamics Simulations

into Rational Drug Design: A Case Study on Influenza a Neuraminidases 179

**Figure 11.** Escape and rebinding of oseltamivir through the electrostatic binding funnel in H5N1 neuraminidase during simFEQ5. Shown here are snapshots of simFEQ5, in which oseltamivir first diffused out of the SA active site through interaction with the electrostatic binding funnel within the first 25ns of simulation. Between 28 and 35ns, oseltamivir diffuses and approaches the periphery of the binding pocket away from the binding funnel, but is prohibited from entering due to electrostatic repulsion (45ns). However, between 45 and 50ns, oseltamivir was observed to approach and enter the neuraminidase binding pocket through the binding funnel, adopting a stable position within the SA

The result has shed light on the important role of the electrostatic surface potentials in directing the diffusion of oseltamivir into the SA binding site. From our simulations, it is clear that the negatively charged funnel serves as a prominent binding and unbinding pathway for oseltamivir in the wild type systems investigated with SMD and followup equilibrium simulations in simSMD1 and simFEQ1-10. It turns out also that the binding funnel may possibly play a crucial role in drug resistance caused by mutations. The conspicuous location of residue 294, which maps directly onto this negatively charged pathway, may play a key role in the N294S mutation for disrupting the proper guidance of the drug into its binding pocket. Thus, it is possible that the 274 and 294 mutations may confer drug resistance by not only disrupting the end-point interactions of oseltamivir but also its entry into the binding pocket by interfering with the binding funnel. While sources

binding pocket through hydrogen bonds with Y406, R292, D151, E119, and R118.

**Figure 10.** Oseltamivir in "flipped" position (position of oseltamivir at 7.5ns in simSMD1) in comparison with its stable equilibrium position,


**Table 2.** Summary of FEQ1-10 simulations starting from "flipped" position of oseltamivir taken from simSMD1 at 7.5 ns.

This image cannot currently be displayed.

**Figure 10.** Oseltamivir in "flipped" position (position of oseltamivir at 7.5ns in simSMD1) in

simFEQ1 Drug escape via binding funnel 15

simFEQ2 Drug escape via binding funnel 10

simFEQ3 Drug escape via binding funnel 15

simFEQ4 Drug interaction with binding funnel but escape via 430-cavity 50

simFEQ5 Drug escape and rebind into SA pocket via binding funnel 100

simFEQ6 Drug returned to "preflip" position 10

simFEQ7 Drug returned to "preflip" position 15

simFEQ8 Drug returned to "preflip" position 50

simFEQ9 Drug returned to "preflip" position 50

simFEQ10 Drug returned to "preflip" position 50

**Table 2.** Summary of FEQ1-10 simulations starting from "flipped" position of oseltamivir taken from

Name Result Time (ns)

comparison with its stable equilibrium position,

simSMD1 at 7.5 ns.

**Figure 11.** Escape and rebinding of oseltamivir through the electrostatic binding funnel in H5N1 neuraminidase during simFEQ5. Shown here are snapshots of simFEQ5, in which oseltamivir first diffused out of the SA active site through interaction with the electrostatic binding funnel within the first 25ns of simulation. Between 28 and 35ns, oseltamivir diffuses and approaches the periphery of the binding pocket away from the binding funnel, but is prohibited from entering due to electrostatic repulsion (45ns). However, between 45 and 50ns, oseltamivir was observed to approach and enter the neuraminidase binding pocket through the binding funnel, adopting a stable position within the SA binding pocket through hydrogen bonds with Y406, R292, D151, E119, and R118.

The result has shed light on the important role of the electrostatic surface potentials in directing the diffusion of oseltamivir into the SA binding site. From our simulations, it is clear that the negatively charged funnel serves as a prominent binding and unbinding pathway for oseltamivir in the wild type systems investigated with SMD and followup equilibrium simulations in simSMD1 and simFEQ1-10. It turns out also that the binding funnel may possibly play a crucial role in drug resistance caused by mutations. The conspicuous location of residue 294, which maps directly onto this negatively charged pathway, may play a key role in the N294S mutation for disrupting the proper guidance of the drug into its binding pocket. Thus, it is possible that the 274 and 294 mutations may confer drug resistance by not only disrupting the end-point interactions of oseltamivir but also its entry into the binding pocket by interfering with the binding funnel. While sources of end-point interactions, including hydrogen bonds, hydrophobic packing and solvent infiltration, of oseltamivir resistance have been thoroughly studied, little is known about the kinetics of the drug binding in mutants at atomic level. The idea that drug resistant mutants actually disrupt entry of oseltamivir into the SA active site of neuraminidase through disruption of an electrostatic binding funnel is in part supported by experiments which have noted reduced drug binding kinetics in H5N1 H274Y and N294S mutants. Even though the oseltamivir-resistant mutations were seen located in or adjacent to the funnel, clearly additional study is still needed for a full understanding of how the H274Y and N294S mutations weaken the binding of the drug. Furthermore, the mutations might not only affect the electrostatic gradient but also the geometry of the drug binding and unbinding path. Future or followup studies should therefore focus on the specific drug entry/exit pathways for oseltamivir-resistant mutant systems, sampling timescales great enough (such as in the case of simFEQ5 described above) to capture both the binding.

Incorporating Molecular Dynamics Simulations

into Rational Drug Design: A Case Study on Influenza a Neuraminidases 181

The author wish to thank the Institute for Computational Science and Technology at Ho Chi Minh City and NAFOSTED (The National Foundation for Science and Technology

Abed, Y.; Baz, M.; Boivin, G., Impact of neuraminidase mutations conferring influenza resistance to neuraminidase inhibitors in the N1 and N2 genetic backgrounds. *Antiviral* 

Andrew C. Kruse, Jianxin Hu, Albert C. Pan, Daniel H. Arlow, Daniel M. Rosenbaum, Erica Rosemond, Hillary F. Green, Tong Liu, Pil Seok Chae, Ron O. Dror, David E. Shaw, William I. Weis, Jürgen Wess, and Brian K. Kobilka, "Structure and Dynamics of the M3 Muscarinic Acetylcholine Receptor," *Nature*, vol. 482, no. 7386,

Amaro, R. E.; Baron, R.; McCammon, J. A. An improved relaxed complex scheme for receptor flexibility in computer-aided drug design. *J. Comput.-Aided Mol. Des.* 2008, 22

Berman HM, B. T., Bhat TN, Bluhm WF, Bourne PE, Burkhardt K, Feng Z, Gilliland GL, Iype L, Jain S, Fagan P, Marvin J, Padilla D, Ravichandran V, Schneider B, Thanki N, Weissig H, Westbrook JD, Zardecki C., The Protein Data Bank, Acta Crystallogr D Biol

Brooks BR, Bruccoleri RE, Olafson BD, States DJ, Swaminathan S, Karplus M: CHARMM - a program for macromolecular energy, minimization, and dynamics calculations. *J* 

Case DA, Cheatham TE, Darden T, Gohlke H, Luo R, Merz KM Jr, Onufriev A, Simmerling C, Wang B, Woods RJ: The AMBER biomolecular simulation programs. J Comput Chem

Cheng, L. S.; Amaro, R. E.; Xu, D.; Li, W. W.; Arzberger, P. W.; McCammon, J. A., Ensemble-Based Virtual Screening Reveals Potential Novel Antiviral Compounds for Avian

Christen M, Hünenberger PH, Bakowies D, Baron R, Bürgi R, Geerke DP, Heinz TN, Kastenholz MA, Kräutler V, Oostenbrink C, Peter C, Trzesniak D, van Gunsteren WF: The GROMOS software for biomolecular simulation: GROMOS05. *J Comput Chem* 2005,

Collins, P. J.; Haire, L. F.; Lin, Y. P.; Liu, J.; Russell, R. J.; Walker, P. A.; Skehel, J. J.; Martin, S. R.; Hay, A. J.; Gamblin, S. J., Crystal structures of oseltamivir-resistant influenza virus

neuraminidase mutants. *Nature* (London, U. K.) 2008, 453 (7199), 1258-1261.

Development) for partially support the research presented in this chapter.

Butler, D., Swine flu goes global. *Nature*. 2009, (458), 1082-1083.

Crystallogr. 2002 Jun;58(Pt 6 No 1):899-907. Epub 2002 May 29.

Influenza Neuraminidase. *J. Med. Chem.* 2008, 51 (13), 3878-3894.

**Acknowledgement** 

*Ther*. 2006, 11 (8), 971-976.

*Comput Chem* 1983, 4:187-217.

2005, 26:1668-1688

26:1719-1751.

2012, pp. 552–556

(9), 693-705.

**6. References** 
