Protein Simulation with CHARMM-GUI and NAMD

Protein simulation systems were prepared with the CHARMM-GUI.^{28 (link)} Briefly, protein structures taken from corresponding protein data bank^{29 (link)} files were solvated in pre-equilibrated cubic TIP3P water boxes of suitable sizes and counter-ions were added to keep systems neutral as detailed in Table 1. Periodic boundary conditions were applied and Lennard-Jones (LJ) interactions were truncated at 12 Å with a force switch smoothing function from 10 Å to 12 Å. The non-bonded interaction lists were generated with a distance cutoff of 16 Å and updated heuristically. Electrostatic interactions were calculated using the particle mesh Ewald method³⁰ with a real space cutoff of 12 Å on an approximately 1 Å grid with 6th order spline. Covalent bonds to hydrogen atoms were constrained by SHAKE.³¹ After a 200 step Steepest Descent (SD) minimization with the protein fixed and another 200 steps without the protein fixed, the systems were first heated to 300 K and then subjected to a 100 ps NVT simulation followed by a 100 ps NPT simulation. The minimization, heating and initial equilibrium was performed with CHARMM,^{32 (link)} and the resultant structures were used to start simulations in NAMD.^{33 (link)} After a 1 ns NPT simulation as equilibration, the production simulations were run for 100 ns in the NVT ensemble (see Table 1). For HEWL NPT ensembles were generated to better compare with previous work that found CMAP helps to better reproduced order parameter S^{2,34 (link)} and simulations were extended to 200 ns to reduce the uncertainty of the computed S². Langevin thermostat with a damping factor of 5 ps⁻¹ was used for NVT simulation and the Nosé-Hoover Langevin piston method with a barostat oscillation time scale of 200 fs was further applied for the NPT simulation at 300 K and 1 atm. The time step equals 2 fs and coordinates were stored every 10 ps. For each protein the above simulation protocol was applied with the C36 and C22/CMAP FFs, while for ubiquitin an additional 1.2 μs trajectories with C36 was generated. This long simulation is used to check the convergence and also to examine whether computed NMR data deteriorate over a longer simulation time, as it was reported that RDCs significantly deviate from experimental values after approximately 500 ns simulations with the C22 FF.^{22 (link)}