Whole-cell modeling of E. coli confirms that in vitro tRNA aminoacylation measurements are insufficient to support cell growth and predicts a positive feedback mechanism regulating arginine biosynthesis

Abstract In Escherichia coli, inconsistencies between in vitro tRNA aminoacylation measurements and in vivo protein synthesis demands were postulated almost 40 years ago, but have proven difficult to confirm. Whole-cell modeling can test whether a cell behaves in a physiologically correct manner when parameterized with in vitro measurements by providing a holistic representation of cellular processes in vivo. Here, a mechanistic model of tRNA aminoacylation, codon-based polypeptide elongation, and N-terminal methionine cleavage was incorporated into a developing whole-cell model of E. coli. Subsequent analysis confirmed the insufficiency of aminoacyl-tRNA synthetase kinetic measurements for cellular proteome maintenance, and estimated aminoacyl-tRNA synthetase kcats that were on average 7.6-fold higher. Simulating cell growth with perturbed kcats demonstrated the global impact of these in vitro measurements on cellular phenotypes. For example, an insufficient kcat for HisRS caused protein synthesis to be less robust to the natural variability in aminoacyl-tRNA synthetase expression in single cells. More surprisingly, insufficient ArgRS activity led to catastrophic impacts on arginine biosynthesis due to underexpressed N-acetylglutamate synthase, where translation depends on repeated CGG codons. Overall, the expanded E. coli model deepens understanding of how translation operates in an in vivo context.


Model Design
The cycle of tRNA usage -aminoacylation by aminoacyl-tRNA synthetases and tRNA release after amino acid transfer by elongating ribosomes -is a system that consists of both deterministic and stochastic features. While tRNA aminoacylation can be modeled deterministically as mass action kinetics with ordinary differential equations (ODEs), the ribosome's utilization of aminoacyl-tRNAs for synthesizing polypeptides depends on the identity of the codon at the ribosome A site, which necessarily depends on the mRNA sequence and may be more appropriately represented by a stochastic model such as the Gillespie algorithm.
For example, aminoacyl-tRNA utilization by ribosomes can be represented as the following one-step reaction: aminoacyl-tRNA i + ribosome j → tRNA i + ribosome j+1 + H + (1) where aminoacyl-tRNA i is the aminoacylated form of the i-th tRNA isoacceptor and ribosome j is an active ribosome with an open A site at the j-th codon on the mRNA transcript. Notably, i indicates the identity of the tRNA isoacceptor whereas j indicates a position along the mRNA transcript. Consequently, as the ribosome processes along the mRNA, the identity of the codon at the j-th position changes and the tRNA isoacceptors that are eligible to interact with the j-th codon also change.
In addition to the tRNA aminoacylation cycle, the whole-cell modeling framework also consists of qualities that lends itself well to both deterministic and stochastic modeling approaches. On one hand, the simulation time step (which is on the order of seconds) enables the translation process to read 2 to 5 × 10 5 codons per second by 1 to 3 × 10 4 active ribosomes, each processing at approximately 17.5 amino acids per second -a fast reaction rate that can be appropriately modeled deterministically as bulk changes. On the other hand, the positions of active ribosomes on mRNA transcripts are detailed to the codon residing at the A site (previously detailed to the residue of the nascent polypeptide), and the sequence order of codons on mRNAs is obeyed -meaning that ribosomes cannot progress to the next codon until the current one has been successfully read. This codon-by-codon procession of ribosomes along mRNAs, which depends on the identity of the j-th codon and the availability of the i-th aminoacyl-tRNA, can be suitably modeled by a stochastic simulation approach.
Thus, to be compatible with both the deterministic and stochastic qualities of the tRNA aminoacylation cycle and the whole-cell modeling framework, the following three-step strategy was developed for use during simulations: 1. Calculate the kinetic limitations of aminoacyl-tRNA synthetases, assuming a constant codon reading rate, during the time ∆t: (a) Estimate the codon reading rate as a bulk change during ∆t, based on the upcoming codon sequences of mRNAs being translated by active ribosomes and the anticipated ribosome elongation rate.
(b) Simulate the tRNA aminoacylation cycle as a deterministic ODE model during time ∆t, where tRNA aminoacylation and aminoacyl-tRNA utilization by ribosomes (estimated as the codon reading rate in step 1(a)) are modeled by mass action kinetics.
(c) Retrieve the number of times each codon type was read during ∆t from the ODE solution.
2. Process ribosomes along mRNAs according to the kinetic limitations and the sequence order: (a) Process ribosomes codon-by-codon -so as to obey the sequence order of codons on mRNAs -as far as possible without exceeding the number of codons read in step 1(c).
(b) Retrieve the number of times each codon type was read after all ribosomes are processed as far as possible.
3. Reconcile disagreements, if any, between the number of codons read in step 1(c) and step 2(b): (a) Perform a random search for the optimal combination of ribosome positions along mRNAs that minimize disagreements while obeying the kinetic limitation determined in step 1(c).
(b) Any disagreements that remain by this point are kinetic over-estimations (where the ODE model estimated more codon reading events than the sequence order permitted); reconcile these disagreements by undo-ing the appropriate number of codon reading events in the ODE solution (ie. returning an unaminoacylated tRNA to its aminoacylated form), which may require undo-ing the number of aminoacylation events (ie. returning an aminoacyl-tRNA to its unaminoacylated form) if the pool of aminoacyl-tRNAs is too small.
The result of these steps is a representation of tRNA aminoacylation cycles and codon-by-codon ribosome procession that obeys both the kinetic constraints of aminoacyl-tRNA synthetases and the sequence order of codons encountered by ribosomes along mRNAs.

tRNA Aminoacylation by Aminoacyl-tRNA Synthetases
Aminoacylation of each tRNA isoacceptor was represented as a one-step reaction: tRNA + amino acid + ATP → aminoacyl-tRNA + AMP + PPi (2) By assuming 1) ATP saturation (a high affinity between the aminoacyl-tRNA synthetase and ATP, such that ATP binding is not limiting), 2) random-ordered substrate-binding, and 3) competition between tRNA isoacceptors that share the same aminoacyl-tRNA synthetase, the following reaction scheme can be described for an aminoacyl-tRNA synthetase with two competing tRNA isoacceptors: Figure 2: tRNA aminoacylation described using Michaelis-Menten enzyme kinetics. tRNA aminoacylation reaction scheme assumes ATP saturation, random-ordered substrate-binding, and competition between two tRNA isoacceptors sharing the same aminoacyl-tRNA synthetase.
where E is the aminoacyl-tRNA synthetase, A is the amino acid, T is the unaminoacylated tRNA, and P is the aminoacyl-tRNA. The Michaelis-Menten constants K A and K T describe the affinity between the aminoacyl-tRNA synthetase and its substrates -amino acid and tRNA, respectively.
According to Michaelis-Menten enzyme kinetics, the rate of formation of the i-th aminoacyl-tRNA can be described as: where k cat is the rate of catalysis of the aminoacyl-tRNA synthetase.
Convenience kinetics [2] can be used to generalize the rate of aminoacyl-tRNA formation for aminoacyl-tRNA synthetases with n tRNA isoacceptors: v aminoacylation,i = k cat [E] [A] K A + [A] [Ti] K T ,i 1 + n j=1 [Tj ] K T ,j (4)

