(Note: These tutorials are meant to provide illustrative examples of how to use the AMBER software suite to carry out simulations that can be run on a simple workstation in a reasonable period of time. They do not necessarily provide the optimal choice of parameters or methods for the particular application area.)
Copyright Vinícius Cruzeiro. 2018


Section 4: Running and analyzing production simulations at different pH and Redox Potential values

We will now see how we can run and analyze simulations at different values of pH and redox potential. This will allow us to make more accurate Eo and pKa predictions, and also to evaluate the pH-effect on Eo and the redox potential-effect on pKa values. For simplicity, we will only make calculations in implicit solvent here, however, the same procedure can be followed in explicit solvent.

Running the simulations

We now need to replicate the input file used in the production phase shown in Section 3 for different pH and redox potential values. This can be easily done by using a simple bash script that makes use of the sed Linux command. The first step is to change the value of solvph by PH and solve by REDOX in a reference input file. PH and REDOX are pointers to be replaced by the script with values of pH and redox potential. This is how the reference input file should look like:

prod.implicit.mdin.ref
Production Stage in Implicit Solvent
&cntrl
  icnstph=1,       ! Switch for constant pH (=1, implicit solvent)
  ntcnstph=5,      ! Protonation state change attempt every ntcnstph steps
  solvph=PH,       ! Solvent pH
  icnste=1,        ! Switch for constant redox potential (=1, implicit solvent)
  ntcnste=5,       ! Protonation state change attempt every ntcnste steps
  solve=REDOX,     ! Solvent redox potential (in Volts)
  ntt=3,           ! Temperature scaling (=3, Langevin dynamics)
  gamma_ln=10.0,   ! Collision frequency of the Langevin dynamics in ps-1
  ntc=2,           ! SHAKE constraints (=2, hydrogen bond lengths constrained)
  ntf=2,           ! Force evaluation (=2, hydrogen bond interactions omitted)
  ntb=0,           ! Boundaries (=0, non-periodic)
  cut=1000.0,      ! Cutoff
  dt=0.002,        ! The time step in picoseconds
  nstlim=10000000, ! Number of MD steps to be performed
  ig=-1,           ! Random seed (=-1, get a number from current date and time)
  ntwr=10000,      ! Restart file written every ntwr steps
  ntwx=10000,      ! Trajectory file written every ntwx steps
  ntpr=10000,      ! The mdout and mdinfo files written every ntpr steps
  ioutfm=1,        ! Trajectory file format (=1, Binary NetCDF)
  igb=2,           ! GB model
  saltcon=0.1,     ! Salt concentration
  irest=1,         ! Flag to restart the simulation
  ntx=5,           ! Initial condition (=5, coord. and veloc. read from the inpcrd file)
/

We show below the simple bash script cphemd_geninputs.sh that helps generating the input files. For simplicity, only 6 values of pH and 6 values of redox potential are being considered. However, the script can be easily modifed to cover more values.

cphemd_geninputs.sh
#!/bin/bash

# Reference input file
infile=prod.implicit.mdin.ref

for PH in {5.0,5.5,6.0,6.5,7.0,7.5}
do
  for REDOX in {-0.263,-0.233,-0.203,-0.173,-0.143,-0.113}
  do
    newfile=${infile%.ref}.pH_$PH.E_$REDOX
    cp $infile $newfile
    sed -i "s/PH/$PH/g" $newfile
    sed -i "s/REDOX/$REDOX/g" $newfile
  done
done

This script will generate 36 input files in the format prod.implicit.mdin.pH_*.E_* (note that the suffix .ref was removed).

We now need to execute the simulations. As the simulations are independent, they can be executed simultaneously in parallel to save time. In case you want to execute the simulations in serial, the following script can be used:

cphemd_run.sh
#!/bin/bash

# Reference input file
infile=prod.implicit.mdin.ref
# Topology file
prmfile=mp8_is.prmtop
# Input rst7 file
crdfile=mp8_is.equi.rst7
# CPIN file
cpin=mp8_is.cpin
# CEIN file
cein=mp8_is.cein
# Prefix for the output files
prefix=mp8_is.prod

for PH in {5.0,5.5,6.0,6.5,7.0,7.5}
do
  for REDOX in {-0.263,-0.233,-0.203,-0.173,-0.143,-0.113}
  do
    newfile=${infile%.ref}.pH_$PH.E_$REDOX
    pmemd.cuda -O -i $newfile -p $prmfile -c $crdfile \
               -x $prefix.nc.pH_$PH.E_$REDOX -inf $prefix.mdinfo.pH_$PH.E_$REDOX \
               -o $prefix.mdout.pH_$PH.E_$REDOX -r $prefix.rst7.pH_$PH.E_$REDOX \
               -cpin $cpin -cpout $prefix.cpout.pH_$PH.E_$REDOX \
               -cprestrt $prefix.cprestrt.pH_$PH.E_$REDOX \
               -cein $cein -ceout $prefix.ceout.pH_$PH.E_$REDOX \
               -cerestrt $prefix.cerestrt.pH_$PH.E_$REDOX
  done
