Protein Simulation in Crystal Environments

Simulations of the proteins in their crystal environments (Table 1), which were used previously during optimization of the C22/CMAP force field ^{40 (link)}, were performed using CHARMM on full unit cells with added waters and counterions to fill the vacuum space. Once the full unit cell was constructed based on the coordinates in the protein databank, a box of water with dimensions that encompassed the full unit cell was overlaid onto the crystal coordinates while preserving crystal waters, ions, and ligands. Water molecules with oxygen within 2.8 – 4.0 Å of any of the crystallographic non-hydrogen atoms were removed, as described below, as well as those occupying space beyond the full unit cell. To neutralize the total charge of each system, sodium or chloride ions were added to the system at random locations at least 3.0 Å from any crystallographic non-hydrogen atom or previously added ions and 0.5 Å from any water oxygen. Final selection of the water molecule deletion distance was performed by initially applying a 2.8 Å criteria to all systems followed by system equilibration and an NPT production run of 5 ns following which the lattice parameters were analyzed. The deletion distances were then increased and the equilibration and 5 ns production NPT simulation were repeated until the final lattice parameters were in satisfactory agreement with experimental data. The final water deletion distances and unit cell parameters from the full 40 ns production simulations are presented in Table S2 of the SI. For the minimization and MD simulations, electrostatic interactions were treated with PME using a real space cutoff of 10 Å. The LJ interactions were included with force switching from 8 Å to 10 Å, while the list of nonbonded atoms was kept for interatomic distances of up to 14 Å and updated heuristically. Each crystal system was first minimized with 100 steps of steepest-decent (SD) with non-water, non-ion crystallographic atoms held fixed followed by 200 steps of SD with harmonic positional restraints of 5 kcal/mol/Ų on solute non-hydrogen atoms. The minimized system was then subject to an equilibration phase consisting of 100 ps of NVT simulation⁴¹ in the presence of harmonic positional restraints followed by 5 ns (100 ps for 135L and 3ICB) of fully relaxed NVT simulation with a time step of 2 fs. During the simulations all covalent bonds involving hydrogens were constrained using SHAKE⁴² . Production phase simulations were conducted for 40 ns in the isothermal and isobaric NPT ensemble⁴³ . The only symmetry enforced was translational (i.e. periodic boundaries). Reference temperatures were set to match the crystallographic conditions (Table S2) and maintained by the Nosé-Hoover thermostat with a thermal piston mass of 1,000 kcal ps²/mol while a pressure mass of 600 amu was used with the Langevin piston. The first 5 ns of the production simulations were considered as equilibration and therefore discarded from analysis, which was performed on coordinate sets saved every 5 ps. The boundaries for α helices and β strands were obtained from a consensus of author annotations and structural assignments calculated by DSSP^{44 (link)} and STRIDE^{45 (link)} from the crystal structures.