**4.1. Characteristics of the electrostatic surface potential of N1 neuraminidases shed light for further study on drug binding and drug resistance**

Up to this point, all proposed mechanisms for oseltamivir resistance have focused on the effects of the mutations on the SA binding site but little has been known about the effects of 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 and probe possible unbinding pathways. Such simulations have

Incorporating Molecular Dynamics Simulations

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

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

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

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

water box 23 and ionized by NaCl (0.152M) to mimic physiological conditions.

corresponds to an RMSD value of (pkBT/kSMD)1/2 < 0.6 Å.

**4.3. Results** 

**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 binding pocket is shown.
