Algorithm for optimized mRNA design improves stability and immunogenicity

Figures

• 

  Fig. 1. Overview of mRNA coding region design for both stability and codon optimality using SARS-CoV-2 spike protein as an example. a, Due to codon degeneracy and combinatorial explosion, there are around 2.4 × 10632 possible mRNA sequences encoding the spike protein. Enumerating every possible sequence would take around 10616 billion years. The pink and blue paths represent the wild-type and the optimally stable (lowest free energy) sequences, respectively. nt, nucleotides. b, The secondary structures of wild-type (left) and optimally stable (right) spike mRNAs. The wild-type mRNA is mostly single-stranded and thus prone to degradation in loop regions (red), whereas the optimally stable mRNA is mostly double-stranded. Optimization using LinearDesign takes around 11 min. c, The application of DFA and lattice parsing in computational linguistics (left) and its adaptation to mRNA design (right). An mRNA DFA (analogous to a word lattice) compactly encodes all mRNA candidates, which are folded simultaneously by lattice parsing to find the optimal mRNA (Fig. 2). d, Two-dimensional visualization of the mRNA design space, with stability (represented by MFE) on the x axis and codon optimality…

• 

  Fig. 2. Illustration of the LinearDesign algorithm. a, Codon DFAs. b, An mRNA DFA (bottom) and lattice parsing on that DFA (top). In the DFA, the optimal mRNA sequence under a simplified energy model is shown as the blue path, together with its optimal structure shown in the dot-bracket format (dots indicate unpaired, and brackets indicate base pairs). In lattice parsing, the brown and black arcs also depict base pairs (two GC pairs and two AU pairs), and the trapezoidal shaded areas depict the decomposition of the optimal structure. Among all mRNA sequences encoded in the DFA, lattice parsing finds the optimal sequence with its optimal structure, achieving the lowest free energy under this energy model, where GC and AU pairs have energies of −3 and −2 kcal mol−1, respectively (Extended Data Fig. 1e). Note that here we use the simplified energy model for illustration, but our implementation uses the nearest-neighbour energy model. c, Another illustration of the optimal sequence and secondary structure in b. d, Joint optimization between stability and codon optimality by integrating the codon optimality in weighted DFAs. Top, the codon frequencies for threonine and serine. The relative…

• 

  Fig. 3. Computational characteristics of the LinearDesign algorithm. a, Runtime analysis of mRNA design for proteins in UniProt (Supplementary Table 1). Overall, our exact search scales quadratically with sequence length (Supplementary Figs. 7 and 8), and our MFE + CAI mode (with λ = 4) is about 15% slower than the MFE-only version. Moreover, beam search (b = 500) significantly speeds up the design of long sequences, with minor search errors (Supplementary Fig. 9). b,c, Two-dimensional (MFE–CAI) visualizations of designs for the SARS-CoV-2 spike (b) and VZV gE (c) proteins, respectively (both using human codon preference). The blue curves form the feasibility limit (optimal boundary), by varying λ from 0 to ∞ (see Extended Data Fig. 4 for λ of (–∞, 0]). The GC percentage is shown in parentheses. The human genome favours GC-rich codons; therefore, codon optimization (pink arrows) also improves stability, but only marginally, as the two optimization directions (codon versus stability) are largely orthogonal. By contrast, with an AU-rich codon preference (such as in yeast), codon optimization decreases stability (Extended Data Fig. 4b). d, Secondary structures of the mRNA designs for…

• 

  Fig. 4. Experimental evaluation of LinearDesign-generated mRNA sequences encoding SARS-CoV-2 spike protein. a, Summary of chemical stability of and protein expression from spike mRNA designs A–G and the corresponding immune response (induction of anti-spike IgG) in mice compared to the codon-optimized benchmark H. The vaccines of mRNA-1273 and BNT-162b2 are annotated with daggers, because they use modified nucleotides, but their MFEs here are calculated with the standard energy model. b, Non-denaturing agarose gel characterization of mRNA, showing the correlation of gel mobility with minimum free energy. For gel source data, see Supplementary Fig. 1a. c, Chemical stability of mRNAs upon incubation in 10 mM Mg2+ buffer at 37 °C. Data are from three independent experiments. Seq., sequence. d, Protein expression levels from mRNAs 48 h after transfection into HEK293 cells, as determined by flow cytometry. Mean fluorescence intensity (MFI) values are derived from three independent experiments. Kruskal–Wallis ANOVA with Dunn’s multiple comparisons with the H group. e–g, C57BL/6 mice (n = 6) were immunized intramuscularly with two doses of formulated mRNA with a two-week interval. e, End-point…

• 

  Fig. 5. Experimental evaluation of LinearDesign-generated mRNAs encoding VZV gE protein. a, Summary of chemical stability of and protein expression from VZV gE mRNA designs and the corresponding immune response (induction of anti-gE IgG) in mice. The ‘sweet spot’ region is highlighted with light blue shading. b, Non-denaturing agarose gel characterization of mRNA showing the correlation of gel mobility with minimum free energy; for gel source data, see Supplementary Fig. 1b. c, Chemical stability of mRNAs upon incubation in 10 mM Mg2+ buffer at 37 °C. Data are from three independent experiments. d, Protein expression levels from mRNAs 48 h after transfection into HEK293 cells, as determined by flow cytometry. MFI values are derived from three independent experiments. Kruskal–Wallis ANOVA with Dunn’s multiple comparisons with the gE-Ther group. e, C57BL/6 mice (n = 5) were immunized intramuscularly with two doses of formulated mRNA with a two-week interval. End-point titre of anti-gE IgG is shown. Two-tailed Mann–Whitney U test. Data are mean ± s.d. (c,d) or geometric mean ± geometric s.d. (e). See Extended Data Fig. 8 for extra experimental results, Supplementary Fig. 11 for predicted…

