A step-by-step explainer · Rosetta ddg_monomer

Inside a ΔΔG prediction

Bench's DDG tool takes a structure and a point mutation and hands back one number. This page opens up the four things that happen in between — repacking, the ramped repulsion, fifty repeats averaged three at a time, and the honest meaning of "r = 0.73" — because every one of them is a decision that changes the answer.

Companion to the Levitate Bench DDG documentation. For what ΔG and ΔΔG are, start with ΔG and ΔΔG.

The whole thing in one sentence

A predicted ΔΔG is the difference between two scores of the same protein — and because no single score means anything on its own, the entire method is machinery for making sure both structures were relaxed equally hard before you subtract.

Everything below is one pipeline. Rosetta's ddg_monomer application (driven by ddGMover) takes your PDB and a mutation, produces a relaxed model of each sequence, scores both, and reports the difference in Rosetta Energy Units:

Input
Structure + mutation
One PDB, one substitution such as A26W.
Step 1
Repack side chains
Monte-Carlo search over a rotamer library, in both sequences.
Step 2
Minimize ×3
Gradient minimization at 10%, then 33%, then 100% repulsive weight.
Step 3
50 repeats
Same starting structure, different random seeds; average the best three.
Output
ΔΔG in REU
Positive = destabilizing. Not kcal/mol without a fitted rescaling.

Before anything else — the two energies that are missing

ΔΔG of folding is defined from four free energies: folded and unfolded, wild type and mutant. ddg_monomer computes two. It never builds an unfolded state at all. The entire method rests on the approximation that a single substitution changes the unfolded ensemble's energy by the same amount in both sequences, so that term subtracts away:

\[ \Delta\Delta G \;=\; \underbrace{\big(G^{\text{mut}}_{\text{folded}}-G^{\text{WT}}_{\text{folded}}\big)}_{\text{what Rosetta scores}} \;-\; \underbrace{\big(G^{\text{mut}}_{\text{unfolded}}-G^{\text{WT}}_{\text{unfolded}}\big)}_{\text{assumed}\;\approx\;0} \]

As the Bench docs put it: "the energy of an unfolded protein is nearly identical with a single point mutation, so the dominant factor in DDG is the energy difference of the folded state." That is usually true and occasionally badly false — proline substitutions restrict the unfolded chain, glycine substitutions loosen it, and charge changes alter unfolded-state electrostatics. When a ddg_monomer prediction is spectacularly wrong for a proline or glycine mutation, this cancelled term is the first suspect. The four-corner cycle this comes from is step 4 of the companion page.

A sign-convention collision worth catching now

The Bench documentation writes free energy in the unfolding direction: \(\Delta G_{\text{folding}} = (H_{\text{unf}} - H_{\text{fold}}) - T(S_{\text{unf}} - S_{\text{fold}})\). That is \(G_{\text{unfolded}} - G_{\text{folded}}\), so larger = more stable. To keep "positive ΔΔG = destabilizing" with that sign, the subtraction has to run wild type minus mutant — which is exactly what the docs write: \(\Delta G_{\text{mutation}} = (\Delta G_{\text{folding}})_{\text{WT}} - (\Delta G_{\text{folding}})_{\text{mutant}}\).

The companion page uses the opposite ΔG sign (\(G_{\text{folded}} - G_{\text{unfolded}}\), negative = stable) and therefore the opposite subtraction order, mutant minus wild type. Both land on the same convention for the headline number: positive ΔΔG means the mutation made the protein worse, which is also what ddg_monomer prints. Two flipped signs cancelling is not a coincidence — it is people insisting on that one shared meaning.

Step 1

You cannot just swap the side chain

The naive way to score a mutation: delete the wild-type side chain, graft the new one on in the same orientation, score. This fails immediately, and it fails in a way that always points the same direction.

A protein core is packed at roughly the density of an organic crystal. Drop a tryptophan where an alanine was and its atoms land on top of the neighbours. The repulsive term explodes — tens of Rosetta Energy Units from a handful of atom pairs that are 0.5 Å too close. Nature's answer to that is not "this mutation is impossible"; it is that the surrounding side chains rotate a little and make room. So the score you want is not the score of the naive graft, but the score after the neighbourhood has re-accommodated.