done

The following tar file contains the CPOUT and CEOUT files generated in the simulations: couts.tar.gz.

Analysing the simulations

Theoretical background

From the Hendersen-Hasselbalch equation we can devise an expression for the fraction of protonated species of a given pH-active residue:

f_{prot} = \frac 1 {1 + 10 ^ {n (pH - pK_a)}}       [1]

where n is the Hill coefficient.

Equivalently, from the Nernst equation we can devise an expression for the fraction of reduced species of a given redox-active residue (see reference [1] for more information):

f_{red} = \frac 1 {1 + e ^ {n \frac{\nu F}{k_b T} (E - E^o)}}       [2]

where n is the Hill coefficient, \nu is the number of electrons being reduced in the redox reaction, F is the Faraday constant, k_b is the Boltzmann constant, T is the temperature, and E is the redox potential.

Using fitpkaeo.py

The fitpkaeo.py AmberTool allows us to automatically get fitted values of pKa and Eo using equations 1 and 2. This is obtained by passing information from cphstats or cestats for multiple CPOUT or CEOUT files to fitpkaeo.py through pipe. fitpkaeo.py is able to even generate a graph containing the fractions from the simulations against the fitted curve.

For example, to fit the Eo value of residue HEH at pH 7.0 and generate the graph, we can use the command:

for ceout in mp8_is.prod.ceout.pH_7.0.E_-0.*; do cestats -i mp8_is.cein $ceout; done | fitpkaeo.py --graph --graph-path fracredu.png

This will generate the following output:

  HEH 14: Eo =    -0.202811 (+- 0.000296 ) V; Hill coefficient =  1.028 (+- 0.011 ); Temperature = 300.00 K
Graph generated to  fracredu.png


and the following figure:

As can be seen, we successfully obtain a Eo value very close to our target value of -203 mV at pH 7.0. The values inside the parenthesis represent fitting errors of the properties shown. The low fitting errors are justified in the figure. As we can see, the fractions of reduced species obtained through the simulations for different redox potential values at pH 7.0 lie very well on top of the fitted curve.

Let's now fit the pKa values of all pH-active residues at redox potential -203 mV. We can do this with the following command, which also generates a graph:

for cpout in mp8_is.prod.cpout.pH_*.E_-0.203; do cphstats -i mp8_is.cpin $cpout; done | fitpkaeo.py --graph --graph-path fracprot.png

This command will generate the following output:

  GL4 12: pKa =   4.53952 (+- 0.00894 ) ; Hill coefficient =  1.130 (+- 0.017 )
  PRN 15: pKa =   6.07351 (+- 0.04300 ) ; Hill coefficient =  0.817 (+- 0.063 )
  PRN 16: pKa =   6.16959 (+- 0.01876 ) ; Hill coefficient =  0.921 (+- 0.034 )
Graph generated to  fracprot.png


and the following figure:

pH-dependence of Eo predictions

We can use the following command to generate Eo fittings at different pH values:

for PH in {5.0,5.5,6.0,6.5,7.0,7.5}
do
  echo pH: $PH
  for ceout in mp8_is.prod.ceout.pH_$PH.E_*; do cestats -i mp8_is.cein $ceout; done | fitpkaeo.py
done


This will generate the following output:

pH: 5.0
  HEH 14: Eo =    -0.186541 (+- 0.000713 ) V; Hill coefficient =  0.984 (+- 0.024 ); Temperature = 300.00 K
pH: 5.5
  HEH 14: Eo =    -0.189420 (+- 0.001411 ) V; Hill coefficient =  1.057 (+- 0.055 ); Temperature = 300.00 K
pH: 6.0
  HEH 14: Eo =    -0.198200 (+- 0.000624 ) V; Hill coefficient =  1.020 (+- 0.023 ); Temperature = 300.00 K
pH: 6.5
  HEH 14: Eo =    -0.199265 (+- 0.000638 ) V; Hill coefficient =  1.025 (+- 0.024 ); Temperature = 300.00 K
pH: 7.0
  HEH 14: Eo =    -0.202811 (+- 0.000296 ) V; Hill coefficient =  1.028 (+- 0.011 ); Temperature = 300.00 K
pH: 7.5
  HEH 14: Eo =    -0.205623 (+- 0.000310 ) V; Hill coefficient =  0.963 (+- 0.010 ); Temperature = 300.00 K

We can plot this data to obtain the Eo of HEH as a function of pH.

As discussed in reference [1], Eo should decrease with pH. As the figure shows, this behavior is correctly described by our simulations.

Redox potential-dependence of pKa predictions

The following command prints the fittings of the pKa values of all pH-active residues for different redox potential values:

for REDOX in {-0.263,-0.233,-0.203,-0.173,-0.143,-0.113}
do
  echo Redox Potential: $REDOX V
  for cpout in mp8_is.prod.cpout.pH_*.E_$REDOX; do cphstats -i mp8_is.cpin $cpout; done | fitpkaeo.py
