Apply: Sigmoid Model -- Pnear Folding Funnel
apply_sigmoid_model_Pnear.RmdModeling Folding Funnels
A common task in molecular modeling is to predict the conformation of the folded state for a given molecular system. For example, the Rosetta ab initio, or protein-protein-interface docking protocols. To turn the simulation into a prediction requires predicting the relative free energy of the folded state relative a reference.
The Rosetta score function can score individual conformations, but doesn’t capture the free energy of the state. Typically, a researcher will run a series of trajectories and generate a score vs. RMSD plot and look for a “folding funnel” e.g. lower energies for conformations that are closer to a target folded state. Here, RMSD is the root-mean squared deviation measuring the euclidean distance of pairs of atom defined by the application (for example just the backbone for sequence design or interface atoms for docking).
Pnear score
To quantify the quality of the folding funnel, recently, there has been interest in using the Pnear score, which is defined by
Pnear = Sum_i[exp(-RMSD[i]^2/lambda^2)*exp(-score[i]/k_BT)] /
Sum_i[exp(-score[i]/k_BT)]where (RMSD[i], score[i]) is the score RMSD and score values for a conformation i. The parameter lambda is measured in Angstroms indicating the breadth of the Gaussian used to define “native-like-ness”. The bigger the value, the more permissive the calculation is to structures that deviate from native. Typical values for peptides range from 1.5 to 2.0, and for proteins from 2.0 to perhaps 4.0. And finally the parameter k_BT is measured in in energy units, determines how large an energy gap must be in order for a sequence to be said to favor the native state. The default value, 0.62, should correspond to physiological temperature for ref2015 or any other scorefunction with units of kcal/mol.
Two state model
Thinking of the folded and unfolded states as a two-state model and RSMD as a reaction coordinate or “collective variable”, then the energy gap can be modeled by a sigmoidal Boltzmann distribution.
## Warning in xkcd::theme_xkcd(): Not xkcd fonts installed! See
## vignette("xkcd-intro")## Warning in theme_xkcd(): Not xkcd fonts installed! See vignette("xkcd-intro")
plot of chunk sigmoidal-cartoon
For a principled molecular dynamics or Monte Carlo simulation that maintains detailed balance, it is in theory possible to use thermodynamic integration to quantify the energy gap between the two states. However, this is often not computationally feasible for proteins of moderate size or in a protein design or screening context where many different molecules need to be evaluated given a limited computational budget. So, Instead, we will assume that the at least locally around the folded state, the degrees of freedom increase exponentially so that the log of the RMSD defines a linear reaction coordinate.
If we simulate, trajectory points from the sigmoid on the log(RMSD) scale, with a Normally distributed error we can generate synthetic score-vs-rmsd plots
## Warning in GeomIndicator::geom_indicator(mapping = ggplot2::aes(indicator = paste0("Pnear: ", : All aesthetics have length 1, but the data has 300 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
plot of chunk simulate-score-vs-rmsd-data
A nice thing about having the parametric model to generate score-vs-rmsd plots, is that it allows us to measure measure the sensitivity of the Pnear to differently shaped score-vs-rmsd plots. For example we can scan both the radius of
plot of chunk simulate-data-scan
Another question we can use this model to investigate is how reproducible is the Pnear score?
plot of chunk simulated-replicates
Antibody SnugDock Case study
As a case study, we can look at the real score-vs-rmsd plots from the Antibody SnugDock scientific benchmark. It is evaluates the SnugDock protocol over 6 Antibody protein targets
We can use the fit the sigmoid model to the log(RMSD)
using the BayesPharma package, which relies on
BRMS and Stan
Check the model parameter fit and estimated parameters:
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: response ~ sigmoid(ec50, hill, top, bottom, log_dose)
## ec50 ~ 0 + target
## hill ~ 0 + target
## top ~ 0 + target
## bottom ~ 0 + target
## Data: data (Number of observations: 3003)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## ec50_target1ahw 1.70 0.11 1.45 1.87 1.00 4597 2714
## ec50_target1jps 1.37 0.13 1.07 1.60 1.00 4199 3043
## ec50_target1mlc 2.39 0.52 0.96 3.17 1.00 2741 2016
## ec50_target1ztx 0.75 0.08 0.59 0.89 1.00 5111 3856
## ec50_target2aep 1.11 0.25 0.58 1.51 1.00 4927 3196
## ec50_target2jel 1.65 0.06 1.53 1.76 1.00 5928 5263
## hill_target1ahw 1.69 0.43 0.91 2.62 1.00 3637 2596
## hill_target1jps 1.50 0.37 0.88 2.30 1.00 3547 3000
## hill_target1mlc 1.05 0.55 0.27 2.29 1.00 2972 2934
## hill_target1ztx 2.72 0.55 1.76 3.91 1.00 6168 5216
## hill_target2aep 2.00 0.66 0.75 3.36 1.00 4343 2235
## hill_target2jel 3.19 0.58 2.13 4.41 1.00 7193 5504
## top_target1ahw -10.57 0.67 -11.59 -8.95 1.00 4882 2870
## top_target1jps -9.69 0.59 -10.63 -8.31 1.00 4673 3178
## top_target1mlc -1.42 5.83 -9.88 12.31 1.00 6140 4915
## top_target1ztx -17.44 0.28 -17.99 -16.89 1.00 9799 5483
## top_target2aep -16.33 0.47 -16.94 -15.47 1.00 5453 2790
## top_target2jel -11.09 0.34 -11.73 -10.41 1.00 8342 5971
## bottom_target1ahw -25.95 2.09 -31.02 -23.10 1.00 3745 2320
## bottom_target1jps -30.53 3.80 -39.70 -24.89 1.00 3415 2998
## bottom_target1mlc -18.38 4.30 -31.43 -14.75 1.00 2561 2009
## bottom_target1ztx -38.76 3.46 -46.52 -32.95 1.00 5696 4004
## bottom_target2aep -30.40 6.00 -44.30 -21.31 1.00 7151 4817
## bottom_target2jel -19.45 0.97 -21.52 -17.75 1.00 6116 4923
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 5.74 0.07 5.60 5.89 1.00 9878 5783
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Excitingly, using leave-one-out cross-validation, the sigmoid model fits the data very well
##
## Computed from 8000 by 3003 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -9521.3 39.0
## p_loo 21.5 1.1
## looic 19042.6 78.0
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.6]).
##
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.Visualize the fit as draws from the expected mean and median quantile intvervals on the log(RMSD) scale:
plot of chunk plot-data-model
And on the original RMSD scale:
plot of chunk plot-data-model-rmsd