which forcefield suits the best for divalent metal ions like calcium in protein in gromacs system and why.
the calcium makes bidentate bonds with protein and water.
Dear Akanksha,
Attached please find a paper entitled "Accurate Description of Calcium Solvation in
Concentrated Aqueous Solutions " which describes methods for calculating calcium in protein systems.
The following are computational methods used:
Computational details
Classical molecular dynamics simulations were carried out for aqueous solutions of CaCl2. The molar ratio used here was 4 moles of CaCl2 to 55.55 moles of water, for convenience hereafter referred to as 4 m. Two different water models were employed, namely SPC/E41 and, to test transferability of the developed calcium force field, also TIP4P/2005. 42 The results for TIP4P/2005 water not shown in the main text, can be found in the SI (Table 8-12 and Fig. 6 and 7). The unit cell contained 5750 water molecules, 414 calcium dications, and 818 chloride anions. Initially, the Lennard–Jones parameter from Ref. 28 were used for chloride and the parameters from the Gromos 53a6 force field27 for calcium, see Table 1.
The mixed Lennard–Jones parameters were derived from self-parameters using the Lorentz– Berthelot mixing rules. Additionally, the geometric averages for both sigma and epsilon parameters, as implemented, e. g., in the Optimized-Potentials for Liquid Simulations 43,44 (OPLS) force field were used. The results of the latter calculations are provided in the Supporting Information (SI) (see Table 3-7 and Fig. 2-5). Additionally, in the same tables and figures of the SI, results for a smaller calcium ion (A1), for a chloride ion from literature 45 (A2) and for recently published parameters 31 (A3) can be found. The exact values, including the ones for the counter ions of set A1 and A2, are given in Table 2 of the SI.
All simulations were performed using the Gromacs 4.5.4 program package 46 with a time step of 1 fs employing 3D periodic boundary conditions. Simulations were carried out under constant pressure of 1 bar, controlled by the Parrinello–Rahman barostat 47 with a coupling 4
Table 1: Force field parameters of the ions employed in this study.
Ca
σ (Å) (kJ/mol) charge
full 2.8196 0.5072 +2.0
ECC 2.8196 0.5072 +1.5
ECCR 2.5376 0.5072 +1.5
Cl
σ (Å) (kJ/mol) charge
full 4.4499 0.4184 -1.0
ECC 4.4499 0.4184 -0.75
ECCR 3.7824 0.4184 -0.75
constant of 0.5 ps and constant temperature of 298 K, controlled by the velocity rescaling thermostat 48 with a coupling time of 0.5 ps. Simulations of a length of at least 30 ns were performed for models with both scaled and unscaled charges. Data from the last 10 ns were used for further analysis, e. g., calculations of the radial distribution functions (RDFs), and comparison with neutron scattering data. For the description of the water molecules constraints for bonds and angles were applied using the Lincs algorithm. 49 The electrostatic and van der Waals interactions were truncated at 1.2 nm and the long range electrostatic interactions were accounted for by the particle mesh Ewald method. 50 Since the original ionic sizes were parametrized together with the full charge model, we refined them after charge scaling according to the experimental number density and the neutron scattering data. 19 Decreasing the van der Waals radius of calcium by 10 % and chloride by 15 % resulted in number densities within 5 % of the experimental values, 8,19 see Table 2. The results from this size rescaling are denoted as electronic continuum correction with rescaling (ECCR) throughout the paper. Note, that changing the size of the chloride ion influences the density more significantly than changing the size of the calcium (see SI), since there are as twice as
many chloride ions in the simulation box than calcium ions, moreover, they have a larger van der Waals radius.
Subsequent simulations at constant temperature and volume with a length of 50 ns were performed in order to obtain equilibrated input structures for viscosity calculations and 5 calculations of the residence times of water molecules in the calcium solvent shell. The last snapshot of each of these equilibrium runs was taken as input for non-equilibrium simulations to calculate viscosities. 51,52 For systems employing the SPC/E water model, three acceleration rates of 0.001, 0.0025, and 0.005 nm/ps2 were tested. The corresponding trajectories have a length of at least 20 ns and the last 10 ns were taken for the viscosity calculations.
Based on the results, an acceleration rate of 0.0025 nm/ps2 was determined as an appropriate choice and used also for the simulations with the TIP4P/2005 water model (see Table 10).
Hoping this will be helpful,
Rafik
- Badges
- Science topic