To defeat the current COVID-19 pandemic, which has already claimed 100,000+ deaths as of early April, a messenger RNA (mRNA) vaccine has emerged as a promising approach thanks to its rapid and scalable production and non-infectious and non-integrating properties. However, designing an mRNA sequence to achieve high stability and protein production remains a challenging problem. Recently, it is discovered that greater secondary structure folding stability and optimal codon usage synergize to increase protein expression. 1 The design problem can therefore be formulated as finding the mRNA sequence(s) that are good in both secondary structure stability and codon optimality among the exponentially many synonymous sequences that encode the same protein.
Each amino acid is translated by a codon, which is 3 adjacent mRNA nucleotides. For example, the start codon AUG translates into methionine, the first amino acid in any protein sequence. But due to redundancies in the genetic code (4 3 = 64 triplet codons for 21 amino acids), most amino acids can be translated from multiple codons. This fact makes the search space of mRNA design increase exponentially with protein length, e.g., for the spike protein of SARS-CoV-2 (the virus that causes the COVID-19 pandemic), which contains 1,273 amino acids (plus the stop codon which is part of the mRNA but not part of a protein), there are about 10 632 mRNA candidates. The mRNA design problem, therefore, aims to exploit the redundancies in the genetic code to find more stable and productive mRNA sequences than the wild type in nature.
Our key idea is to show that this design problem can be reduced to a classical notion in formal language theory and computational linguistics, namely the intersection between a Stochastic Context Free Grammar (SCFG) and a Deterministic Finite Automaton (DFA), which dates back to 1961. 2 Here the SCFG represents the folding free energy model and the DFA represents the set of all possible synonymous mRNA sequences that code a given protein. While the use of SCFG in RNA folding is well-known, 3 the use of DFA to encode the mRNA search space and solving the design problem via the intersection of SCFG and DFA are novel. This intersection can be done in O(n 3 ) time, where n is the mRNA sequence length, but it could still be too slow for large n (e.g., n = 3 × (1273 + 1) = 3, 822 nucleotides for the spike protein of SARS-CoV-2), so we further developed a linear-time approximate version, LinearDesign, inspired by our recent work, LinearFold. 4 We also develop two algorithms for incorporating the codon optimality into the design, one based on k-best parsing to find alternative sequences and one directly incorporating codon optimality into the dynamic programming.
After the completion of our algorithm design and implementation of our Python prototype and C++ code, we became aware of two earlier, independent, papers that tackled the same problem of "most stable RNA design" using dynamic programming. 5, 6 The second paper "CDSfold" published in 2016 did not cite the first one, published in 2003. Our work is different in three aspects. First, we reduced the mRNA design problem to "CFG-DFA intersection", a classical problem in formal language theory and computational linguistics, which is more general and can be applied to other scenarios with alternative inputs, whereas the previous algorithms were ad hoc. Second, we integrated codon optimality in the optimization. Third, we further developed a linear-time approximate version with greatly reduced runtime for long sequences and small sacrifices in search quality.
The mRNA design problem can be formulated as follows: given a protein sequence p = p1 . . . pm where each pi is an amino acid, we search, among all possible mRNA sequences r that translate into protein p, the best mRNA sequence r (p), defined as the sequence that has the structure with minimum folding free energy change: r (p) = argmin r∈RNA(p) MFE(r) [1] MFE(r) = min s∈structures(r) ∆G • (r, s) [2] where RNA(p) = {r | protein(r) = p} is the search space, structures(r) is the set of all possible secondary structure for RNA The thick arrows indicate the best mRNA sequence after intersecting this DFA with the context-free grammar (see Fig. 3 ).
sequence r, and ∆G • (r, s) is the free energy change of structure s for RNA r according to an energy model. Note that the mRNA sequence length is n = 3(m + 1) due to the final stop codon, which is not translated into the protein sequence. Next we first show how to represent the exponentially large search space RNA(p) compactly via DFAs, and then discuss how to do this argmin search (over the product of two exponentially large spaces) efficiently via dynamic programming, which can be reduced to the CFG intersection with DFA.
A. DFA representation for amino acid codons and mRNA search space. In the fields of formal language theory and computational linguistics, the DFA graph is typically used to encode ambiguities in languages. 7 We notice that the ambiguity of codon choice for amino acid is similar as to the language ambiguity problem, and can be represented with a DFA graph too. We first illustrate how to represent the amino acid codons using DFAs.
Informally, a DFA is a directed graph with labeled edges and distinct start and end nodes. For our purpose each edge is labeled by a nucleotide, so that each start-to-end path represents a triplet codon. Formally, a DFA is a 5-tuple Q, Σ, δ, q0, F where Q is the set of nodes, Σ is the alphabet (here Σ = {A, C, G, U}), q0 is the start node, F is the set of end nodes (in this work the end node is unique, i.e., |F | = 1), and δ is the transition function that takes a node q and a symbol a ∈ Σ and returns the next node q , i.e., δ(q, a) = q encodes a labeled edge q a → q . Fig. 1 illustrates how DFAs represent 4 different types of amino acids and their codons. All these DFAs have (0, 0) and (3, 0) as their start and end nodes, respectively. Fig. 1A is the DFA representation for methionine, which has only one path. Fig. 1B is the DFA for amino acid valine, whose codon has a choice at the third nucleotide (most amino acids are of this type). Fig. 1C represents the most complex case, serine, leucine, and arginine, which have 6 codons each, and the branching happens at the start node. Fig. 1D is the DFA for the stop codon. It is special because branching happens at the second step, at node (1, 0).
After building DFAs for each amino acid, we can concatenate them into a single DFA D(p) for a protein sequence p, which represents all possible mRNA sequences that translate into that protein
by stitching the end node of each DFA with the start node of the next. The new end node of the protein DFA is (3(m + 1), 0). Fig. 2 gives a DFA of the protein sequence "methionine leucine".
We also define out_edges(q) to be set of outgoing edges from node q, and in_edges(q) to be the set of incoming edges: 
A stochastic context-free grammar (SCFG) is a context-free grammar in which each rule is augmented with a weight. More formally, an SCFG is a 4-tuple N, Σ, P, S 0,0 S 9,0 where N is the set of non-terminals, Σ is the set of terminals (identical to the alphabet in the DFA, in this case {A, C, G, U}), P is the set of weight-associated context-free writing rules, and S ∈ N is the start symbol. Each rule in P has the form A w → (N ∪ Σ) * where A ∈ N is a non-terminal that can be rewritten according to this rule into a sequence of non-terminals and terminals ( * means repeating zero or more times) and w ∈ R is the weight associated with this rule.
SCFGs are used to represent RNA folding. The weight of a derivation (parse tree, or a secondary structure in this case) is the sum of weights of the productions used in that derivation. 3 For example, for the very simple Nussinov-Jacobson model 8 (which simplifies the energy model to the number of base pairs), we can define this SCFG G:
Here the first line is the bifurcating case, the second line is the base-pairing case (with weight −1, and the negative score mirrors the free energy minimization problem), and the third line is the unpaired cases (note S 0 → N N N makes sure the minimum hairpin length is 3, i.e., no sharp turn is allowed).
The standard RNA secondary structure prediction problem under a Nussinov model can be cast a parsing problem: given the above SCFG G and an input RNA sequence, what is the minimum-weight derivation in G that can generate sequence? For example, for input CCAAAGG, the best derivation is showed in Fig. 4 .
The mRNA design problem is now a simple extension of the above single-sequence folding problem to the case of multiple inputs: instead of finding the minimum energy structure (minimum weight derivation) for a given sequence, we find the minimum energy structure (and its corresponding sequence) among all possible structures for all possible sequences. This can be solved by intersecting the SCFG G on the protein DFA D, which results in a bigger SCFG and find the best derivation in G . 2, 9 In the intersected grammar G , each nonterminal has the form q1Aq2, where A ∈ N is an original nonterminal in G and q1 and q2 are two nodes in D; and the new start symbol is q0Sqn where S is the original start symbol in G and q0 and qn are the start and end nodes in D. The bifurcation rule S 0 → S S will become q1Sq3 0 → q1Sq2 q2Sq3 for all (q1, q2, q3) node triplets in D. We can see that this intersection construction generalizes the widely-used CKY algorithm 10, 11, 12 where the triple of states (q1, q2, q3) generalizes the triple of string indices (i, k, j). The CKY algorithm is a special case of intersection when the DFA only encodes one string, e.g., when the protein is made of amino acids that have only one codon (methionion and tryptophan). Similarly, in computational linguistics, this intersection construction is widely used for word-lattice parsing in speech recognition where the word-lattice is a DFA to account for ambiguity in word identity. 13 The terminal rule N 0 → A will become
This intersected grammar G will have N |Q| 2 nonterminals and P |Q| 3 rules in the worst case (|Q| is the number of nodes in D). This resembles the space and time complexities of the CKY algorithm for i = 0 . . . (n − l) do 13: j ← i + l 14: for each qi ∈ nodes(D, i) do 15: for each qj ∈ nodes(D, j) do 16: if j − i > 4 then pairing (no sharp turn) 17: for each (qi+1, nuci) ∈ out_edges(D, qi) do 18: for each (qj−1, nucj−1) ∈ in_edges(D, qj) do 19: if match(nuci, nucj−1) then 20:
UPDATE(qi, qj, score, (nuci, qi+1, qj−1, nucj−1)) 22: for k = (i + 1) . . . (j − 1) do bifurcation 23: for each q k ∈ nodes(D, k) do (where |Q| = n). Indeed, intersection is a generalization of CKY from fixed input (RNA folding) to lattice input (mRNA design). Now we just need to find the best (minimum weight) derivation in G , from which we can read off the best mRNA sequence and its corresponding structure. Fig. 3 shows one of the best derivations for the DFA in Fig. 2 .
C. Bottom-up dynamic programming on Nussinov model. We describe how to do dynamic programming based on CFG intersection with DFA. First, we use bottom-up dynamic programming on the Nussinov-Jacobson energy model as an introduction. Fig. 5 gives the pseudocode for this simplified version. We first build up the given protein's DFA graph, and initialize two hash tables, best to store the best score for each state, and back to store the best backpointer for each state. The base cases (singleton rule) are best[qi, qi+1] ← 0 and back[qi, qi+1] ← nuci for each state (qi, qi+1), where nuci is the edge between qi and qi+1. Next, for each state (qi, qj) it goes through the pairing rule and bifurcation rules, and updates if a better score is found. After filling out the hash tables bottom-up, we can backtrace the best mRNA sequence stored with the backpointers. See The algorithm based on bottom-up dynamic programming runs in cubic time, however, it is still slow for long sequences. Inspired by our previous work, LinearFold, 4 we further developed a linear-time 1: function LINEARDESIGN(p, b) b is beam size 11: for j = 1 . . . n do 12: for each qj−1 ∈ nodes(D, j − 1) do 13: for each qi such that [qi, qj−1] ∈ best do 14: for each (qj, nucj−1) ∈ out_edges(D, qj−1) do 15: backpointer ← (qj−1, nucj−1) 16 :
if j − (i − 1) > 4 then 18: for each (qi−1, nuci−1) ∈ in_edges(D, qi) do 19 : approximation algorithm for mRNA design. We apply beam pruning, 14 a classical pruning technique, to significantly narrow down the search space without sacrificing too much search quality. Fig. 6 gives the pseudocode of simplified LinearDesign algorithm for the Nussinov model, based on left-to-right dynamic programming and beam pruning. LinearDesign replaces bottom-up dynamic programming with a left-to-right parsing. At each step j (the jth position of mRNA sequence), we only keep the top b states with the lowest free energy and prune out the less promissing states, since they are unlikely to be the optimal sequence. Here b, the beam size, is a user-adjustable parameter to balance runtime and search quality. This approximation leads to a significant runtime reduce from O(n 3 ) to O(nb 2 ). Notice that we use b = 100 as default in LinearFold, 4 but in LinearDesign we usually use a larger b because the search space is much larger.
E. Implementation on Turner model. Our real system uses a leftto-right dynamic programming with beam pruning on the Turner nearest neighbor free energy model. 15, 16 We implement the thermodynamic parameters following Vienna RNAfold, 17 except for the dangling ends and special hairpins. Dangling ends refer to stabilizing interactions for multiloops and external loops, which require knowledge of the nucleotide sequence outside of the state (qi, qj). Though it could be integrated in LinearDesign, the implementation becomes more involved. Special hairpins are hairpin loop sequences of length 3, 4, or 6 unpaired nucleotides with folding free energies stored in lookup tables, rather than estimated from features like other sequences. Special hairpins can be also integrated in LinearDesign with some preprocessing. We will include both dangling end and special hairpin stabilities in future versions.
Since CAI is also important for mRNA functional half life, 1 
where m is the protein length, logwi(r) is the measurement of deviation from the optimal codon (0 is optimal) for the ith amino acid (given an mRNA candidate), and λ is a hyperparameter that balances MFE and CAI.
We integrate this equation into LinearDesign dynamic process, i.e., each DFA graph's edge will have a cost so that the combined cost of traversing a local path (choosing a codon) equals logwi. Each edge cost is the "best" of the paths (i.e., codons) that uses this edge. with their CAI, each has a corresponding path in the DFA graph (see Fig. 7B ). For example, codon UCU has a CAI of 0.18, while the best codon AGC has a CAI of 0.24. The edge costs in "UCU path" are 0, 0, log(0.18/0.24) = −0.29, therefore, logwi of UCU is -0.29. The best codon AGC has a logwi of 0, meaning that choosing AGC would not have a cost.
Considering LinearDesign is doing left-to-right dynamic programming with beam pruning, and at each step j states with lower scores are pruned, it is better to incur edge costs as early as possible in a path, which ensures the states with better CAI paths are more likely to survive in each step. Fig. 7C rearranges the edge cost to fit better for LinearDesign. Note that different edge costs in Fig. 7B and C would not affect exact search.
G. The top k best mRNA candidates. An alternative solution for finding mRNA candidates with both stable secondary structure and high CAI is to provide the top k mRNA candidates with the lowest MFE, and post-process them by features such as CAI. Although this is not as principled as the algorithm described in subsection 2F, this solution is easier to implement, and is more flexible in the sense that users are free to add other customized filters.
Inspired by the k-best parsing algorithm, 18 we developed an efficient algorithm to find suboptimal candidates in a sorted order. During the dynamic programming process (forward-phase), instead of just saving the single best prestate as the backpointer for each state, we store alternative prestates that all transit to this state. Then in the backtrace process (backforward-phase), starting from the last state, we query the second best. The answer is one of the two cases: (1) the second best is from another prestate S1; or (2) the second best is from the same prestate S0. In the former case we can find the second best by backtracing the best path going through prestate S1, while in the latter case we keep querying for the second best from the prestate S0. Recursively, we would compute and get as many solutions as needed.
To our knowledge, our algorithm is the first one to output suboptimal results in the mRNA design problem. Two previous algorithms explored searching suboptimal secondary structure for RNA folding problem. 19, 20, 21 Zuker's algorithm is to find diverse suboptimal secondary structures, and Wuchty's algorithm is to find all secondary structures in a given free energy gap. Our algorithm is different from these two in the sense that: (1) ours is for mRNA design problem; and (2) ours can output all top k best candidates in a sorted order.
Combining the k-best algorithm and linearization, LinearDesign is able to quickly design a large set of mRNA candidates, which provide a set of alternative designs for follow up with wet lab experiments.
Some studies have shown that protein translation level will drop if the 5'-end leader region has more secondary structure. 22, 23, 24, 25, 1 Considering this practical issue, LinearDesign can be used to design an ORF with an absence of base pairing in the 5'-end leader region by utilizing a simple strategy.
Instead of designing the most stable sequence for the whole coding region, we leave the 5'-end leader region (e.g., the first 15 nt) unchanged from the wildtype, since the wildtype usually has less structure in this region. Then we use LinearDesign for the rest of the coding region. Because the designed region will be composed of strong base pairs (generally maximizing GC content), it is unlikely for a global structure change when concatenating with the wildtype 5'-end head region, which we observe is often depleted of GC content. Refolding using secondary structure prediction tools for the concatenated sequence, we get its corresponding secondary structure and observe that the first 15 nucleotides are unpaired.
If a wildtype sequence is not available, we can alternatively enumerate all possible sequences in the 5'-end leader region. Because each amino acid has 3 codons on average and the start codon is fixed, the enumeration space of the first 15 nt in the 5'-end region is small (3 4 = 81), which makes the enumeration feasible.
A. Efficiency and scalability. To estimate the run-time complexity of LinearDesign, we use 100 protein sequences from Uniprot 26 following CDSfold, with length from 78 to 2,828 nt (not including the stop codon). We found that there are three sequences whose lengths reported in CDSfold do not match the ones currently in Uniprot, so we removed these sequences, resulting in a dataset with 97 diverse protein sequences.
We compared LinearDesign in exact (infinite beam size) and approximate modes (beam size b = 1, 000 and b = 100). Because CDSfold code is not currently available, we directly use the runtime reported in CDSfold paper as a comparison allowing us to compare the computational complexity. Note that CDSfold results and our results were run in different machines. We run LinearDesign on a machine with 2 Intel Xeon E5-2660 v3 CPUs (2.60 GHz), while CDSfold is run on the Chimera cluster system at AIST, which was reported in the paper as 176 Intel Xeon E5550 CPUs (2.53 GHz). B. Search quality of linear-time approximation. Since the Lin-earDesign algorithm is significantly faster than its exact counterpart, we envision that the LinearDesign algorithm will be used for long sequences. To ensure the quality of the approximation used in Lin-earDesign (with beam pruning mode), we compared the energy gap between the exact (b = ∞) and approximate algorithm (b = 1, 000).
For this analysis, we used the same dataset as in subsection 3A. Fig. 9 represents the folding free energy differences between the exact search and approximation. Fig. 9A compares the free energy changes of the mRNA sequences designed with exact search and with the b = 1, 000 approximation. The x-axis is the free energy of exact search, while the y-axis corresponds to the free energy of approximation. We see that all plots are on or close to the diagonal, which confirms that the folding free energy differences are 0 or small. Fig. 9B shows the trend of free energy differences increases linearly with mRNA sequence length. They y-axis is the percentage of free energy change gap, which is the free energy change gap (∆∆G • ) divided by the total free energy change of the MFE structure (∆G • ). The percentage of free energy change gap is small for all sequences in the dataset. For b = 1, 000, all sequences have gaps within 1%. Even for the longest sequence (8,484 nt), the gap is 0.8%. For b = 100, most gaps are within 5%, and the largest gap is 7%. We also investigate the percentage of free energy change gap against beam size for two specific protein sequences, the spike protein of SARS-CoV-2 and EGFP (GenBank KM042177.1), in Fig. 9C . The purple curve shows the result of the spike protein. Starting with a small beam size, b = 20, the gap is 10.6%. With increasing beam size, the gap shrinks quickly to less than 6% at b = 100. Further increasing b up to 500, the gap drops to 1%. With a beam size of b = 2, 000 the gap drops to 0, which indicates that the approximate result is the same as the exact search. The EGFP result (the green curve) has the same shape as the spike protein, but the gap drops faster and down to less than 2% at b = 100. This is because EGFP is shorter (239 amino acids), thus the approximation is close to the exact result even with smaller beam size. At b = 200, the gap increases 0.25% compared with b = 100. This happens because more states are kept when enlarging the beam size, among which states with higher scores at length j survive. Their extension to longer lengths (like offspring) beat states and fill the beams, but all these states become worse at the end, resulting a small drop of the search quality. The fluctuation happens in this sequence, but the jump is small and the gap quickly decreases to 0 at b = 400 and above.
C. Example for the coding mRNA of the spike protein. The spike protein of SARS-CoV-2, which has 1,273 amino acid residues, is the target of mRNA vaccines (https://clinicaltrials. gov/ct2/show/NCT04283461).Therefore, we use the spike protein of SARS-CoV-2 as an example, and compare our designed sequences with the wildtype sequence and random generated sequences.
We use the mRNA sequence from the reference genome of SARS-CoV-2 (https://www.ncbi.nlm.nih.gov/nuccore/ MN908947) as the wildtype sequence, which contains 3,822 nucleotides (including the stop codon). Additionally, we use two different strategies to generate random sequences (5,000 sequences for each strategy) as another baseline. One of the two strategies, named "pure random", is to randomly (with equal probabilities) choose a codon for each amino acid in the spike protein, and form a mRNA sequence by concatenation. The other strategy, called "codon-biased random", is to choose a codon the probability proportional to its usage frequency. We also run LinearDesign (in both exact mode and approximation mode) to evaluate the best sequences we can achieve. Since previous studies show that both the folding free energy change of mRNA secondary structure and codon optimality influence the vaccine effectiveness, 1 we do a two-dimensional comparison, MFE and codon adaptation index (CAI), 27 between mRNA sequences. Fig. 10 shows the results. The wildtype sequence, denoted in a red circle, folds into a structure with the minimum free energy change of -967.80 kcal/mol, and has a low CAI of 0.655. Most "pure random" sequences, denoted in blue cloud, have similar free energy changes (-987.90 kcal/mol on average) as the wildtype, but with higher CAI (0.671 on average). This may be because SARS-CoV-2 just recently infected human cells and does not have enough mutations to optimize for human codons. Compared with the wildtype and "pure random" sequences, "codon-biased random" sequences, denoted in A UCUU 1,000 Fig. 11 . The secondary structures of the wildtype and designed mRNA sequences that translate into the spike protein of SARS-CoV-2. A: Wildtype sequence and its secondary structure. B: The best mRNA sequence (with lowest free energy change) designed by LinearDesign with b = 1, 000. C: The best mRNA sequence designed by LinearDesign with b = +∞, i.e. an exact search. D: A zoom in of C showing the 5'-end, which is base paired with the 3' end. E: The best mRNA sequence designed by LinearDesign with b = +∞, while using the wildtype sequence for the first 15 nt nucleotide at 5'-end. These 15 nucleotides do not base pair with the designed sequence and therefore remain unstructured. F: A zoom in of E showing just the 5'-end and the 3' end.
orange cloud, have much higher CAI (0.768 on average), and slight improvement on MFE (-1063.23 kcal/mol on average). We also notice that both "pure random" and "codon-biased random" sequences are packed in cloud-shaped small regions. This is because the search space of possible mRNAs is extremely huge and most of the random sequences have similar MFE and CAI. On the left (with much lower MFE) we plot our designed sequences. The blue plots are optimized by MFE only. The leftmost one is the sequence designed in exact search mode, which has the lowest MFE of -2,477.70 kcal/mol and a CAI of 0.726. The MFE gap between our best designed sequence and the wildtype, as well as random sequences, is large (more than 1,300 kcal/mol). With only 0.56% MFE loss from the exact search sequence, the designed sequence with beam size b = 1, 000 achieves an MFE of -2,463.8 kcal/mol and a higher CAI of 0.751. Compared to the exact mode, which takes 1.6 hours for designing the sequence, the approximation with b = 1, 000 only takes 16 minites, resulting in a 6× speed-up. For b = 100 and b = 20, the MFE are still lower than -2,200 kcal/mol, with CAI both at around 0.735 and 0.725, respectively. Our designed sequences, for both exact search sequence and approximate search sequences, are much better than random ones and the wildtype in terms of MFE.
We show that the top 10, 000, 000 suboptimal sequences for exact mode (the dark-red cloud on the right of exact design sequence). The MFE of the sequences are very close to the optimal one, and the free energy gaps are within 20 kcal/mol. Some of the sequences have higher CAI, e.g., some have the CAI higher than 0.730. This shows that our k best algorithm can be used to select sequences with low MFE and relative higher CAI as vaccine candidates.
Further, we show the results of MFE and CAI joint optimization.The light-blue curve is joint optimization design using exact mode. Each point on the curve is with a different λ, which balances the importance of MFE and CAI. We see that the curve is on the topleft of the figure, indicating that the sequences on the curve have both stable secondary structures and higher CAI. In fact, this curve is the accessible boundary of all possible sequences, i.e., no sequences can achieve the region beyond (to the top-left) the curve. The points on the curve are good candidates for an mRNA vaccine. For example, the point with λ = 100, has the free energy change of -2,414.6 kcal/mol and CAI of 0.823, which is only 2.5% away from the optimal MFE sequence but with 0.097 increase in CAI. We observe that with λ > 100, the sequences on the curve have better CAI than codon-biased random sequences. Shifting right from the light-blue curve with a small margin, the magenta curve is the results of joint optimization using b = 1, 000. This curve shows that the approximation quality is good with b = 1, 000. We also designed a CAI-greedy sequence, which greedily chooses the best codon for each amino acid, leading to a special sequence with CAI=1. We see that the two curves both point to the CAI-greedy design, comfirming that the designed sequences achieve better CAI but sacrifice MFE with increasing λ, and reach CAI-greedy design with a large λ (e.g., λ = 3, 000).
Note that the MFE of the wildtype, the CAI-greedy design and random sequences are calculated by Vienna RNAfold 17 with "-d0" (disable stabilizing interactions for multiloops and external loops), to make fair comparisons with our designed sequences. All the sequences can be refolded using RNAfold without "-d0", by which the points will shift to the left. Figure 11 shows the secondary structures of the wildtype sequences, our designed sequences with b = 1, 000 and b = +∞, as well as designed sequences with less structures at the 5'-end leader region. We can see that the secondary structures of the wildtype (Fig. 11A) have a large number of loops, and our designed sequences ( Fig. 11B and C) have longer helices and fewer loops, which makes the structure more stable. The designed sequence with b = 1, 000 (Fig. 11B) has similar free energy changes as the one with b = +∞ (Fig. 11C) , but it has a multiloop in the middle.
Additionally, we investigate the effectiveness of our strategy for less structure design at the 5'-end leader region. Fig. 11E is the whole secondary structure of the designed sequence with the goal of leaving the 5' end unpaired (the first 15 nucleotides kept identical to the wildtype and the remaining nucleotides designed by LinearDesign), and we zoom in the 5'-end in Fig. 11F . As a comparison, we also zoom in the 5'-end of designed sequence without constraint in Fig. 11D . This demonstrates our strategy can keep the 5' end unstructured, whereas designing the complete sequence results in basepairing at the 5' end.
The mRNA design problem is of utmost importance, especially for mRNA vaccines during the current COVID-19 pandemic. We reduced this problem into a classical problem in formal language theory and computational linguistics, namely the intersection of a CFG (encoding the energy model) with a DFA (encoding the mRNA search space). This reduction provides a natural O(n 3 )-time CKY-style bottom-up algorithm, where n is the mRNA sequence length, but this algorithm might still be too slow for long proteins such as the spike protein of SARS-CoV-2, a promising candidate for an mRNA vaccine. Inspired by our recent work of LinearFold, we then developed a left-to-right algorithm, LinearDesign, which employs beam search to reduce the runtime to O(n), with the cost of exact search. LinearDesign is orders of magnitude faster than exact search (i.e., b = +∞) and suffers only a small loss in folding free energy. For example, for this spike protein, LinearDesign can finish in 16 minutes while exact search takes 1.6 hours, and the free energy difference is only 0.6%. We also developed two algorithms for incorporating codon optimality (CAI) into the consideration, one using k-best algorithms to compute suboptimal sequences and one directly integrating CAI into dynamic programming. Our work provides efficient computational tools to speed up and improve mRNA vaccine development.