Direct RNA sequence design under codon constraints using expressive tensor-based secondary structure models

Abstract: Nucleic acid sequence design via codon optimization is a fundamental task with applications across synthetic biology, mRNA therapeutics, and vaccine design. Given a target protein, it is a major open challenge to navigate the combinatorially large design space of codon sequences mapping to its amino acid sequence. Computational approaches generally seek to optimize simple objectives based on the codon sequence, possibly together with more complicated contributions based on secondary structure analysis. In this work, we demonstrate a direct and efficient algorithm to sample sequences from a suitable Boltzmann distribution defined in terms of the codon sequence and a fully detailed secondary structure free energy model, as well as related algorithms for exact computation of statistical quantities such as free energies, base pairing probabilities, and base and codon marginals. These algorithms draw upon a recently developed tensor-based formulation of secondary structure thermodynamics and demonstrate, for the first time, that global sequence design can be accomplished with respect to a highly accurate free energy model. Moreover, the algorithms can leverage any available CPU and GPU resources in parallel for massive computational speedups.

• The paper introduces exact O(n³) dynamic-programming algorithms that combine codon-constrained sequence design with expressive tensor-based RNA secondary-structure energy models.
• The framework supports Boltzmann sampling, base-pair and codon marginals, and minimum-energy designs while incorporating objectives such as CAI and CPB through tensor-train representations.
• Experiments on proteins up to 2,411 amino acids show practical performance, including GPU speedups of up to 100×, while demonstrating that simplified energy models can produce substantially poorer designs under full-model evaluation.

Overview

This paper addresses the codon optimization problem: given a target amino acid sequence $\psi$, design an mRNA sequence $\phi$ that codes for $\psi$ while optimizing secondary structure stability and, optionally, sequence-based objectives such as the Codon Adaptation Index (CAI) or Codon Pair Bias (CPB). The central contribution is a family of exact, $O(n^3)$ dynamic programming algorithms that sample RNA sequences from a Boltzmann distribution defined jointly over sequences, unpseudoknotted secondary structures, and hidden tensor configurations, using a fully detailed tensor-based free energy model. The authors state that this is the first demonstration that global sequence design can be performed with respect to a highly accurate free energy model via dynamic programming, and this claim appears well supported by their positioning against prior work, which relied either on simplified free energy models (e.g., CDSfold [Terai2015], LinearDesign [Zhang2023-algorithm], DERNA [Gu2024DERNA]) or on heuristic local search with unbounded convergence.

Problem formulation

The authors define an ensemble $D(\psi)$ over triples $(\phi, s, \zeta)$ of coding RNA sequences $\phi \in \Phi(\psi)$, compatible unpseudoknotted secondary structures $s$, and hidden configurations $\zeta$ underlying each structure in the tensor model, with probability mass proportional to $e^{-q(\phi,s,\zeta)}$. Marginalizing over these degrees of freedom yields physically meaningful quantities: the Gibbs free energy of a structure $\phi$0, of a sequence $\phi$1, and of the entire coding ensemble $\phi$2. Soft sequence-based constraints are incorporated through a biasing function $\phi$3 added to the energy, defining a modified distribution $\phi$4 without altering the algorithmic structure. Beyond sampling, the framework supports exact computation of equilibrium base pair probabilities averaged over all coding sequences, per-position base and codon marginal probabilities, and minimization of $\phi$5 via semiring substitution.

Tensor train encoding of codon constraints

The key representational device is a tensor train formulation of the hard and soft codon constraints. Each position $\phi$6 carries a 3-index tensor $\phi$7 whose contraction trace defines the sequence-level energy $\phi$8. Two concrete constructions are given: a codon pair tensor train, in which transfer indices track neighboring codon identities and support arbitrary pairwise codon bonuses (e.g., CPB), and a leaner single codon tensor train supporting per-codon bonuses (e.g., CAI). Both exactly implement the hard constraint that $\phi$9 codes for $\psi$0, so that the unconstrained partition function counts $\psi$1. Theoretical cost prefactors relative to fixed-sequence analysis are modest — approximately $\psi$2 for the single codon formulation and $\psi$3 for the codon pair formulation — and lower in practice.

Combining sequence and structure dynamic programs

