Electrostatics

Initial helical conformations were defined as all amino acids having (φ, ψ)=(−60°, −40°). Initial extended conformations were defined as all (φ, ψ)=(180°, 180°). Native conformations, as appropriate, were defined for each system as below. Explicit solvation was achieved with truncated octahedra of TIP3P water¹⁶ with a minimum 8.0 Å buffer between solute atoms and box boundary. All structures were built via the LEaP module of Ambertools. Except where otherwise indicated, equilibration was performed with a weak-coupling (Berendsen) thermostat³³ and barostat targeted to 1 bar with isotropic position scaling as follows. With 100 kcal mol⁻¹ Å⁻² positional restraints on protein heavy atoms, structures were minimized for up to 10000 cycles and then heated at constant volume from 100 K to 300 K over 100 ps, followed by another 100 ps at 300 K. The pressure was equilibrated for 100 ps and then 250 ps with time constants of 100 fs and then 500 fs on coupling of pressure and temperature to 1 bar and 300 K, and 100 kcal mol⁻¹ Å⁻² and then 10 kcal mol⁻¹ Å⁻² Cartesian positional restraints on protein heavy atoms. The system was again minimized, with 10 kcal mol⁻¹ Å⁻² force constant Cartesian restraints on only the protein main chain N, C_α, and C for up to 10000 cycles. Three 100 ps simulations with temperature and pressure time constants of 500 fs were performed, with backbone restraints of 10 kcal mol⁻¹ Å⁻², 1 kcal mol⁻¹ Å⁻², and then 0.1 kcal mol⁻¹ Å⁻². Finally, the system was simulated unrestrained with pressure and temperature time constants of 1 ps for 500 ps with a 2 fs time step, removing center-of-mass translation and rotation every picosecond.
SHAKE³⁴ was performed on all bonds including hydrogen with the AMBER default tolerance of 10⁻⁵ Å for NVT and 10⁻⁶ Å for NVE. Non-bonded interactions were calculated directly up to 8 Å. Beyond 8 Å, electrostatic interactions were treated with cubic spline switching and the particle-mesh Ewald approximation³⁵ in explicit solvent, with direct sum tolerances of 10⁻⁵ for NVT or 10⁻⁶ for NVE. A continuum model correction for energy and pressure was applied to long-range van der Waals interactions. The production timesteps were 2 fs for NVT and 1 fs for NVE.

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)}

One approach for simulating a small part of a large system (e.g.,
the enzyme active site region of a large protein) uses a solvent boundary
potential (SBP). In SBP simulations, the macromolecular system is separated
into an inner and an outer region. In the outer region, part of the
macromolecule may be included explicitly in a fixed configuration, while the
solvent is represented implicitly as a continuous medium. In the inner
region, the solvent molecules and all or part of the macromolecule are
included explicitly and are allowed to move using molecular or stochastic
dynamics. The SBP aims to “mimic” the average
influence of the surroundings, which are not included explicitly in the
simulation.27 ,28 There are several implementations of the SBP
method in CHARMM. The earliest implementation, called the stochastic
boundary potential (SBOU), uses a soft nonpolar restraining potential to
help maintain a constant solvent density in the inner or
“simulation” region while the molecules in a shell
or buffer region are propagated using Langevin dynamics.27 By virtue of its simplicity, this treatment
remains attractive and it is sufficient for many applications.320 (link),321 (link) To improve the treatment of systems with irregular
boundaries in which part of the protein is in the outer region, a refinement
of the method has been developed that first scales the exposed charges to
account for solvent shielding and then corrects for the scaling by
post-processing.307
The Spherical Solvent Boundary Potential (SSBP), which is part of
the Miscellaneous Mean Field Potential (MMFP) module (see Section III F), is
designed to simulate a molecular solute completely surrounded by an
isotropic bulk aqueous phase with a spherical boundary.28 In SSBP the radius of the spherical region is
allowed to fluctuate dynamically and the influence of long-range
electrostatic interactions is incorporated by including the dielectric
reaction field response of the solvent.28 ,29 This approach has
been used to study several systems.322 –325
Because SSBP incorporates the long-range electrostatic reaction field
contribution, the method is particularly useful in free energy calculations
that involve introducing charges.322 –325
Like the SBOU charge-scaling method,307 the Generalized Solvent Boundary Potential (GSBP) is
designed for irregular boundaries when part of the protein is outside the
simulation region.29 However, unlike
SBOU, GSBP includes long-range electrostatic effects and reaction fields. In
the GSBP approach, the influence of the outer region is represented in terms
of a solvent-shielded static field and a reaction field expressed in terms
of a basis set expansion of the charge density in the inner region, with the
basis set coefficients corresponding to generalized electrostatic
multipoles.29 ,326 The solvent-shielded static field from the
outer macromolecular atoms and the reaction field matrix representing the
coupling between the generalized multipoles are both invariant with respect
to the configuration of the explicit atoms in the inner region. They are
calculated only once (with the assumption that the size and shape of inner
region does not change during the simulation) using the finite-difference
Poisson-Boltzmann (PB) equation of the PBEQ module. This formulation is an
accurate and computationally efficient hybrid MD/continuum method for
simulating a small region of a large macromolecular system,326 and is also used in QM/MM approaches.281 (link),327 (link)