done


This will generate the following output:

Redox Potential: -0.263 V
  GL4 12: pKa =   4.71920 (+- 0.03529 ) ; Hill coefficient =  1.311 (+- 0.116 )
  PRN 15: pKa =   6.18241 (+- 0.06711 ) ; Hill coefficient =  0.766 (+- 0.088 )
  PRN 16: pKa =   6.15350 (+- 0.03770 ) ; Hill coefficient =  0.766 (+- 0.049 )
Redox Potential: -0.233 V
  GL4 12: pKa =   4.76361 (+- 0.03594 ) ; Hill coefficient =  1.363 (+- 0.137 )
  PRN 15: pKa =   6.09123 (+- 0.04907 ) ; Hill coefficient =  0.835 (+- 0.075 )
  PRN 16: pKa =   6.25073 (+- 0.04109 ) ; Hill coefficient =  0.740 (+- 0.051 )
Redox Potential: -0.203 V
  GL4 12: pKa =   4.53952 (+- 0.00894 ) ; Hill coefficient =  1.130 (+- 0.017 )
  PRN 15: pKa =   6.07351 (+- 0.04300 ) ; Hill coefficient =  0.817 (+- 0.063 )
  PRN 16: pKa =   6.16959 (+- 0.01876 ) ; Hill coefficient =  0.921 (+- 0.034 )
Redox Potential: -0.173 V
  GL4 12: pKa =   4.37096 (+- 0.11188 ) ; Hill coefficient =  0.856 (+- 0.117 )
  PRN 15: pKa =   5.96110 (+- 0.03184 ) ; Hill coefficient =  0.835 (+- 0.049 )
  PRN 16: pKa =   6.11422 (+- 0.01552 ) ; Hill coefficient =  0.851 (+- 0.024 )
Redox Potential: -0.143 V
  GL4 12: pKa =   4.18668 (+- 0.01406 ) ; Hill coefficient =  0.900 (+- 0.013 )
  PRN 15: pKa =   5.97212 (+- 0.02218 ) ; Hill coefficient =  0.756 (+- 0.029 )
  PRN 16: pKa =   6.07389 (+- 0.04610 ) ; Hill coefficient =  0.730 (+- 0.056 )
Redox Potential: -0.113 V
  GL4 12: pKa =   3.94276 (+- 0.26421 ) ; Hill coefficient =  0.616 (+- 0.121 )
  PRN 15: pKa =   6.00046 (+- 0.05143 ) ; Hill coefficient =  0.793 (+- 0.072 )
  PRN 16: pKa =   6.13903 (+- 0.05881 ) ; Hill coefficient =  0.704 (+- 0.067 )

By plotting this data we obtain the pKa of each pH-active residue as a function of redox potential.

We see the pKa values decrease with increases in the redox potential, in agreement with the expectation (see reference [1]).

How can we increase convergence?

From the last two figures, it is clear that our results are still not converged. We would need to run longer to get a better convergence, or, as discussed in reference [1], make use of Redox Potential Replica Exchange MD (E-REMD) to obtain improved convergence. Reference [1] shows that E-REMD significantly improves sampling convergence in comparison to performing C(pH,E)MD only (e.g., without exchange attempts, as done in this section). The procedure to execute E-REMD in Amber is equivalent to running pH-REMD. Therefore, a nice pH-REMD tutorial written by Jason Swails that can be found here can be used as a reference. As opposed to the option "-rem 4" using in pH-REMD, in E-REMD we use "-rem 5". Furthermore, in a follow-up publication in reference [2] it has been shown that the use of multidimensional REMD further improves the convergence of the results in comparison to E-REMD. More information about E-REMD and MultiD-REMD can also be found in Amber's manual.

The converged behavior of Eo versus pH and the pKa of each pH-active residue versus redox potential has been reported in reference [1]. There, simulations were performed longer, using E-REMD, and more values of pH and redox potential were considered.


Click here to go back to the Introduction


Click here to go back to Section 3


References

[1] Vinícius Wilian D. Cruzeiro, Marcos S. Amaral, and Adrian E. Roitberg, "Redox Potential Replica Exchange Molecular Dynamics at Constant pH in AMBER: Implementation and Validation", J. Chem. Phys., 149, 072338 (2018). DOI: https://doi.org/10.1063/1.5027379

[2] Vinícius Wilian D. Cruzeiro, and Adrian E. Roitberg, "Multidimensional Replica Exchange Simulations for Efficient Constant pH and Redox Potential Molecular Dynamics", J. Chem. Theory Comput., 15, 871–881 (2019). DOI: https://doi.org/10.1021/acs.jctc.8b00935

(Note: These tutorials are meant to provide illustrative examples of how to use the AMBER software suite to carry out simulations that can be run on a simple workstation in a reasonable period of time. They do not necessarily provide the optimal choice of parameters or methods for the particular application area.)
Copyright Vinícius Cruzeiro. 2018