Post-processing and analysis

Here we assume that you checked the convergence of your system and that all the run ended successfully (no protein unfolding, no other issues). Here we will not focus on convergence, error calculation and metadynamics technical details, considering that such topics are already covered by the Metadynamics tutorial.

Reweighting

After the end of the simulation, we need to reweight our free energy landscape on apt collective variables. As anticipated in the previous step, we will use the distance from the origin of the reference frame $\rho$ and the number of contacts $c$. We thus prepare step-by-step a reweighting script (you can find it in the GitHub folder, called reweight.dat).

We begin from the reading of the relevant data files (refer to the plumed.dat file if you do not remember what is contained in them):

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
rho: READRead quantities from a colvar file. This action has hidden defaults. More details FILEthe name of the file from which to read these quantities=results/coord_rho.dat VALUESthe values to read from the file=rho IGNORE_FORCES use this flag if the forces added by any bias can be safely ignored IGNORE_TIME ignore the time in the colvar file
c: READRead quantities from a colvar file. This action has hidden defaults. More details FILEthe name of the file from which to read these quantities=results/coord_rho.dat VALUESthe values to read from the file=c IGNORE_FORCES use this flag if the forces added by any bias can be safely ignored IGNORE_TIME ignore the time in the colvar file
metad: READRead quantities from a colvar file. This action has hidden defaults. More details FILEthe name of the file from which to read these quantities=results/metad_data.dat VALUESthe values to read from the file=metad.* IGNORE_FORCES use this flag if the forces added by any bias can be safely ignored IGNORE_TIME ignore the time in the colvar file
restr_rmsd: READRead quantities from a colvar file. This action has hidden defaults. More details FILEthe name of the file from which to read these quantities=results/rmsd_restraint.dat VALUESthe values to read from the file=restr_rmsd.* IGNORE_FORCES use this flag if the forces added by any bias can be safely ignored IGNORE_TIME ignore the time in the colvar file

After loading the data, we have to perform the reweighting of the metadynamics potential via the Tiwary-Parrinello estimator. Concurrently, we remove the (almost negligible) contribution of the RMSD restraining, obtaining the final weights for the histogram

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# --- Click here to reveal hidden parts of input file ---- 
weights: REWEIGHT_BIASCalculate weights for ensemble averages that negate the effect the bias has on the region of phase space explored More details TEMPthe system temperature=300 ARG the biases that must be taken into account when reweighting=metad.rbias,restr_rmsd.bias

Having the weights, we can compute the histogram, considering the data from 200 ns to 1 µs to ignore the out-of-equilibrium portion of the simulation:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# --- Click here to reveal hidden parts of input file ---- 
HISTOGRAMAccumulate the average probability density along a few CVs from a trajectory. More details ...
  ARGthe quantities that are being used to construct the histogram=rho,c
  GRID_MIN the lower bounds for the grid=0.,0
  GRID_MAX the upper bounds for the grid=3.0,160
  GRID_BINthe number of bins for the grid=300,160
  KERNEL the kernel function you are using=DISCRETE
  LOGWEIGHTSthe logarithm of the quantity to use as the weights when calculating averages=weights
  LABELa label for the action so that its output can be referenced in the input to other actions=histo
  UPDATE_FROMOnly update this action from this time=200000
  UPDATE_UNTILOnly update this action until this time=1000000
... HISTOGRAM

This histogram can be converted to a free energy landscape

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# --- Click here to reveal hidden parts of input file ---- 
ff: CONVERT_TO_FESConvert a histogram to a free energy surface. More details GRIDYou should use ARG instead of this keyword which was used in older versions of PLUMED and is provided for back compatibility only=histo TEMPthe temperature at which you are operating=300

and finally printed

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# --- Click here to reveal hidden parts of input file ---- 
DUMPGRIDOutput the function on the grid to a file with the PLUMED grid format. More details GRIDYou should use ARG instead of this keyword which was used in older versions of PLUMED and is provided for back compatibility only=ff FILE the file on which to write the grid=reweighted_fes.dat

From this data we can plot the 2D free energy landscape reweighted on $\rho$ and $c$:

Alt text
Free energy surface projected along the new CVs.

∆G calculation

In the GitHub folder you can find a simple python script to get the free energy difference between two basins, deltaG.py. You have to input the path of the reweighted FES, the minima and the maxima in x and y for both the basins and it returns the free energy difference. For my run, we defined the two minima and obtained this result:

$> python3 deltaG.py reweighted_fes.dat 0.45 0.55 105 115 2. 2.8 0 0.5
The free energy difference between basin A and basin B is -3.53 kcal/mol

We now have $\Delta G_{\text{MetaD}}$, the last step is the calculation of the entropic correction.

Entropic correction calculation

We have to compute the volume of the protein included in the sphere. To do so, we can again use VMD. After loading the structure, we open the Tk console and select all the atom that are within 2.8 nm from the center of mass we computed in the previous step of the tutorial:

set sel [atomselect top "sqrt((x-35.28876876831055)^2 + (y-34.06196594238281)^2 + (z-32.041622161865234)^2) < 28.0"]

having set up this selection, we can write an output PDB file.

$sel writepdb atoms_in_sphere.pdb

Having this file, we can use gmx sasa to have an estimate of the volume of these atoms:

gmx sasa -f atoms_in_sphere.pdb -s atoms_in_sphere.pdb -tv volume.xvg

obtaining a volume of 27.294 nm$^3$.

Putting this value and $\rho_{\text{s}}=2.8$ nm in the formula we get the entropic correction

\[\begin{align*} RT \log\left( \frac{V^{0}}{\frac{4}{3}\pi \rho_{\text{s}}^3 -V_{\text{host}}}\right) &= -2.18 \text{ kcal/mol} \end{align*}\]

And we have a final binding free energy value $\Delta G^{0}=-(3.53+2.18)$ kcal/mol $= -5.71$ kcal/mol, which is close to the experimental value of $-5.2$ kcal/mol obtained by Morton et al. in this work.

Back to VMetad home