Represented Molecules
E. coli K-12 MG1655 has 23 documented aminoacyl-tRNA synthetase genes [3], of which 22 are modeled and 1 (lysU ) is excluded. Lysine tRNAs are reported to be aminoacylated by two aminoacyl-tRNA synthetases: lysS -which is modeled -and lysU, which is expressed at elevated temperatures and was thereby excluded from our model of cell growth at 37°C.

Codon-to-tRNA Anticodon Interactions
The interactions of mRNA codons to tRNA anticodons were described by:

Initiator and Elongator Methionine tRNAs
• E. coli has two types of methionine tRNAs: the initiator methionyl-tRNA used for translation initiation and the elongator methionyl-tRNA used for translation elongation. Initiator methionyl-tRNAs are formylated by Fmt, and are thereby distinguished from elongator methionyl-tRNAs by the cell's translation machinery [3]. • The four reported initiator methionine tRNAs are: metV, metW, metY, and metZ; the two reported elongator methionine tRNAs are: metT and metU [3].
To enable matrix operations, the codon-to-tRNA anticodon interactions for each amino acid family have been summarized into the following Boolean matrices:

Represented Codons
Of the 64 possible codons, all 61 sense codons were mechanistically modeled and the role of the 3 stop codons were assumed by instantaneously disassembling active ribosomes and releasing the finished protein product after the final sense codon is processed.
Additionally, start codons were distinguished from elongating AUG codons in order to assign initiator and elongator Methionine tRNAs to their corresponding codons.

Utilization of Aminoacyl-tRNAs by Ribosomes
The utilization of aminoacyl-tRNAs by ribosomes (representing the release of tRNAs after amino acid transfer by elongating ribosomes) was represented as a one-step reaction: where aminoacyl-tRNA i is the aminoacylated form of the i-th tRNA isoacceptor, and ribosome j is an active ribosome with an open A site at the j-th codon on the mRNA transcript.
For codons that can be read by multiple aminoacyl-tRNA isoacceptors, the distribution of codon reading load across the different aminoacyl-tRNA isoacceptors was assumed to match the relative abundances of those aminoacyl-tRNAs. To do so, an (n × m) matrix T was described for each amino acid family: where n is the number of tRNA isoacceptors, m is the number of codons, and T * i is the abundance of the i-th aminoacyl-tRNA isoacceptor in the amino acid family.
Then, the aminoacyl-tRNA abundances were masked by codon-to-tRNA anticodon interactions via element-wise multiplication: where M is one of the matrices described by Equations 5 through 24 (depending on the amino acid family).
Finally, the columns of D were normalized to yield the codon reading load distributionD, which has elements:D During optimization of aminoacyl-tRNA synthetase kinetic parameters (Section 2) and simulation (Section 3.1), a vector describing the codon reading rate, v codon is prepared. The dot product ofD and v codon yields the utilization rate of each aminoacyl-tRNA isoacceptor by active ribosomes.

Optimization of Aminoacyl-tRNA Synthetase Kinetic Parameters
Parameter optimization of the aminoacyl-tRNA synthetase kinetic parameters was described as objective minimization problems, where each aminoacyl-tRNA synthetase was treated as an independent optimization problem.
The developments related to optimization of aminoacyl-tRNA synthetase kinetic parameters occurred in the following files of the E. coli Model repository:  The optimized aminoacyl-tRNA synthetase kinetic parameters are listed in Table S2.

Approach
For each aminoacyl-tRNA synthetase E, the rate of tRNA aminoacylation v of the i-th tRNA was described by Michaelis-Menten enzyme kinetics as described in Section 1.1, Equation 4.
where [E] is the concentration of the aminoacyl-tRNA synthetase, [A] is the concentration of the amino acid, [T ] is the concentration of total tRNA (ie. both aminoacylated and unaminoacylated forms), f is the fraction of [T ] that exists in the unaminoacylated form, K A and K T are the Michaelis-Menten constants describing the binding affinity between the aminoacyl-tRNA synthetase and its substrates -amino acid and tRNArespectively, and n is the number of tRNA isoacceptors that aminoacyl-tRNA synthetase E interacts with.
During optimization, the free parameters that were solved for were:

Objective Function
The objective function for each aminoacyl-tRNA synthetase problem takes the following form: where e ss is the steady-state error, R is the regularization term, and B is the bounds penalty. The weights, w r and w b , describe the significance of R and B, respectively.

