Impact of calcium on N1 influenza neuraminidase dynamics and binding free energy

The highly pathogenic influenza strains H5N1 and H1N1 are currently treated with inhibitors of the viral surface protein neuraminidase (N1). Crystal structures of N1 indicate a conserved, high affinity calcium binding site located near the active site. The specific role of this calcium in the enzyme mechanism is unknown, though it has been shown to be important for enzymatic activity and thermostability. We report molecular dynamics (MD) simulations of calcium-bound and calcium-free N1 complexes with the inhibitor oseltamivir (marketed as the drug Tamiflu), independently using both the AMBER FF99SB and GROMOS96 force fields, to give structural insight into calcium stabilization of key framework residues. Y347, which demonstrates similar sampling patterns in the simulations of both force fields, is implicated as an important N1 residue that can “clamp” the ligand into a favorable binding pose. Free energy perturbation and thermodynamic integration calculations, using two different force fields, support the importance of Y347 and indicate a +3 to +5 kcal/mol change in the binding free energy of oseltamivir in the absence of calcium. With the important role of structure-based drug design for neuraminidase inhibitors and the growing literature on emerging strains and subtypes, inclusion of this calcium for active site stability is particularly crucial for computational efforts such as homology modeling, virtual screening, and free energy methods. Proteins 2010. © 2010 Wiley-Liss, Inc.


Methods
Both calcium-bound and calcium-free N1 monomer simulations were performed using the GRO-MOS05 software for biomolecular simulation 1 and the GROMOS96 force field (45A3 parameter set). 2 Parameters for oseltamivir were derived from existing building blocks 3 (Table S2). Amino acid charges were defined to reproduce an apparent pH 7. The systems were solvated in boxes of SPC water molecules, with a 12 Å barrier to the periodic boundary of the cube, and neutralized with sodium ions. For the ion-bound simulations, the calcium was parametrized in the classical force field, with all ion parameters taken from those developed by Åqvist. 4 After 2,000 steps of steepest descent energy minimization, the system was then brought to the reference temperature of 300 K in six consecutive 25 ps MD periods (50 K increments), gradually Both monomer complexes of N1-oseltamivir, with and without the bound calcium ion, were simulated in ten independent trajectories, each 4 ns, for a total of 40 ns of simulation for each complex. Comprising the ten simulations were five simulations generated from the chain B monomer in the holo Loop 150 "open" crystal structure (PDBID: 2HU0) and five simulations started from the chain A monomer in the holo Loop 150 "closed" crystal structure (PDBID: 2HU4). 10 As the calcium density was not present in these crystal structures, overlap with the apo N1 structure 2HTY aided in positioning of the ion in the protein. After minimization, these structures were each initialized with random velocities assigned from a Maxwell-Boltzmann distribution at 5 K to generate the independent trajectories.
Calcium-bound 100 ns N1 tetramer simulations were performed with the PMEMD module in AMBER 10 11 and the AMBER FF99SB force field. 12 Atomic coordinates were taken from the holo, open Loop 150 crystal structure (2HU0), with the calcium inserted from overlap with the apo 2HTY structure and parameterized in the classical force field. Protonation states for histidines and other titratable groups were determined at pH 6.5 by the PDB2PQR 13 web server and manually verified. The tetrameric 2HU0 crystal structure has a single oseltamivir molecule bound in the active site of chain B. In order to introduce the oseltamivir within each of the other chains, chain B was aligned to chain A, C, and D using VMD 14 and the resulting transformation matrix was applied to the oseltamivir molecule. Oseltamivir was parameterized according to quantum chemical calculations, which included performing a geometry optimization with Gaussian03 15 at the Hartree-Fock/6-31G* level. The resulting atomic partial charges were then determined according to the RESP method, 16 and the atom types were assigned by the Antechamber module of AMBERTools 1.2. The GAFF 17 force field within AMBER was employed to generate the bond, angle, and dihedral parameters. As no water molecules were reported in the 2HU0 structure, we structurally aligned the 2HTY and 2HU0 systems and kept all crystallographic water molecules that did not clash with oseltamivir in the binding pocket. The system was built using the AMBER9 program Leap and the Amber99SB force field. Each monomer chain contained 8 disulfide bonds, which were properly enforced using the CYX notation in AMBER. A box of TIP3P 18 waters was added to solvate each system, resulting in a rectangular box of dimensions 124 x 127 x 77 Å 3 . The system was neutralized with the addition of sodium counter ions and a 150 mM NaCl salt bath was introduced.
The constructed N1 ion-bound tetramer complex was first subjected to 2000 steps of steep-est descent, followed by 5000 steps of conjugate gradient minimization using 5 kcal·mol −1 · Å −2 harmonic restraints on all non-hydrogen protein atoms. Then, 5000 steps of conjugate gradient minimization with just the backbone atoms restrained cleaned up the initial hydrogenated complex. A further 25,000 conjugate gradient minimization steps were then performed on the entire complex, without restraints, in order to alleviate any steric clashes prior to performing molecular dynamics. Following minimization, the system was linearly heated to 310 K in the NVT ensemble using a Langevin thermostat, with a collision frequency of 1.0 ps −1 , and harmonic restraints of 4 kcal·mol −1 · Å −2 on the backbone atoms. Further, three 250 ps periods were run in the NPT ensemble with the restraint force constant being reduced by 1 kcal·mol −1 · Å −2 each time. A final 250 ps of NPT dynamics was run without restraints. Production runs were then made for 100 ns duration in the NVT ensemble; temperature was controlled with a Langevin thermostat (1.0 ps −1 collision frequency), and pressure was controlled using a Berendsen Barostat 5 with a coupling constant of 1 ps and a target pressure of 1 atm. The time step used was 2 fs and all hydrogen atoms were constrained using the SHAKE algorithm. 6 Long range electrostatics were included using the Particle Mesh Ewald algorithm 19 with a 4th order B-spline interpolation and a grid spacing of <1.0 Å and a direct space cutoff of 8 Å. The trajectories for each monomer of the tetramer were extracted and concatenated to approximate 400 ns of monomer N1 sampling.
The ion-free tetramer system was simulated with atomic coordinates and parameters identical to that described above for the AMBER10 simulations, except for calcium presence, and using the AMBER ff99SB force field and the Desmond MD engine developed by D. E. Shaw Research. 20 The Maestro modeling suite was utilized for system construction in an 128 x 130 x 80 Å 3 orthorhombic box, for a minimum distance of 12 Å between protein heavy atoms and box edges.
Sodium and chloride ions were added to neutralize the system charge and create an approximately 150 mM NaCl solution, as in the AMBER10 simulation. Following 10,000 steps of steepest-decent minimization, the systems were equilibrated with restraints of 10 kcal· mol −1 ·Å −2 for 50 ps followed by 200 ps in which the restraints were continuously and slowly removed; then, unrestrained molecular dynamics were performed for 99.75 ns. Numerical integration was performed with a 2 fs