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.
## ℹ Did you mean to use `annotate()`?
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.44 1.87 1.00 3855 2019
## ec50_target1jps 1.37 0.13 1.09 1.59 1.00 4683 3670
## ec50_target1mlc 2.39 0.53 0.91 3.18 1.00 3115 2020
## ec50_target1ztx 0.75 0.08 0.59 0.89 1.00 5633 4243
## ec50_target2aep 1.13 0.29 0.63 1.59 1.00 3537 1060
## ec50_target2jel 1.65 0.06 1.53 1.76 1.00 6256 4982
## hill_target1ahw 1.68 0.44 0.89 2.60 1.00 3088 1821
## hill_target1jps 1.50 0.36 0.89 2.31 1.00 3901 3482
## hill_target1mlc 1.04 0.56 0.25 2.32 1.00 2889 3080
## hill_target1ztx 2.71 0.55 1.78 3.91 1.00 6778 5448
## hill_target2aep 1.98 0.68 0.55 3.36 1.00 1535 476
## hill_target2jel 3.19 0.59 2.09 4.43 1.00 7529 5076
## top_target1ahw -10.55 0.75 -11.56 -8.81 1.00 3866 2206
## top_target1jps -9.69 0.59 -10.62 -8.33 1.00 5050 3681
## top_target1mlc -1.23 5.88 -9.86 12.05 1.00 5941 4949
## top_target1ztx -17.44 0.28 -17.98 -16.90 1.00 8988 5820
## top_target2aep -16.20 1.25 -16.96 -14.77 1.01 1477 477
## top_target2jel -11.09 0.34 -11.73 -10.43 1.00 9507 5865
## bottom_target1ahw -26.02 2.26 -31.59 -23.09 1.00 3274 1812
## bottom_target1jps -30.51 3.77 -39.49 -24.92 1.00 4264 3586
## bottom_target1mlc -18.54 4.61 -32.33 -14.78 1.00 2644 1743
## bottom_target1ztx -38.89 3.49 -46.84 -33.10 1.00 5294 4182
## bottom_target2aep -30.28 5.92 -43.99 -21.21 1.00 5628 5171
## bottom_target2jel -19.47 0.98 -21.53 -17.74 1.00 6849 4766
##
## 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 10025 5933
##
## 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.4 1.2
## looic 19042.6 78.0
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.1, 1.4]).
##
## 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