• 

  Extended Data Fig. 1. Illustrations of the optimization problems in mRNA design, DFA representations, single sequence folding as natural language parsing, and lattice parsing. a–b, Visualization of mRNA design as optimization problems for stability (objective 1, in a) and joint stability and codon optimality (objectives 1 and 2, in b). c–h show how lattice parsing solves the first optimization problem (see Fig. 2d for the second). c, Codon DFAs. d, An mRNA DFA made of three codon DFAs. The thick paths depict the optimal mRNA sequences under the simplified energy model in e, AUGCU⋆UGA, where ⋆ could be any nucleotide. e, Stochastic context-free grammar (SCFG) for a simplified folding free energy model. Each rule has a cost (i.e., energy term, the lower the better), and the dotted arcs represent base pairs in RNA secondary structure. f, Single-sequence folding is equivalent to context-free parsing with an SCFG; the parse tree represents the best secondary structure for the input mRNA sequence. g, We extend single-sequence parsing (top) to lattice parsing (bottom) by replacing the input string with a DFA, where each string index becomes a DFA state, and a span becomes a path between two states. h,…

• 

  Extended Data Fig. 2. Word lattice and lattice parsing in natural language processing, and correspondence between linguistics and biology. a, An example of word lattice (sentence DFA) for speech recognition. b, Simplified language grammar. c, Single sentence parsing with between-word indices, which is a special case of word lattice parsing. d, Illustration of word lattice parsing for speech recognition with given word lattice and language grammar; the dashed blue arcs above the DFA depict the best parsing structure for the optimal sentence “I like this meal”, while the dashed light-blue arcs below the DFA represent the best parsing structure for a non-optimal sentence “alike this veal”. e, Correspondence between computational linguistics (left) and computational biology (right). See also Fig. 1.

• 

  Extended Data Fig. 3. Examples of the DFA representations for extended codons, modified nucleotides, and coding constraints. a, Alternative genetic codes of serine, tryptophan, and threonine. b, Avoiding certain codon. On the left it shows the original DFA of serine (up), in which the red dashed arrows indicating UCA and UCG are chosen to be avoided, resulting in a new DFA (down). On the right it shows removing the rare amber STOP codon (UAG). c, Avoiding a specific adjacent codon pair. d, Extended serine DFA can include chemically modified nucleotides pseudouridine (Ψ), 6-Methyladenosine (m6A) and 5-methylcytosine (m5C).

• 

  Extended Data Fig. 4. Two dimensional (MFE-CAI) visualizations of mRNA designs for the Spike protein using human codon preference (a) and yeast codon preference (b) with positive and negative λ’s. GC% are shown in parentheses. The human genome prefers GC-rich codons that lead to higher CAI designs are with higher GC%, while the yeast genome prefers AU-rich codons that exhibit an opposite relationship between CAI and GC%. See also Fig. 3 for more in silico results of LinearDesign.

• 

  Extended Data Fig. 5. Extra experimental results of LinearDesign-generated mRNAs encoding the Spike protein. a, In-solution stability of sequences A–H in PBS buffer containing 20 mM Mg2+ at 37 °C over the course of 24 h. The degradation experiments were performed in triplicate independently and the data were presented as mean ± s.d. and fitted with a one-phase decay curve. b, Protein expression of mRNAs following transfection into HEK293 cells for 24 h was determined by flow cytometry. MFI values derived from three independent experiments are shown. Kruskal–Wallis analysis of variance (ANOVA) with Dunn’s multiple comparisons test to H group was performed for statistical analysis. ***P < 0.001. See Fig. 4c, d for similar experiments but with 10 mM Mg2+ and 48 h, respectively.

• 

  Extended Data Fig. 6. In-solution stability and protein expression of sequences A, C, H and BNT, for a head-to-head in vitro comparison between LinearDesign and BioNTech-Pfizer mRNA sequences. a-b, In-solution stability of mRNAs in PBS buffer containing 20 mM Mg2+ or 10 mM Mg2+ at 37 °C. Data are from three independent experiments and were presented as mean ± s.d. and fitted with one-phase decay curve. c-d, Protein expression of mRNAs was determined 24 h or 48 h following transfection into HEK293 cells. MFI value is presented as mean ± s.d. Each group has three independent assays and 10,000 live cells were collected for analysis in each assay. Kruskal–Wallis analysis of variance (ANOVA) with Dunn’s multiple comparisons test to BNT group was performed for statistical analysis. **P < 0.01, ***P < 0.001.

• 

  Extended Data Fig. 7. Antibody (Ab) responses induced by sequences A, C, H and BNT-based mRNA vaccines, for a head-to-head in vivo comparison between LinearDesign and BioNTech-Pfizer mRNA sequences. C57BL/6 mice (n = 5) were immunized i.m. with two doses of mRNA vaccines at a 2-week interval. Seven days after boost immunization, levels of anti-Spike IgG (a) and neutralizing Abs (b) against pseudotyped SARS-CoV-2 were measured. Data were presented as geometric mean ± geometric s.d. A two-tailed Mann-Whitney U test was used for statistical analysis. *P < 0.05. See Source Data for details. Source Data