Steady-State Error
The steady-state error, e ss , describes the relative error between the rates of tRNA aminoacylation v A and aminoacyl-tRNA utilization by ribosomes v D , which are anticipated to be equal at steady state according to the tRNA aminoacylation cycle depicted in Figure 1.
The rate of aminoacylation of the i-th tRNA isoacceptor is described by Equation 31 (reproduced below):     Table 7: Total (both aminoacylated and unaminoacylated forms) concentrations of tRNA isoacceptors used as constants in the optimization of aminoacyl-tRNA synthetase kinetic parameters. Data were retrieved from the Parca of the E. coli Model [1].
To describe the utilization of aminoacyl-tRNAs by ribosomes, the anticipated rate of codon reading was calculated by describing the proteome in units of number of codons read N by time t. At the start of the cell cycle, the proteome can be described as the result of having read N 0 codons by a certain time t = 0: By the end of the cell cycle (t = τ , where τ is the doubling time), the proteome will have doubled: Together, equations 34 and 35 describe exponential growth, where the change in the number of codons read as time goes on can be described by the growth rate r: and the rate of change in the number of codons read by a certain time t is: To describe the rate of codon reading, v codon , the number of codons read during time step ∆t is described as the difference in the number of codons read by time t and by a small step later t + ∆t divided by the time step: To express the rate of codon reading in units molar concentration per time, the right-hand side of Equation 38 is divided by Avogadro's number (N A ) and the cell's volume at time t (V (t)): Since the cell's volume is not represented directly in the E. coli Model, but rather calculated as the cell's mass (M (t), which is represented directly) divided by the cell's density (ρ), the rate of codon reading is described as: Substituting Equation 37 yields: where N 0 is the number of codons read to create the initial proteome and M 0 is the initial mass of the cell.
To generalize Equation 41 for the reading rate of the i-th codon type, N 0 is replaced with the number of times codon i is read to create the initial proteome:  Table 8: Codon reading rates used as constants in the optimization of aminoacyl-tRNA synthetase kinetic parameters. Data were calculated according to Equation 42, where N i,0 , M 0 , ρ, N A , and τ were retrieved from the Parameter Calculator (Parca) of the E. coli Model [1].
To convert the reading rate of codons into the rate of aminoacyl-tRNA utilization by ribosomes, the aminoacyl-tRNA abundance matrix T is populated as described by Equation 26.
where T i is the total amount of the i-th tRNA isoacceptor (ie. both aminoacylated and unaminoacylated forms), (1 − f i ) is the fraction of [T i ] that exists in the aminoacylated form, and n is the number of tRNA isoacceptors that interact with the aminoacyl-tRNA synthetase described by this objective function.
Then, the appropriate codon-to-tRNA anticodon mapping matrix M is retrieved from Equations 5 through 24 (depending on the amino acid family) and used to calculate D according to Equation 27: where ⊗ is the element-wise multiplication operation.
Normalizing the columns of D according to Equation 28 yields the codon reading load distribution matrix D, which has elements:D where n is the number of tRNA isoacceptors that interact with the aminoacyl-tRNA synthetase described by this objective function.
Finally, the codon reading rates (from Equation 41) for the codons listed along the columns of matrix M are assembled into a vector v codons . The dot product ofD and v codons yields the rate of utilization each aminoacyl-tRNA isoacceptor by ribosomes.
from which the utilization rate of the i-th aminoacyl-tRNA isoacceptor is the i-th element of v D .
Taken together, the sum of the square of the relative errors between the tRNA aminoacylation and aminoacyl-tRNA utilization rates of all tRNA isoacceptors yields the steady-state error: where n is the number of tRNA isoacceptors that interact with the aminoacyl-tRNA synthetase described by this objective function.
When testing candidate kinetic parameter solutions in simulations, it was observed that the best solutions continued to achieve steady-state (i.e. v A matched v D ) even at low aminoacyl-tRNA synthetase concentrations. Consequently, a second term was incorporated to account for the minimum aminoacyl-tRNA synthetase concentration, [E min ]: where v A,avg,i (previously v A,i in Equation 47) is the rate of aminoacyl of the i-th tRNA isoacceptor at the average aminoacyl-tRNA synthetase concentration ([E avg ]) and v A,min,i is the rate of aminoacylation of the i-th tRNA isoacceptor at [E min ]: where where c = 2 is the number of growth conditions used in all parameter optimization problems (the E. coli Model's representation of M9 Minimal Media supplemented with 0.4% Glucose, and M9 Minimal Media supplemented with 0.4% Glucose and amino acids).

Regularization
As an underdetermined system, the aminoacyl-tRNA synthetase kinetic parameter optimization problem generated many solutions, some with needlessly large estimates of kinetic parameters. Thus, a regularization term was incorporated to prefer smaller kinetic parameters that could achieve the same output.
Inspecting Equation 49 revealed that the K T 's (the Michaelis-Menten constant describing the binding affinity between tRNAs and their aminoacyl-tRNA synthetases) would be the most appropriate parameters to apply the regularization to because their primary role is to distribute the aminoacylation capacity across the different tRNA isoacceptors.
Thus, the sum of the K T 's was scaled with weight w r = 10 −9 and minimized, as described in Equation 32, where R is: where n is the number of tRNA isoacceptors that interact with the aminoacyl-tRNA synthetase described by this objective function.

Bounds Penalty
As a fraction, f must range between 0 and 1; this requirement is described by the constraints of the optimization: where n is the number of tRNA isoacceptors that interact with the aminoacyl-tRNA synthetase described by this objective function.
As an underdetermined system, the optimization problem generated many solutions, some with f values close to 0 or 1, which led to unstable simulations. Thus, a bounds penalty was incorporated to prefer f values that could achieve the same output with greater stability. A 5% buffer from either extreme was described the term w b · B in Equation 32, where w b = 10 −9 is the weight and B is: where n is the number of tRNA isoacceptors that interact with the aminoacyl-tRNA synthetase described by this objective function, and B(f ) has the shape: Figure 3: Profile of bounds penalty applied to f , the fraction of unaminoacylated tRNA.

Complete Form
Taken together, the complete form of the objective function for the set of kinetic parameters that describe one aminoacyl-tRNA synthetase is: where n is the number of tRNA isoacceptors that interact with the aminoacyl-tRNA synthetase described by this objective function, c = 2 is the number of growth conditions used in all parameter optimization problems (the E. coli Model's representation of M9 Minimal Media supplemented with 0.4% Glucose, and M9 Minimal Media supplemented with 0.4% Glucose and amino acids), v A,avg and v A,min are the rates of tRNA aminoacylation at the average and minimum aminoacyl-tRNA synthetase concentrations, respectively (as described by Equation 49), v D is the rate of aminoacyl-tRNA utilization by ribosomes (as described by Equation 46), f avg and f min are the fractions of the total tRNA T that exist in the unaminoacylated form at [E avg ] and [E min ], respectively, and K T is the Michaelis-Menten constant describing the binding affinity between the aminoacyl-tRNA synthetase and the unaminoacylated tRNAs.

Random Initialization
To explore the solution space, optimization of the minimization problem described in Equation 54 was performed many times with random initialization of the initial solution. For simplicity, tRNA isoacceptors that have the same codon interaction patterns (according to Equations 5 through 24) are assigned a single, shared K T .
x is a random number selected from a uniform distribution between 1 and 6. tRNA isoacceptors that have the same codon interaction patterns (according to Equations 5 through 24) are assigned a single, shared K T .
/* Initialize f 's */ for f avg,i do f avg,i ← a random number selected from a uniform distribution between 0.1 and 0.9 for f min,i do f min,i ← a random number selected from a uniform distribution between 0.1 and 0.9 /* Estimate a compatible k cat in each growth condition */ for condition, c do /* Describe saturation of substrates */ sat A = [Ac] /* Initialize the k cat at the maximum estimate */ k cat ← the maximum of k cat,1 and k cat,2 Result: One set of kinetic parameters has been randomly initialized.

Parametric Sweep on the Minimum Aminoacyl-tRNA Synthetase Concentration
To assess the impact of the minimum aminoacyl-tRNA synthetase concentration, a parametric sweep was performed on [E min ]. To inform the range of this sweep, a sample set of 200 simulations (100 in M9 Minimal Media + 0.4% Glucose, and 100 in M9 Minimal Media + 0.4% Glucose + Amino Acids) was run in the prior E. coli Model [1], which assumes unlimited tRNA aminoacylation, to gauge the dynamic range of aminoacyl-tRNA synthetases. Each set of 100 simulations was composed of 20-generation runs initialized at 5 random seeds. The commands used to run these sample sets were: python runscripts/manual/runSim.py --no-trna-charging --translation-supply -v condition 0 0 -g 20 -i 5 python runscripts/manual/runSim.py --no-trna-charging --translation-supply -v condition 2 2 -g 20 -i 5  Table 9) indicate the values of the minimum aminoacyl-tRNA synthetase concentration [E min ] used in Equation 49. Note that Figure 4 and Figure S3 are the same.
In order to prepare the tRNA aminoacylation system to remain robust at low levels of aminoacyl-tRNA synthetase (which can be variable between cells), the [E min ] was allowed to be equal to or less than the observed minimum (labeled '5' in Figure 4). To do so, the distance between the observed minimum ('5') and 0 µM (which is the lowest feasible value) was divided into five equal parts for each aminoacyl-tRNA synthetase. Of these parts, the values corresponding to the highest four values (labeled '2', '3', '4', and '5' in Figure 4) were identified and used as the four values of [E min ] in the parametric sweep on the minimum aminoacyl-tRNA synthetase concentration ( Table 9). The aminoacyl-tRNA synthetase concentrations corresponding to '1' was found to lead to extremely low intracellular concentrations that could not generate viable aminoacyl-tRNA synthetase kinetic solutions, and was therefore not used.
Next, the aminoacyl-tRNA synthetase kinetic solutions with the minimum objective value was identified for each aminoacyl-tRNA synthetase at each [E min ] sweep value. The k cat s corresponding to these minimumobjective solutions are depicted Figure 3H of the main text.
Finally, to select the best performing aminoacyl-tRNA synthetase kinetic parameters, the behavior of these minimum-objective solutions were assessed in test simulations of a nutritional downshift experiment from rich (M9 Minimal Media + 0.4% Glucose + Amino Acids) to poor (M9 Minimal Media + 0.4% Glucose) media. These tests revealed that the aminoacyl-tRNA synthetase kinetic parameter set corresponding to [E min ] sweep value '4' demonstrated balanced exploration of aminoacylated fraction values across all tRNAs. For example, the aminoacyl-tRNA synthetase kinetic solutions generated from '2' caused aminoacylation fractions to remain unresponsive to the downshift test, while those generated from '5' caused aminoacylation fractions to fall close to 0%.

Comprehensive Algorithm
Optimization was performed using SciPy's minimize function with the Powell method. At least 100 iterations were performed for each aminoacyl-tRNA synthetase problem. for each value of E x,min,c (listed in Table 9) do for iteration do Generate a random initial solution (Algorithm 1. Call SciPy minimize, using the objective function described by Equation 54 and the constant parameters. Save solution with its objective value. Identify the minimum-objective solution.
Result: Candidate solutions of aminoacyl-tRNA synthetase kinetic parameters are generated.

Simulation
The model described in Section 1 was implemented into the polypeptide elongation model of the E. coli Model published in 2020 [1]. A summary-level flowchart of the computational steps that take place within each time step is shown in Figure 1C of the main text, where the model described in Section 1 impacts the polypeptide elongation model during the estimation of ribosome steps, determination of feasibility, reconciliation of kinetic and sequence solutions, and the update of molecule abundances 3.4.

Estimation of Ribosome Steps
To estimate the anticipated codon reading rate during each time step, the number of times each codon type appears in the upcoming codon sequence of all ribosomes is computed.

Determination of Feasibility
To determine the number of feasible ribosome steps, the tRNA aminoacylation and codon reading cycle was described as a system of ODEs and the initial value problem was solved at each time step by the Runge-Kutta method. The resulting solution yields the kinetically feasible number of codon reading events that can occur in each time step.

Algorithm 4: Determination of Feasibility
/* Each i-th aminoacyl-tRNA synthetase interacts with its amino acid. */ Input : n E,i , Number of available i-th aminoacyl-tRNA synthetase in the cell Input : n A,i , Number of available amino acids cognate to the i-th aminoacyl-tRNA synthetase in the cell Input : k cat,i , Rate of catalysis of i-th aminoacyl-tRNA synthetase Input : K A,i , Michaelis-Menten constant describing the binding affinity between the i-th aminoacyl-tRNA synthetase and its cognate amino acid /* Each j-th tRNA exists in the unaminoacylated and aminoacylated forms.
*/ Input : n T d ,j , Number of available unaminoacylated j-th tRNA in the cell Input : n T a ,j , Number of available aminoacylated j-th tRNA in the cell Input : K T,j , Michaelis-Menten constant describing the binding affinity between the j-th unaminoacylated tRNA and its corresponding aminoacyl-tRNA synthetase /* Codon reading rate of each codon type is described in one vector. */ Input : v codons /* Cellular parameters. */ Input : v, Cell volume Input : N A , Avogadro's number Input : ∆t, Time step /* Convert Michaelis-Menten constants to units of number of molecules. */ Dnm where N is the number of rows in D. /* Describe tRNA aminoacylation and codon reading as a system of ODEs */ for each aminoacyl-tRNA synthetase, i do 1. Describe the rate of aminoacylation for each j-th tRNA (that is a substrate of the i-th aminoacyl-tRNA synthetase) according to Equation 4.
where n is the number of tRNA isoacceptors for the i-th aminoacyl-tRNA synthetase. 2. Assemble the codon reading rate for the codons read by aminoacyl-tRNA synthetase i from v codons into a vector, v codon,i . 3. Describe the rate of ribosome utilization of all aminoacyl-tRNA isoacceptors of the i-th aminoacyl-tRNA synthetase according to Equation 29. v utilization =D i · v codon,i /* Retrieve the kinetically feasible number of codon reading events */ Call the solve ivp program from SciPy to solve the initial value problem during ∆t. Stochastically round the resulting number of codon reading events that occurred during ∆t for each codon type. Result: The kinetically-feasible number of reading events for each codon type is estimated.
By using the abundance of aminoacyl-tRNAs to calculate the codon reading load distribution matrixD, the distribution of codon reading load across different aminoacyl-tRNA isoacceptors is proportional to the relative abundances of aminoacyl-tRNAs.

Motivation
In the updated model reported in this study, two different calculations each yield an estimation of the number of codons read during a given time step: 1. a kinetics-based model of the tRNA aminoacylation cycle that obeys the kinetic capacity of aminoacyl-tRNA synthetases to produce aminoacyl-tRNAs (described in Section 3.2), and 2. a sequence-based model of ribosome procession that obeys the codon sequence order that ribosomes encounter on mRNAs (which is the method used in the prior E. coli Model [1]).
To demonstrate each of their outputs during time step ∆t in a two-codon system: 1. The kinetics-based model yields the reading rate of each codon. For example, • codon A was read at v A codons per second • codon B was read at v B codons per second Multiplying by the time step t results in the number of times each codon type was read: • codon A was read n A times, where n A = v A · ∆t • codon B was read n B times, where n B = v B · ∆t 2. The sequence-based model yields the number and order of codon reading events. For example, • Ribosome 1 read codon A, then codon B, then codon A; • Ribosome 2 read codon B, then codon A, then codon A; • ... and so on until Ribosome n read codon A, then codon B, then codon B.
Collectively, the ribosomes read codon A n A times and codon B n B times.
Each model approaches the calculation of the number of codons read during t from different perspectives -the kinetic capacity of aminoacyl-tRNA synthetases and the sequence order of codons -and thereby, their estimations of n A and n B may not agree. In the E. coli Modeling framework, both perspectives are valuable and must be honored. As a result, the purpose of the Reconciliation Program is to reconcile any disagreements between the kinetics-based and sequence-based solutions.

Tools
To reconcile disagreements in the number of codons read, two points of flexibility were identified and used as tools for the reconciliation task: the flexibility offered by the codon sequence and the tRNA pools.

Flexibility Offered by the Codon Sequence
For each ribosome, the position of the ribosome and the codon sequence of the mRNA could be used to survey the identities of the surrounding codons. For instance, if the sequence-based model overestimated (compared to the kinetics-based model) the number times codon A was read by 1 instance -and if codon A was the final codon read by the i-th ribosome -then the disagreement could be resolved by reversing the last codon reading event by ribosome i. Flexibility offered by the tRNA pools For each tRNA isoacceptor, the pools of aminoacylated and unaminoacylated tRNAs could be used to absorb differences in codon reading events. For instance, the profile of ribosome positions on mRNAs from the previous example may result in a maximum of n B,max instances of codon B readings -no more because that would require reading codon A to reach the next instance of codon B (according to the sequence order of codons) and thereby compromising the kinetic limit of codon A readings.
If the kinetics-based solution n B is greater than n B,max , then the kinetics-based model has overestimated (compared to the sequence-based model) the number of times codon B was read. If the overestimation occurred by 1 codon B reading event -and if at least 1 copy of the tRNA isoacceptor that interacts with codon B exists in the unaminoacylated form -then the disagreement can be resolved by returning the unaminoacylated tRNA to its aminoacylated form, which reverses the final codon reading event in the kinetics-based model.
If no copies of the tRNA isoacceptor exist in the unaminoacylated form, then the assumption of a constant codon reading rate in the kinetics-based model (Algorithm 3) must have been violated. This can be seen in a linear representation of the tRNA aminoacylation cycle. Figure 6: Flexibility offered by tRNA pools. In this scenario, the kinetics-based model proposes that an unaminoacylated tRNA (yellow rectangle) can be aminoacylated three times and de-aminoacylated two times during a particular time step. In contrast, the sequence-based model has calculated that the final codon-reading event was not feasible according to the sequence order of codons on mRNAs (green X). This overestimation can be resolved by reversing the last codon B reading event (green X) by the ribosome (blue), and eliminating all downstream reactions, such as the final tRNA aminoacylation reaction by the aminoacyl-tRNA synthetase (E).
In Figure 6, the final tRNA aminoacylation event assumed the existence of an unaminoacylated tRNA. However, the reversal of the final codon reading event nullifies the existence of the unaminoacylated tRNA and all downstream reactions.
The result of resolving kinetics-based model overestimations through reversing codon-reading events (and downstream reactions) is a solution of the number of tRNA aminoacylation and codon reading events that obey the order of codons encountered by ribosomes during time step ∆t.

Strategy
Whereas utilizing the flexibility offered by the codon sequence preserves the kinetic solution, adjusting the tRNA pools can result in a reduced kinetic solution than originally intended if tRNA aminoacylation events are reversed (Section 3.3.2). As such, the two tools were applied in a preferential manner: resolve as many disagreements as possible using the flexibility of codon sequences, then resolve any remaining disagreements using the flexibility of tRNA pools.
Stage 1: Using the flexibility of the codon sequence Disagreements in the number of codons read by the kinetics-based and sequence-based models can occur in both orientations: the sequence-based solution may overestimate the reading of some codons and underestimate the reading of other codons. When resolving these disagreements using the flexibility of the codon sequence, resolving sequence-based underestimations first (by processing ribosomes towards the 5' end) followed by the overestimations (by processing ribosomes towards the 3' end) ensures that the kinetic limit of codon reading events is always honored. In the scenario depicted in Figure 7, the kinetics-based solution determined that two reading events each for codons A (green) and B (yellow) occurred during time step ∆t. Meanwhile, the sequence-based solution determined that three codon A reading events and one codon B reading event occurred. Resolving the sequence-based underestimation first (codon B is read 1 time in the sequence solution, which is less than its 2 reading events in the kinetic solution), followed by the sequence-based overestimation next (codon A is read 3 times in the sequence solution, which is greater than its 2 readings in the kinetic solution) ensures that the resulting solution honors the kinetic limit of codon reading events (codon A is read 2 times in both solutions and codon B is read 1 less time in the sequence solution than the kinetic solution).

Selecting the Best Compromise
When using the flexibility of the codon sequence to reconcile disagreements in the number of codons read between the two models, two key characteristics of codon sequences made it possible to arrive at many different solutions: 1. Tens of thousands of ribosomes: At each time step, there are 10,000s of active ribosomes processing on mRNAs, each with a unique sequence of codons flanking its current location. As a result, for each disagreement of codon i, there are many candidate ribosomes that can undergo the process depicted in Figure 7.
2. Presence of 'in-between' codons: In some cases, the next available instance of codon i is not immediately accessible (as it was for the ribosome in Figure 7, which took +1 and -1 steps). Rather, the ribosome may need to take multiple steps (as depicted in Figure 8) to arrive at the next instance of codon i. The presence of the 'in-between' codons -the two instances of codon C (pink), for examplemay introduce further disagreements between the two models. Consequently, encountering further disagreements while pursuing global reconciliation is anticipated. The presence of 'in-between' codon C (pink) introduces new disagreements that did not exist at the start of the reconciliation.
As a result of these characteristics, each random selection of a ribosome to use (among many candidate ribosomes that could be used for resolving a codon disagreement) affects the landscape of ribosome positions and thereby impacts the options that are available for resolving the next codon disagreement. To resolve as many disagreements as possible in Stage 1, an iterative approach that enables sampling of the reconciliation solution space was developed. To compare candidate solutions, it was necessary to design a metric that would measure a solution's degree of success. This metric, called the compromise C, reflects the total number of disagreements remaining in a solution: where k i is the number of times codon i is read according to the kinetics-based solution, s i is the number of times codon i is read according to the sequence-based solution, and n is the number of codon types. A compromise of C = 0 indicates a completely reconciled state.
Stage 2: Using the flexibility of tRNA pools Empirically, it was found that 75% of initial disagreements could be resolved in Stage 1. To resolve the remaining disagreements, the flexibility offered by the tRNA pools was used.
Due to the design of resolving sequence-based underestimations first (by processing ribosomes towards the 5' end) followed by the overestimations (by processing ribosomes towards the 3' end), any disagreements remaining by Stage 2 were sequence-based underestimations. These underestimations were resolved by reversing the appropriate number of final codon reading events (and any downstream reactions) as depicted in Figure 6.
The Reconciliation Program algorithm has been summarized as a flowchart in Figure S1, and its context in the polypeptide elongation model is indicated in Figure 1C of the main text.

Update of Molecules
Due to the fast turnover of aminoacyl-tRNAs, each tRNA undergoes multiple cycles of aminoacylation during the simulated time step. As a result, net changes to the tRNA pools were reflected when updating the number of molecules at the end of each time step.
To do so, the net number of aminoacylation events ∆n that occurred during time step ∆t was determined by subtracting the number of aminoacyl-tRNAs at the end of the time step n 1 from that at the beginning of the time step n 0 : ∆n = n 1 − n 0 where ∆n > 0 indicates net aminoacylation and ∆n < 0 indicates net utilization of aminoacyl-tRNAs.
Each net aminoacylation event increments the number of aminoacyl-tRNAs, decrements the number of unaminoacylated tRNAs and amino acids, and converts one molecule of ATP to PPi and AMP.
where i indicates the tRNA isoacceptor type and its corresponding amino acid, and ∆n i > 0.
Each net aminoacyl-tRNA utilization event indicates an incorporation of an amino acid residue from the aminoacyl-tRNA to the nascent polypeptide and releases a proton molecule. For clarity, the net number of aminoacyl-tRNA utilization events for the i-th aminoacyl-tRNA is described as ∆d i , where ∆d i = −∆n i .
As a result of reflecting net changes to the tRNA pools, a portion of the total number of times the i-th codon is read (n codon,i ) is represented as aminoacyl-tRNA utilization events (∆d i ) while the remaining portion is represented as direct incorporation of residues (n direct,i ) from the amino acid pool as in the prior E. coli Model [1].
Each direct incorporation of residue produces a water molecule:

N-terminal Cleavage of Initial Methionines
In each time step, the maximal number of cleavage events by MAP (n M AP,max ) was determined by assuming the maximal reaction rate according to Michaelis-Menten enzyme kinetics (described in Section 1.4 during the time step duration ∆t.
where v M AP,max is described in Equation 30, N A is Avogadro's number, and V is the cell volume.
At each time step, the maximal number of cleavage events by MAP (n M AP,max ) was compared with the candidate number of cleavage events (n candidates ) to determine how many cleavage events could occur. MAP substrates that are successfully cleaved proceed to translation termination. Input : n candidates , number of terminating ribosomes that are synthesizing a MAP substrate (Table  3 lists all documented substrates of MAP.) Input : V , Cell volume Input : N A , Avogadro's number Input : ∆t, Time step /* Determine kinetic capacity of MAP */ v M AP,max = k cat · n M AP n M AP,max = v M AP,max · V · N A · ∆t /* Compare the maximal and candidate numbers of cleavage events */ if n candidates ≤ n M AP,max then Proceed all terminating ribosomes to termination. else Determine the number of substrates that cannot be cleaved, n cannot cleave . n cannot cleave = n candidates − n M AP,max Randomly select n cannot cleave ribosomes. Prevent the n cannot cleave ribosomes from termination in this time step. Result: The number of successfully terminating ribosomes is limited by the maximal kinetic capacity of MAP.

Simulations
Prior to running simulations, the Parameter Calculator (Parca) was run to initialize simulation data, including the parameter optimization program described in Section 2. The command used to run the Parca was: python runscripts/manual/runParca.py --optimize-trna-charging-kinetics sims where sims is the name of the simulation output directory, and the --optimize-trna-charging-kinetics option enables the parameter optimization program.

Methods
Several figure panels investigate comparisons made across different simulations. The table below summarizes which simulations were used to perform the analysis presented in each figure panel.  Table 12: Simulations used to perform the analyses presented in each figure of the main text. Figure 2 • Panel A: Doubling times of cells simulated in the Prior Model (Variant 1, gray) and the Updated Model when using measured aminoacyl-tRNA synthetase k cat s (Variant 2, blue) were obtained by measuring the time between cell division events. The anticipated doubling time of 44 minutes was informed [10,11]. The command that generated the panel was:

Analysis Methods for Panels in
python models/ecoli/analysis/variant/trna synthetase kinetics variant doubling.py sims • Panel B: Mass accumulation during the cell cycle in a representative cell simulated in the Updated Model when using measured aminoacyl-tRNA synthetase k cat s was obtained from Variant 2, Seed 11, Generation 0, which was selected for exhibiting a doubling time that is closest to the population average. Protein mass includes all protein monomers and complexes. Non-protein mass includes DNA, RNA, and small molecules. The command that generated the panels was: python models/ecoli/analysis/cohort/mass cell cycle.py -v 2 sims • Panel C: Mass accumulation during each cell cycle was obtained from cells simulated in the Updated Model when using measured aminoacyl-tRNA synthetase k cat s (Variant 2). Protein mass includes all protein monomers and complexes. Non-protein mass includes DNA, RNA, and small molecules. The command that generated the panels was: python models/ecoli/analysis/cohort/mass cell cycle.py -v 2 sims • Panel E: The simulated proteome was obtained from simulations in the Prior Model (Variant 1, y-axis) by taking the average number of molecules across all time steps, adding 1 (to view molecules with an average of 0 molecules on the log-log plot), and taking the log base 10. The measured proteome was obtained from [12]. Of the 4,307 proteins in the model, 2,141 were measured in the proteomics dataset and are therefore included in the analysis. The command that generated the panel was: python models/ecoli/analysis/cohort/trna synthetase average.py -v 1 sims • Panel F: The simulated distribution of the number of molecules was obtained from simulations in the Prior Model (Variant 1, blue histogram) by recording the molecule abundance at every time step. The measured distribution was obtained from [13], which was performed in slower-growing cells with doubling times of 150 minutes. To be comparable with our faster growing cells (44 minutes doubling time), the measured distributions were scaled up to match the observed mean number of molecules for each aminoacyl-tRNA synthetase while preserving the shape of the measured distributions. The command that generated the panel was: python models/ecoli/analysis/cohort/trna synthetase distribution.py -v 1 sims Figure 3 •  Table 13. The command that generated the panel was:

Analysis Methods for Panels in
python models/ecoli/analysis/variant/trna charging validation.py sims • Panel F The simulated intracellular concentrations of tRNA isoacceptors was obtained from simulations in the Updated Model when using the optimized aminoacyl-tRNA synthetase kinetic parameters (Variant 0, y-axis) by taking the average concentration of molecules across all time steps. The measured intracellular concentrations were obtained from [14]. Each tRNA isoacceptor is represented by a gray dot; the mapping between tRNA isoacceptors reported in [14] to the tRNA isoacceptors modeled in our simulations is reported in Table 14. The command that generated the panel was: python models/ecoli/analysis/variant/trna charging validation.py sims • Panel G: The distribution of total tRNA was obtained my monitoring the number of all tRNA molecules (as a single sum) across all time steps. The doubling time value (x position of the violin plot) was calculated by taking the average of all simulated doubling times. The molecular abundance of tRNAs across the range of doubling times shown on the x-axis was taken from [15]. The command that generated the panel was: python models/ecoli/analysis/variant/trna charging validation.py sims • Panel H: The range of candidate k cat values resulting from the aminoacyl-tRNA synthetase kinetic parameter optimization (blue bars) and the selection of the best candidate (blue dots) are described in Section 2. The range of measured k cat values (gray bars) was observed from Table S1, and the Jakubowski and Goldman-estimated lower limit and the in vitro measurement they referenced were obtained from [16]. The ranking of aminoacyl-tRNA synthetases by their agreement with other reports was calculated by subtracting the higher of the highest measurement (top edge of gray bars) or Jakubowski and Goldman's estimation (gray triangles) by the optimized k cat (blue dots). The command that generated the panel was: python models/ecoli/analysis/parca/trna synthetase kcats.py sims   Pro2  proL  Pro3  proM  Ser1  serT  Ser2  serU  Ser3  serV  Ser5  serW, serX  Thr1  thrV  Thr2  thrW  Thr3  thrT  Thr4  thrU  Trp  trpT  Tyr1  tyrT, tyrV  Tyr2  tyrU  Val1 valT, valU, valX, valY, valZ Val2A valW Val2B valV Table 14: Mapping of tRNA isoacceptors from the measurements reported by Dong and others [14] to this study. Figure 4 • • Panel I: The sensitivity threshold was calculated by determining the highest HisRS concentration that was associated with a ribosome elongation rate of 17 amino acids per second per ribosome or less, which approximately characterizes where the linear decrease in ribosome elongation rate meets the stable elongation of ribosomes at 17.5 amino acids per second per ribosome. The minimum ribosome rate was determined by identifying the minimum y-value. The portion of time with tRNA aminoacylation limitation was calculated as the portion of time steps displaying a ribosome elongation rate that was less than the minimum observed ribosome elongation rate in the simulations performed in the Updated Model when using the optimized aminoacyl-tRNA synthetase kinetic parameters (gray). All columns were calculated from the simulations shown in Figure 4H of the main text. Figure 5 •  -Row 2: The aminoacylated fraction was calculated by dividing the total number of argininecharged arginine tRNAs by the total number of arginine tRNAs (both aminoacylated and not).

Analysis Methods for Panels in
-Row 3: The ribosome A site position along argA mRNA transcripts was monitored during simulations. The codon sequence is numbered from 0 (start codon) to 442 (the final sense codon). Codon identities are colored: CGU in blue, CGG in green, CGC in yellow, AGG in pink, and non-arginine codons in gray).
-Row 4: Ribosome accumulation on argA mRNA transcripts was calculated by dividing the number of active ribosomes on argA mRNAs by the number of argA mRNA transcripts. Note: volume exclusion of polysomes on shared transcripts is not yet represented in this model. For reference, our observation of up to 40 ribosomes per argA transcript in the measured ArgRS k cat simulation (right column) agrees with the length of the transcript (443 amino acids, or 1332 base pairs including the stop codon), which suggests that a maximum of 55 ribosomes would physically fit (when assuming an estimated ribosome footprint of 24 nucleotides [17]).
-Row 5: The premature ribosome termination rate was calculated by dividing the number of premature ribosome termination events from argA mRNAs in each time step by the time step length (2 seconds).
-Row 6: The number of ArgA monomers were monitored at each time step and were directly shown.
-Row 7: The number of ArgA hexamer (acetylglutamate synthase) were monitored at each time step and were directly shown. The period of time when the ArgA hexamer abundance was at 0 molecules was indicated by the red line.
-Row 8: Flux through the N-acetyl transfer reaction in the arginine biosynthesis pathway is shown as the moving average using 10-minute windows. The period of time when the N-acetyl transfer reaction had a flux of 0 mmol/gDCW/h was indicated by the red line.
-Row 9: The intracellular concentration of arginine was monitored at each time step and were directly shown. -Row 13: Total protein mass was monitored at each time step of the simulations and was directly shown.
The command that generated the panel was: python models/ecoli/analysis/variant/argRS kcat impact.py sims 4.2.6 Analysis Methods for Panels in Supporting Figure S1 The flowchart in Figure S1 depicts the algorithm described in Section 3.3, which has been developed in the following file: wholecell/utils/ trna charging.pyx. Figure S2 For Figure S2, please refer to the description of Figure 2F in Section 4.2.1. Figure S3 For Figure S3, the concentrations of aminoacyl-tRNA synthetases were observed in sample runs from the previous model that were generated by the commands:
Then, for each aminoacyl-tRNA synthetase, the distance between the minimum observed concentration and 0 µM was divided into five equal parts, of which the values corresponding to the four highest values (labeled '2' through '5') were used as the minimum aminoacyl-tRNA synthetase concentration described in Section 2.3.2. The command that generated the panel was:   python models/ecoli/analysis/cohort/premature ribosome termination.py -v 6 sims • Panel B: The locations of arginine codons were identified from the known argA codon sequence [3]. The command that generated the panel was: python models/ecoli/analysis/cohort/premature ribosome termination.py -v 6 sims • Panel C: Proteins with tandem CGGs were identified by searching for two CGG codons in-a-row in their codon sequences; the same criteria were used to identify proteins with tandem AGAs and AGGs. Proteins with no arginine codons were identified as proteins that do not use any of the six arginine codons (AGA, AGG, CGA, CGC, CGG, and CGU). The remaining 4,094 proteins in the proteome were classified as the group of proteins that contain arginine codons that are not tandem arrangements of CGGs, AGAs, or AGGs. The command that determined the number of proteins in each classification was: python models/ecoli/analysis/variant/rare arginine codons expression.py sims • Panel D: The groups of protein features investigated in this panel are composed of the same protein categories identified in panel C. For each protein, the fold change of its simulated expression was calculated as 1 + the mean number of molecules per cell in the measured ArgRS simulations divided by 1 + the mean number of molecules per cell in the optimized ArgRS simulations. The "No Arginines" group was taken as the reference distribution and described using its mean and standard deviation. For the statistical significance of the difference between the median of each distribution to the "No Arginines" group, the two-tailed p-value was calculated from the z-score of the median. To focus on the time when premature ribosome terminations were occurring (as read from Fig. 6A, row 3), generations 8 through 11 were used in this analysis (total of n = 30 cells in the optimized ArgRS k cat and n = 40 cells in the measured ArgRS k cat ). The command that generated the panel was: python models/ecoli/analysis/variant/rare arginine codons expression.py sims • Panel E: For each perturbation experiment, the number of molecules per cell was recorded from each time step of the simulations and the distributions are shown. The "CGG" simulations were taken as the reference distributions and described using their means and standard deviations. For the statistical significance of the difference between the means of each "CGU" distribution to its corresponding "CGG" distribution, the two-tailed p-value was calculated from the z-score of the median. For the "Measured ArgRS k cat , CGG" group: to focus on the time when premature ribosome terminations were occurring (as read from Fig. 6A, row 3 Figure S1. Reconciliation Program. The algorithm used in the Reconciliation Program is described in Section 3.3. Notes: a: Identifies sequence-based underestimations by processing ribosomes forwards toward the 5' end. b: Identifies sequence-based overestimations by processing ribosomes backwards toward the 3' end. c: Disagreements remaining in Stage 2 are sequence-based underestimations. The Reconciliation Program is shown in context to the polypeptide elongation model in Figure 1 of the main text. Figure S2. Comparison of distribution of aminoacyl-tRNA synthetase abundances. Distributions of protein abundances estimated from simulations of the Prior Model (n=100 cells, blue) and measured from [13] (yellow) are shown as probability densities. Cells were simulated in aerobic growth in M9 Minimal Media supplemented with 0.4% glucose at 37°C. Simulations of the Prior Model represent 10-generation long lineages initialized at 10 random seeds (total of n=100 cells). CysRS, AlaRS, and ThrRS are highlighted in Figure 2F of the main text.   [18,19] or the standard deviation [20]. All simulations in this panel were performed in the Updated Model when using the optimized aminoacyl-tRNA synthetase kinetic parameters. All cells were simulated in aerobic growth in M9 Minimal Media supplemented with 0.4% glucose at 37°C. Simulations in this figure represent 10-generation long lineages initialized at 15 random seeds (n=150 cells). For each protein, the fold change of its expression was calculated as a ratio between the expression in the measured ArgRS simulations to the optimized ArgRS simulations. Statistical significance of the difference between the median of each distribution to the "No Arginines" group (which was taken as the reference distribution) was determined from the two-tailed p-value, which was calculated from the z-score of the median: "CGG, CGG": z-score for median (0.06) is -2.83, p-value is 4.7 × 10 −3 ; "AGA, AGA": z-score for median (0.41) is -1.55, p-value is 1.2 × 10 −1 ; "AGG, AGG": z-score for median (0.72) is -0.37, p-value is 7.1×10 −1 ; and "Other Arginine Codons": z-score for median (0.44) is -1.42, p-value is 1.6 × 10 −1 . Error bars indicate the range of observed values, and the squares (and their labels) indicate the median. (E) Distribution of the number of ArgA molecules per cell when using the optimized (gray) or measured (blue) ArgRS k cat , and when using the original CGG codons at positions 153 and 154 or using CGU at those positions. Statistical significance of the difference between the mean of each distribution to the "CGG" simulations (which was taken as the reference distribution) was determined from the two-tailed p-value, which was calculated from the z-score of the mean: "Optimized ArgRS k cat ": z-score for mean (569.7) is -0.25, p-value is 0.80; and "Measured ArgRS k cat ": z-score for mean (276.4) is 15.22, p-value is 2.6 × 10 −52 . Error bars indicate the range of observed values, and the squares (and their labels) indicate the mean.   Table S1 -continued. Table S2. Optimized aminoacyl-tRNA synthetase kinetic parameters. For each aminoacyl-tRNA synthetase, the optimization approach described in Section 2 generated several candidate solutions, of which the one corresponding to the minimum objective value (reported here) was taken to be the best solution and used for the rest of the study. These optimized parameters are compared to the estimated lower limit of aminoacyl-tRNA synthetase activity reported in [16].   Table S4. Properties of the 114 genes exhibiting tandemly arranged CGG codons in their codon sequences. The fold change in expression between simulations using the optimized and measured ArgRS k cat are reported here and also depicted as a violin plot in Figure S5C. The gene abbreviation, EcoCyc ID, and gene name were taken from the EcoCyc database [3].  Table S4