Evaluating Force Matching as a Parametrization Strategy for the CHARMM36m Force Field Using Phosphorylation

1. Characterization of different conformations obtained from clustering for simulations containing 2 MP0 molecules (top left), MP1 and MAM (top right), and MP2 with two MAMs (bottom row, left with the first MAM molecule MAM_1 and right with the second MAM molecule MAM_2). A broad spectrum of orientations of the solutes is obtained and later used in the force matching process. In the middle, examples are shown in ball and stick representation. Hydrogens are white, carbon gray, oxygen red, phosphorus orange, and nitrogen blue. The vectors used for the estimation of the relative angles are depicted as yellow arrows through the molecules. The visualized data points are highlighted in yellow.

2. Exemplary frame of one significant conformation indicated by clustering of MP1 and MAM in a 4 × 4 × 4 nm simulation box containing 2003 water molecules. The corresponding cropped 1.6 × 1.6 × 1.6 nm box, containing 112 water molecules, for quantum chemical (QC) calculation is shown in the middle, and the box with the QC-optimized solutes (PBE/TZV2P-MOLOPT-PBE-GTH/DFTD3) is shown on the right.

3. Exemplary visualization of the QC systems used for the net force calculation. Forces obtained at the PBE/TZV2P-MOLOPT-PBE-GTH/DFTD3 level from the boxes with each solute are subtracted from the forces in the entire box to obtain intermolecular forces between the solutes, here MP1 and MAM. In the case of two counterions to MP, both were subtracted at once.

4. Top left: Ac–K+–pS0–NH2 at pH = 2, i.e., where both lysine (LYS) and phosphoserine (SP0) are fully protonated. Bottom left: Ac–K+–pS––NH2 at pH = 4.5, i.e., where lysine (LYS) is fully protonated and phosphoserine (SP1) singly protonated, i.e., carrying the charge −1. Right: Cubic simulation box including 3 dipeptides Ac–K+–pS––NH2, 7 Na+ (purple dots) and 7 Cl– ions (green dots) and 8789 water molecules (cyan) used for the estimation of the NMR relaxation times from MD simulations.

5. Mapping of the here-defined atom types of MP0, MP1, and MP2 representing the amino acid side-chain analogs of phosphorylated serine residues SP0, SP1, and SP2, respectively. Hydrogens are shown as white spheres, carbons as dark gray spheres, oxygens as red spheres, and phosphorus as orange spheres.

6. Simulation setup for osmotic pressure simulations exemplarily shown for the MP2–MAM 1:2 solution in a total concentration of 1.2 mol/L using the MP2–2MAM vc parameters and the CHARMM TIP3P water model. The solutes are shown in ball + stick representation, hydrogen atoms are colored white, oxygens red, phosphorus atoms orange, and nitrogens blue. The two force boundaries which are acting only on solutes, i.e., on MP2 and MAM, are visualized as green surfaces. The z-axis is located in the left–right direction.

7. Comparison of the osmotic concentration of MP1–MAM solutions in simulations to experimental osmotic concentrations of MPA–MAM 1:1 solutions. Osmotic concentrations are plotted against total molar concentrations of both electrolytes. Left, the effect of different water models, i.e., standard CHARMM TIP3P, TIP4P, TIP4P/ε, and OPC on the osmotic activity using the original CHARMM36m force field parameters (original parameters). Right, comparison of the here developed parameters for MP1 (vc and vpc stand for variance over all conformations and over each conformation, respectively, Na+ or MAM indicates the counterion used in the force-matching procedure) to the original CHARMM36m parameters (all simulated together with the CHARMM TIP3P water model) and to the experimentally obtained values. The error bars show the standard error of the mean over 3 independent replicas.

8. Comparison of the osmotic concentration of a MP2–2MAM solution in simulations compared to experimental values of a MPA-MAM 1:2 solution. Osmotic concentrations are plotted against total molar concentrations of both electrolytes. Left, the effect of different water models, i.e., standard CHARMM TIP3P, TIP4P, TIP4P/ε, and OPC on the osmotic concentration using the original CHARMM36m force field parameters (original parameters). Right, comparison of the here developed parameters for MP2 to the original CHARMM36m parameters (all simulated in the CHARMM TIP3P water model) and to the experimentally obtained values. The error bars show standard error of the mean over 3 independent replicas.

9. Left: Comparison of osmotic concentrations from simulation of MP2–1-1MAM (vc MP2–1-1MAM, corresponding to our final parameters) to the MP2–2MAM simulations using original CHARMM parameters, or our MP2–2MAM vc parameters (vc MAM), and experimental values of MPA-MAM 1:2 solutions. All simulations are performed in CHARMM TIP3P water. Osmotic concentrations are plotted against total molar concentrations of both electrolytes. The vc MP2–1-1MAM parameter set is based on the Lennard-Jones parameters from the MP2–2MAM “vc MAM” data set while the charges were modified to OPM2: −0.68, PO2: 1.75, OP2: −0.69, yielding a total charge of −1 for MP2 in the simulation. Right: Comparison of osmotic concentrations from simulation using the MP2–1-1MAM (vc MP2–1-1MAM) parameters for methylphosphate, denoted as final parameters, together with different counterions (Na+ and guanidinium) described by the original CHARMM36m parameters.

10. Comparison of NMR relaxation times T 1 and T 2 between experiments and simulations for Ac–K+-pS0-NH2. SP0 is the residue name of phosphoserine0. For mapping of the atom names, consult Figure . The * in the plot marks relaxation times that could not be measured experimentally. The error bars denote standard error of the mean estimated over 3 independent dipeptides in MD simulations and standard errors from the peak fit of the experimental resonances according to eq . The experiment and simulation were conducted at 298 K.

11. Comparison of NMR relaxation times T 1 and T 2 between experiments and simulations for Ac–K+-pS–-NH2. SP1 denotes the residue name of phosphoserine–1. For mapping of the atom names, refer Figure . The * in the plot marks relaxation times that could not be measured experimentally. The error bars denote standard error of the mean estimated over 3 independent dipeptides in MD simulations and standard errors from the peak fit of the experimental resonances according to eq . The experiment and simulation were conducted at 298 K.

12. RMSD of Cα atoms of the anti-sigma factor antagonist SpoIIAA in its unphosphorylated form (1H4Y wt) and with a phosphorylated serine (pS57)­(1H4X), containing either a singly charged phosphate group (P1) or a doubly charged phosphate group (P2). ‘Original’ denotes the original CHARMM36m parameters, ‘new’ refers to our final parameters for MP1 and MP2 listed in Table . The curves have been smoothed using SciPy’s gaussian_filter1d with sigma = 50. As reference structure, the crystal structure of 1H4Y was used.