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

Similar protocols

Protocol publication

[…] N1 monomer simulations were performed using the GROMOS05 software for biomolecular simulation and the GROMOS96 force field (45A3 parameter set). Molecular Dynamics (MD) set up is described in detail in the Supporting Information Methods section. Both complexes of N1-oseltamivir, with and without the bound calcium ion, were simulated in 10 independent trajectories, each 4 ns, for a total of 40 ns of simulation for each complex. Comprising the 10 simulations were five simulations generated from the chain B monomer in the oseltamivir-bound 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). 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; the calcium was parametrized in the classical force field. To generate the independent trajectories, these structures were each initialized with random velocities assigned from a Maxwell–Boltzmann distribution at 5 K.Two 100 ns N1 tetramer simulations with the AMBER FF99SB force field, were performed, each with atomic coordinates taken from the holo, open Loop 150 crystal structure (2HU0) and with the calcium inserted from overlap with the apo 2HTY structure. Calcium-bound simulations used the PMEMD module in AMBER 10, whereas calcium-free simulations were performed using the Desmond Molecular Dynamics package developed by D.E.Shaw Research. For details of these simulations, please see the Supporting Information Methods section. Although there are a few differences in the specifics of the MD engines used to run the AMBER FF99SB simulations (Supporting Information Methods), the similarity in RMSF (Supporting Information Fig. 1) indicates similar sampling of conformational space. AMBER FF99SB trajectories for each monomer of the tetramer were extracted and concatenated to approximate 400 ns of monomer N1 sampling.Analysis tools in Visual Molecular Dynamics (VMD) were applied for calculation of root mean squared fluctuation (RMSF) and deviation (RMSD), and monitoring of torsion, hydrogen bonds, and salt bridge distances. Hydrogen bond criteria applied a donor–acceptor distance cutoff of 3.5 Å and a 120° cutoff for the donor-hydrogen-acceptor angle. For analysis, structures were extracted every 10 ps for each 400 ns AMBER FF99SB trajectory and 2 ps for the 40 ns GROMOS96 trajectories. These structures were used for the GROMOS RMSD clustering algorithm, applied using the cluster2 program in GROMOS05 software. A rotational and translational fitting was applied to all the Cα carbons, followed by application of the RMSD clustering to all atoms of 41 residues that comprise the binding site (residues 117–119, 133–138, 156, 178–179, 196–200, 223–228, 243–247, 276, 277, 292–295, 344–347, 371, 401, 402). This selection is similar to that used for relaxed complex docking of N1 in Cheng et al., but here we exclude residues in the highly flexible Loop 150 (residues 147–152) and Loop 430 (430–439) regions to focus the algorithm on the active site portion near the calcium binding site. After screening a variety of RMSD cutoff values for cluster generation, using the total number and diversity of clusters as criteria, a value of 1.5 Å was selected for the GROMOS96 trajectories, and a value of 1.2 Å for AMBER FF99SB trajectories. The discrepancy in the chosen RMSD cutoff, as well as the varied cluster populations, can be attributed to the varied potential energy landscapes produced by the different force fields, as well as the use of multiple, shorter GROMOS96 trajectories and fewer, longer AMBER FF99SB simulations.Two different, rigorous binding free energy calculations were used: free energy perturbation (FEP), with the AMBER FF99SB force field, and thermodynamic integration (TI) with the GROMOS96 force field. Both bound and unbound calculations were performed, to give a ΔGprotein from decoupling the bound oseltamivir from N1, and a ΔGwater from decoupling of unbound oseltamivir, free in a box of water. Following the thermodynamic cycle involved in double decoupling,, ΔGprotein is subtracted from ΔGwater to give ΔGbind.The GROMOS05 software was used for both oseltamivir-bound and oseltamivir-unbound TI calculations, repeated in independent trajectories (IT-TI). For the bound state decoupling, six calculations decoupled oseltamivir from the ion-bound complex, and six decoupled the ligand from the ion-free complex, for a total of 12 N1 TI calculations. For each of these complexes, the calculations were performed in triplicate, with different starting structures that capture diverse conformations of Loop 150; three calculations were initialized from a closed Loop 150 N1 conformation, and three were initialized from an open Loop 150 conformation. These starting structures were extracted from the independent MD trajectories described above, after 1 ns of equilibration. Six independent oseltamivir-unbound TI calculations decoupled oseltamivir from a 37 × 37 × 37 Å3 box of SPC water. For each TI, 26 windows span the λ range from 0.00 (ligand coupled) to 1.00 (ligand decoupled), in 0.04 intervals. For ligand–protein decoupling, a 500 ps equilibration period preceded a 500 ps sampling period, from which data was recorded; the unbound decoupling calculations involved 400 ps for both equilibration and sampling periods per λ window.To maintain a defined bound state in the protein calculations, a 0.6 kcal/mol/Å2 harmonic restraint constant was applied to the C4 atom of oseltamivir and the associated volume correction + RTln(CV) for the free energy was calculated to be −5.5 kcal/mol. Soft-core interaction potentials (sLJ = 0.5 and sC = 0.5) were used for ligand atoms involved in the perturbation to avoid end point singularities and to enhance sampling. For each λ window, the statistical inefficiency of the raw data was calculated using the method described and coded by Chodera et al. The decorrelated data was averaged to approximate , which was integrated by the trapezoidal rule over the 26 λ windows for each of N= 6 independent calculations, as: (1) Subtraction of the ΔGprotein from ΔGwater gives a ΔGbind. A standard error was computed for both ΔGprotein and ΔGwater over the N = 6 independent runs and propagated for ΔGbind (listed in ).To examine the advantage of shorter, independent trajectories over a single, long trajectory in the TI calculations, one run for each set of protein TI calculations was selected for the extension of equilibration and sampling times (1 ns equilibration and 3.8 ns sampling per λ), for an 85.8 ns increase in total sampling time for the calculation. The single ΔGprotein produced from each extended run differed from that of the “short” calculation result by <4%; however, the associated ΔGbind was +3.3 kcal/mol less favorable than the experimentally derived target value. This discrepancy is larger than that for the IT-TI result for ΔGbind, which is off the target value by +1.3 kcal/mol (). The extended ion-free result was also significantly more unfavorable than the IT-TI value. This indicates that these single, extended simulations were unable to overcome sampling barriers to generate a complete ensemble average . Combining multiple TI calculations initialized from diverse conformations improves sampling and can result in a more correct ΔGbind., , Free energy calculations on the tetrameric complex were performed with the AMBER FF99SB force field in Desmond using FEP. For both the decoupling of unbound oseltamivir (to calculate ΔGwater) and bound oseltamivir (for ΔGprotein), 21 windows were used, one with full non-bonded interaction, 10 for the annihilation of columbic interactions (with full van der Waals interactions), and 10 with the removal of van der Waals interactions (with no electrostatics). A softcore potential was used with α = 0.5. Decoupling of unbound oseltamivir in water was performed on a system of the ligand in a 24 × 28 × 29 Å3 TIP3P water box with 250 ps of equilibration followed by 1 ns of sampling in each window. Bound oseltamivir decoupling calculations were performed on the tetrameric complex (from 2HU0, described above) with all four ligand molecules simultaneously removed to enhance sampling. Following 5 ns of equilibration of the starting structure, each λ window was individually equilibrated for 1 ns and then sampled for 2 ns.Three sets of protein–ligand calculations were performed, two with the protein coordinates from the 2HU0 crystal structure (one with the bound calcium and one without) and another without calcium and initiated from clustering results in which Y347 is flipped out of the pocket (see below). Reported work functions were post-processed with the Multistate Bennett Acceptance Ratio (MBAR). Energies were decorrelated based on their statistical inefficiencies, as in the GROMOS96 TI simulations, and error bars representing the analytic errors developed in the MBAR implementation, here represented by σ, were calculated with code made publicly available by Shirts and Chodera. The free energy of decoupling the four ligands from the N1 tetramer was divided by four to obtain an average binding free energy to the N1 monomer, with a standard error calculated using N = 4. […]

Pipeline specifications

Software tools GROMOS, VMD
Application Protein structure analysis
Diseases Neural Tube Defects
Chemicals Calcium, Amber, Oseltamivir