The context: therapeutic mRNA design
The COVID-19 mRNA vaccines (BNT162b2 from Pfizer, mRNA-1273 from Moderna) demonstrated that synthetic mRNA is a viable therapeutic platform. But designing an optimal mRNA sequence is not trivial: the same protein can be encoded by trillions of different mRNA sequences (thanks to the degeneracy of the genetic code), and each has different biophysical properties.
The state of the art is LinearDesign (Zhang lab, Nature 2023), which simultaneously optimizes secondary structure (MFE, minimum free energy) and codon usage (CAI). It is an elegant algorithm. But it is missing a dimension.
The missing dimension: nearest-neighbor stacking
When two nucleotides sit next to each other in the mRNA, they interact through stacking forces (van der Waals, hydrophobic interactions). The energy of that interaction depends on the identity of both nucleotides — not just the one you are looking at. This is nearest-neighbor thermodynamics (SantaLucia 1998, Turner 2004).
The stacking ΔG profile sweeps along the mRNA position by position and computes that energy for each pair of adjacent nucleotides. It is local context information that metrics like GC3 (which only looks at the third codon position) cannot capture.
At nucleotide resolution, stacking and base pairing share only 12 % of their variance, and 78 % of the ΔG profile is explained by neither local composition nor secondary structure. Nearest-neighbor stacking captures information that secondary structure does not contain.
Key results
Complementarity, measured at the right scale: averaged over windows the stacking ΔG is mostly composition, so the comparison says little there. At nucleotide resolution, across 13 sequences, stacking shares only 12 % of its variance with base-pairing probability (pooled r −0.345, 95 % CI [−0.378, −0.312]). Optimizing one does not automatically optimize the other.
Expression correlation: ΔG correlates with CSC (Codon Stability Coefficient) at the codon level (r = −0.497) and much more strongly at the intra-gene window level (r = −0.797). The thermodynamic signal predicts mRNA stability.
Universality: In 9 out of 9 genes analyzed, the AT → native → GC optimization maintains the pattern ΔG(AT) > ΔG(native) > ΔG(GC), with a mean effect of -0.416 kcal/mol.
Residual contribution: After controlling for GC3 (which already explains much), ΔG still contributes r = −0.359 (p = 7.5e-27), with ΔR² = +0.044, and the effect is significant in 6 of the 9 genes individually. That ~4 % of explained variance is the fine-tuning when GC3 is already saturated — the fine adjustment that makes the difference between a good vaccine and an excellent one.
The dual optimizer
Two optimization modes implemented as proof-of-concept:
Greedy with lookahead: At each position, evaluates all possible synonymous codons and selects the one that minimizes ΔG locally, with a lookahead window to avoid local optima. Time: 0.1 seconds for the 1,273 codons of the spike protein.
Beam search with integrated MFE: Maintains multiple candidates in parallel (beam_width=5) and every 10 codons evaluates secondary structure with RNAfold. Integrates MFE into the beam pruning decision. Time: 13 seconds for the same sequence.
| Sequence | Mean ΔG | MFE/nt | GC% |
|---|---|---|---|
| Control: S2P_AT_backtrans | -1.181 | -0.117 | 30.1% |
| BNT162b2 (Pfizer) | -1.481 | -0.268 | 56.8% |
| mRNA-1273 (Moderna) | -1.542 | -0.313 | 62.0% |
| Greedy (ΔG only) | -1.593 | -0.345 | 63.9% |
| Beam search (ΔG+MFE) | -1.592 | -0.370 | 63.9% |
The first row is not the native viral sequence: it is an AT-maximal back-translation of the S-2P antigen (the same one used in the vaccines), serving as the optimizer's lower control rather than a biological reference. The real wild-type CDS encodes a different protein and is therefore left out of the comparison.
Both modes reach more favourable ΔG and MFE values than the two commercial vaccines. These are computational measures: none of these sequences has been expressed or measured in the lab, so there is no evidence that the numerical advantage translates into more expression or more stability. Beam search significantly improves MFE compared to greedy (-0.370 vs -0.345) with virtually identical ΔG — demonstrating that integrating secondary structure into the decision adds real value.
Positioning vs LinearDesign
LinearDesign optimizes MFE + CAI through dynamic programming. It is the undisputed state of the art. But it operates in a two-dimensional space.
QMetrika's optimizer adds a third dimension: nearest-neighbor stacking. Given that 78 % of the ΔG variance is explained by neither composition nor base pairing, this dimension is mostly its own — it is not redundant with what LinearDesign already optimizes. The proposal is not to replace LinearDesign, but to complement it with thermodynamic information it is not using.
The difference between a good vaccine and an excellent one may lie in a dimension that current optimizers do not use. Nearest-neighbor stacking is that dimension.
Intellectual property
The optimizer code is not public.
Next steps
Experimental validation (Oct 2026 - Mar 2027): 3-5 constructs in HEK293 — native, greedy-ΔG, beam-ΔG+MFE, GC-max control — measuring protein expression (luciferase reporter) and mRNA stability (half-life). The goal is to correlate predicted ΔG with measured expression and stability.
Pareto frontier: Sweep beam_width (3, 5, 10, 20, 50) and map the ΔG vs MFE trade-off. Is there an optimal point? Or is the frontier smooth and gene-dependent?