The methodological core is the observation that the codon-constrained sequence summation can be fused into the tensor-network secondary structure partition function computation by replacing the folding function $\psi$4 of prior work [Fornace2025-new] with a generalized function $\psi$5 that additionally sums over base identities at each incorporated unpaired base or base pair, propagating hard/soft constraints through the tensor train's $\psi$6 indices. Crucially, to the outer dynamic program the combined sequence-structure degrees of freedom can simply be flattened, so that each block $\psi$7 becomes a rectangular $\psi$8 matrix and all downstream algorithms — recursions, GPU-accelerated triangular matrix multiplication, backtracking — carry over unchanged. This is what makes the approach practical: it inherits the efficiency of the tensor-based analysis machinery rather than requiring new algorithmic development.

Downstream algorithms

Using the evaluation algebra framework of Fornace et al., the paper derives several exact algorithms at no additional asymptotic cost:

• Sampling: direct i.i.d. Boltzmann sampling of joint $\psi$9 states via priority-queue backtracking, with a simultaneous sampling scheme amortizing shared recursion elements across many samples.
• Base pair probabilities: an exact formula for the sequence-averaged pairing matrix $O(n^3)$0, computed from forward matrices $O(n^3)$1 and wraparound "reverse" matrices $O(n^3)$2.
• Marginals: per-position base and codon marginal probabilities in $O(n^3)$3 each after the initial $O(n^3)$4 computation, enabling exact expected codon frequencies and iterative reweighting to match target host codon frequencies (converging in 8–15 iterations with Anderson acceleration).
• Minimization: MFE sequence-structure pairs via tropical (MinSum) semiring substitution, yielding the zero-temperature limit of the ensemble.

An important caveat is that the framework solves $O(n^3)$5 exactly but only approximately solves $O(n^3)$6; the authors argue the sampling problem is more fundamental in any case, both practically and theoretically.

Computational studies

The methods are demonstrated on four proteins of practical interest: NY-ESO-1 (179 aa), ECD_HUMAN (644 aa), SARS-CoV-2 spike (1273 aa), and lipomycin polyketide synthase (2411 aa). Using the $O(n^3)$7 tensor model regressed against the rna06 nearest-neighbor model, sampled ensembles from $O(n^3)$8 are systematically more stable than uniform codon sampling, and surrogate-minimized sequences achieve extremely low free energies. A notable negative result: designs produced under a simplified $O(n^3)$9 model — essentially flat per-base-pair energies — evaluate substantially worse under the full model, leading the authors to conclude that simplified free energy models do not suffice for real-world stability optimization. This directly challenges the adequacy of earlier dynamic-programming design methods built on coarse models.

Performance results are strong. Design costs only about $D(\psi)$0 the cost of analyzing a single fixed sequence (measured prefactors 5.3–6.8 across systems), with GPU acceleration providing up to ~100× speedups; multi-GPU parallelism is demonstrated for the first time in this line of work. A 5000-amino-acid design takes approximately 7 minutes. For the largest system (2409 aa), computing $D(\psi)$1 takes 74 s on 2 GPUs, sampling 2500 sequences adds only 7.7 s, and MFE computation (CPU-only, 64 cores) takes 2567 s — consistent with the authors' observation that partition function computations are better suited to hardware acceleration than min-sum recursions.

Limitations and open questions

Several limitations are acknowledged explicitly. The framework excludes pseudoknotted structures, which may be significant for large natural complexes, though approximate extensions are anticipated. It handles single-stranded coding regions only, omits start/stop codons and noncoding regions (the latter trivially includable), and relies on a free energy model whose parameters were fit to reproduce the rna06 nearest-neighbor ensemble to roughly 10% accuracy rather than to experimental melt data directly — a conservative choice the authors note could be improved without changing the algorithms. Model renormalization (equalizing spectral radii of the $D(\psi)$2 tensors) is required to avoid biasing the sequence sampling prior, and this normalization is heuristic. Some common soft constraints, such as Relative Codon Bias and KL-divergence-based Relative Codon Pair Bias scores, do not fit naturally into the tensor train formulation, although the authors argue their marginal-based frequency matching obviates most such needs. Motif-based constraints are compatible in principle but increase tensor train dimensionality in ways not yet benchmarked. Finally, conditional structural objectives based on equilibrium base pair probabilities remain future work, as does CUDA acceleration of the MFE path.

Conclusion

This paper establishes that codon optimization under a fully detailed, expressive secondary structure free energy model can be performed exactly and efficiently via dynamic programming, by encoding codon constraints as tensor trains and fusing them into tensor-network structure recursions. The resulting framework yields principled Boltzmann sampling, exact statistical marginals, and competitive runtimes with hardware acceleration, offering an ensemble-based alternative to both simplified-model DP methods and heuristic local search. Its extensibility to data-driven soft constraints makes it a plausible foundation for integrating learned codon fitness models with rigorous thermodynamic optimization.