That re-accommodation is repacking. Side chains are not allowed to bend arbitrarily: Rosetta uses a discrete rotamer library, the handful of torsion angles each amino acid actually adopts in known structures, each carrying a small penalty for how rare it is. Choosing one rotamer per residue to minimize the total score is a combinatorial problem — with 7 residues and 9 rotamers each there are \(9^7 \approx 4.8\) million states — so Rosetta searches it with Monte Carlo simulated annealing rather than enumeration.

The widget runs that actual search. Seven residues sit around a small pocket; each has nine rotamers; the score is a real pairwise sum of a repulsive term (quadratic in overlap) plus a contact bonus plus the rotamer penalties. Mutate the amber site from Ala to Trp and watch the clash appear — then repack, and watch three quarters of it disappear.

A seven-residue pocket, repacked by simulated annealing

live search · press the stages in order

Red bonds are steric overlaps — pairs of atoms closer than their contact distance. A well-packed core still carries a small repulsive term, so the wild-type baseline is not zero; what matters is how the term changes. The ΔΔG on the right is measured against the packed wild-type score, so it is the same quantity ddg_monomer reports. Note what changes it: not the mutation, but how much of the structure you let move afterwards.

The formal version — the packer, and why annealing rather than enumeration

Fix the backbone. Each residue \(i\) is assigned a rotamer \(r_i\) from its library, and the score decomposes into one-body and two-body terms:

\[ E(\mathbf{r}) \;=\; \sum_i E_1(r_i) \;+\; \sum_{i \(E_1\) collects the rotamer's intrinsic probability penalty and its interaction with the fixed backbone; \(E_2\) is the pairwise atom-atom energy between two side chains. Because the backbone is fixed, every \(E_2(r_i,r_j)\) can be precomputed into a table once, after which evaluating a candidate assignment is table lookups rather than geometry.

Minimizing \(E\) over \(\mathbf{r}\) is NP-hard in general (it is a discrete Markov-random-field labelling problem), so Rosetta anneals: propose a single-residue rotamer substitution, accept it with probability \(\min(1, e^{-\Delta E/T})\), and lower \(T\) geometrically from a high starting value to near zero. High \(T\) early lets the search cross barriers; low \(T\) late makes it settle. The widget uses the same acceptance rule with \(T\) falling from 3.0 to 0.02.

Kellogg et al. found the best-performing protocol repacks with a softened repulsive term and only then minimizes with the hard one — which is the subject of the next step.

Step 2

Why the repulsion is turned down to 10% first

After repacking, ddg_monomer minimizes: it follows the gradient of the score downhill, letting side-chain torsions and small backbone motions relax continuously. The Bench docs describe the distinctive part precisely — it minimizes three times, with the repulsion term weighted at 10%, then 33%, then 100%.

This looks like a fudge and is not. Gradient minimization can only go downhill; it stops at the first local minimum it reaches. In a packed protein the walls between one packing arrangement and a better one are made almost entirely of the repulsive term — two atoms that would have to pass through each other to reach the better arrangement. At full weight, those walls are hundreds of REU tall and utterly impassable. Structures get pinned wherever the mutation dropped them.

Scaling the repulsion to 10% flattens the walls without touching the attractive landscape, so the structure can slide across a barrier it would otherwise never cross. Restoring the weight to 33% and then 100% then re-imposes real sterics from the better side of the wall. It is the same trick as annealing, run on a continuous coordinate: soften, move, re-harden.

Below, one coordinate stands for the whole conformational change. The smooth part of the score (a broad attractive well, minimum on the right) is what the structure wants; two narrow spikes are steric clashes in the way. The teal ball minimizes at 100% weight throughout. The amber ball follows the ramp. Same start, same gradient descent, same final energy function — different answers.

Gradient minimization through a steric wall: hard-only vs. 10 → 33 → 100%

staged · drag the starting point

The drawn curve is the score the minimizer currently feels, so it visibly flattens and re-sharpens as the ramp advances. Both balls are scored on the same 100%-weight function at the end — the ramp is a search strategy, not a different energy.

The formal version — soft-rep, and what the ramp is really approximating

Rosetta's full-atom repulsion fa_rep is a Lennard-Jones-like term that rises steeply inside the pair's contact distance. Write the total score as

\[ E_w(\mathbf{x}) \;=\; E_{\text{attr}}(\mathbf{x}) \;+\; w\,E_{\text{rep}}(\mathbf{x}), \]

with \(w\) the repulsive weight. Every \(E_w\) has the same global structure of attractive basins; what \(w\) controls is the height of the barriers between them. Minimizing at \(w = 0.1\) is therefore a homotopy: you solve an easy, smooth problem first and use its solution to initialize the next, harder one — the same idea as continuation methods in optimization or graduated non-convexity in vision.

The reason it is 10 → 33 → 100 rather than one jump from 10 to 100 is that each intermediate keeps the structure inside the basin found at the previous weight. Jump straight to full repulsion and a structure sitting in a mildly overlapped arrangement gets a huge gradient kick that can throw it back out.

In the Kellogg et al. row-16 protocol the repack uses a soft-repulsive score function and the subsequent quasi-Newton minimization uses the hard one, with harmonic distance restraints tethering the model to the input structure so that "minimization" cannot quietly drift into a different fold.

Step 3

Fifty repeats, and why the answer is the average of the best three

Repacking is stochastic. Minimization inherits that randomness through its starting point. Run the pipeline twice on the same input and you get two different scores. So ddg_monomer runs it 50 times from the same original structure and — this is the part worth dwelling on — reports the average of the best three.

Why not the average of all 50? Because a failed trajectory, one that got stuck in a bad packing, produces a score that is too high, and there is no corresponding failure mode that produces a score that is too low. The distribution has a sharp lower edge near the true minimum and a long upper tail. Averaging everything averages in the failures.

Why not simply the single best of 50, then? Because the minimum of 50 noisy draws is itself a very noisy statistic: it is determined by one trajectory, so it inherits that trajectory's full variance, and it systematically chases whichever run got the luckiest noise. The mean of the best three is the compromise: it discards the failures like the minimum does, but averages three draws instead of trusting one.

And this matters twice over, because ΔΔG is a difference of two such estimates — the noise in the wild-type run and the noise in the mutant run both land in the final number.

The widget runs the whole experiment. Each trajectory is a real noisy descent on a rugged landscape, and the true answer is known (it is the difference of the two landscapes' global minima, 0.954), so bias and spread can be measured rather than asserted.

Three estimators of the same ΔΔG, measured against the truth

staged · runs real stochastic trajectories

Stage 3 repeats the entire 100-trajectory experiment 30 times and plots where each estimator's ΔΔG lands relative to the truth. Root-mean-square error is the honest summary: it charges an estimator for bias and for scatter at once.

The formal version — bias, variance, and why three

Let \(E_1,\dots,E_{50}\) be the trajectory scores, independent draws from a distribution that piles up just above the true minimum \(E^\star\) and trails off upward. Three estimators of \(E^\star\):

\[ \hat E_{\text{mean}} = \tfrac{1}{50}\textstyle\sum_i E_i, \qquad \hat E_{\min} = \min_i E_i, \qquad \hat E_{3} = \tfrac{1}{3}\textstyle\sum_{k=1}^{3} E_{(k)} \]

where \(E_{(k)}\) is the \(k\)-th order statistic. \(\hat E_{\text{mean}}\) is biased upward by the full mean of the failure tail and that bias does not shrink with more samples. \(\hat E_{\min}\) is asymptotically unbiased for \(E^\star\) but has variance set by a single extreme draw; extreme-value statistics converge slowly, so its spread stays stubbornly wide. \(\hat E_3\) averages three order statistics: the variance falls roughly like the average of three correlated extremes while the bias picks up only the small gap between the first and third best.

For ΔΔG the relevant quantity is the difference \(\hat E^{\text{mut}} - \hat E^{\text{WT}}\). Biases partially cancel when the two landscapes are similarly rugged — but a mutation that makes the local packing problem harder raises the mutant's failure rate, so the cancellation is imperfect, and that residue is a genuine systematic error in the method, not sampling noise you can average away.

Step 4

What "REU" and "r = 0.73" actually entitle you to say

The number that comes out is in Rosetta Energy Units. The Bench docs are careful here, and the care is the point: REUs "are not directly convertible to kcal/mol", though "a conversion factor has been determined that can be used for approximation."

The reason is structural, not a missing calibration. A Rosetta score is a weighted sum of empirical terms — van der Waals attraction and repulsion, implicit solvation, hydrogen bonding, rotamer probability, a reference energy per amino-acid type — whose weights were fitted, largely to reproduce native structures and rankings. Nothing in that procedure forces the sum onto a calorimetric scale. Any REU→kcal/mol factor is a regression slope fitted on a benchmark set, and it is only valid for the score function, protocol, and kind of mutation it was fitted on. Kellogg et al. got theirs by re-training the reference energies against experimental ΔΔG data, which is what put that particular energy function on a kcal/mol-matched scale in the first place.

Then the accuracy claim. The Bench docs quote r = 0.73 across 1,210 mutations, citing Kellogg et al. 2011. In that paper the best sampling protocol reaches r = 0.69 on those 1,210 mutations, rising to 0.73 once the score function's reference weights are re-trained on ΔΔG data. Either number is genuinely good for a structure-based predictor. Neither means what people usually hear.

A correlation of 0.73 explains \(r^2 \approx 53\%\) of the variance. The remaining half is spread across individual predictions — the exact things you use the tool for. The widget makes that concrete: it draws a synthetic benchmark with the correlation you dial in, then lets you set the threshold at which you would call a mutation destabilizing and counts how often that call is wrong.

What a correlation of r = 0.73 does to individual decisions

synthetic data · drag r and the threshold

The 1,210 points are generated, not the Kellogg benchmark — only the correlation and the set size are borrowed, so that the decision arithmetic below is arithmetic you can trust rather than a real result being re-plotted. The Pearson r printed on the canvas is recomputed from the generated points.

The formal version — from r to the error rate on one prediction

Model predicted and experimental values as jointly normal with correlation \(r\), both standardized. Conditioned on the experimental value \(t\), the prediction is

\[ p \mid t \;\sim\; \mathcal{N}\!\big(r\,t,\; 1-r^2\big), \]

so the residual standard deviation of a single prediction is \(\sqrt{1-r^2}\). At \(r = 0.73\) that is \(0.68\) — you have removed about a third of the uncertainty in any individual value, not three-quarters. In the widget's units, where experimental ΔΔG has a spread of 1.53 kcal/mol, a single prediction carries 1.04 kcal/mol of irreducible scatter — the same size as the effect most single mutations have.

This is why ΔΔG predictors are used the way they are: to rank and triage a library of candidate mutations, where errors on individual entries wash out across the set, rather than to decide the fate of one specific substitution. Ranking is a far weaker demand on the model than valuation, and \(r = 0.73\) supports it comfortably.

Putting it together

Return to the kernel: a ΔΔG prediction is a difference of two scores of the same protein. Every step of ddg_monomer is machinery in service of making that subtraction fair.

Where to go next

The companion page ΔG and ΔΔG covers what these numbers mean thermodynamically — populations, sign conventions, and why the difference of differences is the only trustworthy quantity. Free energy perturbation covers the other computational route to the same number: instead of scoring two relaxed structures, morph one amino acid into the other and integrate the work.

Source for the protocol described here: the Levitate Bench DDG documentation and Kellogg EH, Leaver-Fay A, Baker D, Role of conformational sampling in computing mutation-induced changes in protein structure and stability, Proteins 2011;79(3):830–8.