Dynamics and growth rate implications of ribosomes and mRNAs interaction in E. coli.

Understanding how cells grow and adapt under various nutrient conditions is pivotal in the study of biological stoichiometry. Recent studies provide empirical evidence that cells use multiple strategies to maintain an optimal protein production rate under different nutrient conditions. Mathematical models can provide a solid theoretical foundation that can explain experimental observations and generate testable hypotheses to further our understanding of the growth process. In this study, we generalize a modeling framework that centers on the translation process and study its asymptotic behaviors to validate algebraic manipulations involving the steady states. Using experimental results on the growth of E. coli under C-, N-, and P-limited environments, we simulate the expected quantitative measurements to show the feasibility of using the model to explain empirical evidence. Our results support the findings that cells employ multiple strategies to maintain a similar protein production rate across different nutrient limitations. Moreover, we find that the previous study underestimates the significance of certain biological rates, such as the binding rate of ribosomes to mRNA and the transition rate between different ribosomal stages. Furthermore, our simulation shows that the strategies used by cells under C- and P-limitations result in a faster overall growth dynamics than under N-limitation. In conclusion, the general modeling framework provides a valuable platform to study cell growth under different nutrient supply conditions, which also allows straightforward extensions to the coupling of transcription, translation, and energetics to deepen our understanding of the growth process.

Introduction

Growth is a fundamental process of life. The study of cell growth has always been of great interest to the scientific community, especially with increasing technological advances that allow for precise measurements and modifications of the biochemical composition and gene expressions of cells [1]. Protein synthesis is directly linked to cell growth rate, in particular during the exponential growth phase. However, the overall cell growth is a complex process that involves other heavily regulated cell functions that are necessary for cell survival [2], [3]. Thus, understanding how cells distribute their resources for growth and survival under different conditions is crucial to obtain a complete picture of integrated cellular function.

One promising direction examines optimal resource allocation theory to provide a quantitative framework to study cell growth under contrasting nutrient supply conditions [4], [5]. The theory proposes that cell growth is the result of optimal resource allocation, which is usually embedded in simple models of the translation process. Some of the earliest developments hypothesize the so-called “the constant-efficiency (of ribosome action) hypothesis”, which suggests that, since protein production is the most limiting process involved in cell growth, cell should optimize the protein synthesis rate of ribosomes to keep it at an optimal maximum regardless of situation. This implies that for one cell to grow twice as fast as another, it then must have twice the number of ribosomes [2], [6], [7]. This concept is an underlying principle behind “the growth rate hypothesis”, an important component in the study of ecological stoichiometry, which uses mass balance principles to link growth to organismal biochemistry and elemental composition [8], [9].

A version of the constant-efficiency hypothesis takes the form of a simple system of differential equations involving protein translation by ribosomes and self-replicating ribosomes [10]. An important implication of the model is that the protein to RNA ratio is linearly proportional to the growth rate, which is supported experimentally [4], [11]. While the translation rate may be relatively constant across different growth rates, it has been noted that there maybe alternative mechanisms that cells employ under different scenarios, especially during nutrient shifts or different types of nutrient limitation [6], [10], [11].

In a recent study, Li and colleagues carried out experiments to more deeply assess how the growth of E. coli is affected by carbon (C), nitrogen (N), and phosphorus (P) limitation. They find that E. coli uses three different mechanisms to maintain the same growth rate under the three scenarios [12]. In the same year, Iyer et al. (2018) reported how C and N starvation would affect cell growth and also find distinct patterns of transcription and translation regulation [13].

To qualitatively and quantitatively describe complex experimental findings, a theory or mathematical model is often necessary. One popular approach in E. coli studies [14] is to build a quantitative model that incorporates all relevant biological mechanisms and experimental parameters to describe the complex regulation system, as seen in the work by Hu et al. [15]. On the other hand, a simple model can also be built based on fundamental principles. One benefit of simpler models is that it can often be studied thoroughly to extract insights regarding general behaviors of the biological system, which can then be used in the formulation of more complex models. For example, a class of stochastic model called Totally Asymmetric Simple Exclusion Process (TASEP) which describes particles hopping on a one-dimension chain that exhibits a broad range of complex behaviors, has been used widely in biophysical literature to explore core elements of biological processes that involved transportation of biological matters such as translation [16], [17], [18], [19].

Additionally, Scott and colleagues constructed a simple model by arguing that cell growth must be balanced between varying levels of translation inhibition and varying nutrient quality [4]. Furthermore, Scott et al. also employed the approach of partitioning the ribosomes and protein into multiple factions each with its own function, which is further explored in subsequent studies [20]. Perhaps inspired by this idea, Li et al. constructed a model which partitions ribosomes and mRNA into several classes to describe the translation process under different nutrient limitations.

The model by Li et al. produces predictions that are consistent with experimental results to elucidate on biological mechanisms of cell growth and thus is of interest for further exploration. However, the study focuses on the experimental evidence, so the analytical properties and biological insights of the model are not explored in great detail. For example, it is unclear whether the model dynamics satisfy the assumptions used in its construction, or perhaps contain interesting dynamical properties. Furthermore, the authors use the steady state values for their calculations without guaranteeing their existence, uniqueness, positivity, and stability. In addition, to bypass the need for parameter estimations of two unknowns, the authors absorbed the two unknowns into a new parameter called the saturation parameter and use it to characterize the cell growth under different nutrient profiles. While ingenious from an analytical point of view, this approach does not allow for simulation of the dynamics of the model since the values and ranges of two parameters are unknown. Furthermore, the saturation parameter is somewhat ambiguous and unlikely to be experimentally measurable, reducing its usefulness to provide more in-depth biological insights into mechanisms regulating organismal growth.

Therefore, it is critical to accurately estimate the values of the biologically interpretable unknown parameters to obtain further hypothesis on contrasting effects of limitations by different nutrients on cell growth. In order to do so, we first generalize the model formulation by Li and colleagues and examine the interpretation and possible connection between each parameter. We then calibrate the model using the experimental results from Li et al. and discuss their implications on the growth dynamics and growth strategies under different nutrient limitation conditions. By thoroughly examining the model, we obtain valuable insights that may allow for better future model calibrations and developments to reach a more complete theory of cell growth.

Methods

In this section, we will first motivate and discuss the formulation of the macroscopic model by Li and colleagues and how it was used to describe the experimental data. Then, we will describe our approach to analyze the model and study its connections and implications using the experimental measurements by Li et al. [12].

Model motivation

Various formulations of a mathematical expression for cell growth rate exist [4], [5], [6], [10]. In the work of Li and colleagues, the growth rate was looked at from the perspective of protein mass accumulation in the translation process of each cell. Let be the protein synthesis rate, or the amount of proteins being synthesized per second (amino acids/sec), be the protein mass in a cell (g), and be the average mass of amino acid (g). Then the amount of protein mass produced every second is the product of the protein synthesis rate and the average mass of amino acid (. In other words, the growth rate, μ (), is the relative increase in protein mass rate: where 3600 is the conversion factor from seconds to hour. Additionally, assuming the elongation rate in the translation process is a constant, then the protein synthesis rate (or the amount of amino acids being translated per second) is proportional to the number of ribosomes working on the translation process (). where is the elongation rate. Equating the protein synthesis rate in both equations and solve for the growth rate, we obtained an expression of the growth rate in term of the working ribosomes: Ideally, a model describing the growth of cells (in a chemostat) should produce a uniquely globally stable steady state. Thus, if we suppose the dynamics of cell growth indeed reaches a steady state, then the growth rate is proportional to the constant number of working ribosomes. Therefore, we want a model that can capture the dynamics of different fractions of ribosomes during the translation process to explain the impacts of different nutrient limitation conditions on growth rate.

Model formulation

To provide a theoretical framework for different mechanisms of cell growth under various nutrient limitations, Li and colleagues formulated the following kinetic ODEs describing the translation dynamics of ribosomes (R) and mRNA (M). The assumptions of the model will be marked in bolded texts with the superscript denoting the numerical order.

Ribosomes are separated into three groups. signifies the unbound ribosomes, which includes the free ribosomes and ribosomal sub-units 30S and 50S (i.e. separate, un-assembled sub-units are considered as part of ). is the initiating ribosomes, or the ribosomes that successfully bind to the initiation site and occupy no more than the first 10 codons (since one ribosome occupies about 10 codons on an mRNA). signifies the working ribosomes, which are bound to mRNA starting at the or higher numbered codons.

The bound ribosomes are composed of the initiating ribosomes and the working ribosomes. The differentiation between initiating and working ribosomes is supported using a microscopic model and some experimental data (see the supplementary material in [12]); however, this distinction, especially the length of mRNA corresponding to , has not been explored in depth in literature, so it remains an (1) of the model formulation. We will not discuss the microscopic model in detail here; however, the microscopic version resembles a continuous-differential version of the TASEP, and perhaps can be further analyzed within the existing framework of TASEP [18]. The main functional difference is that during the initiating phase, ribosomes can abort translation, while they cannot do so during the working phase.

The mRNA is categorized as free () and bound (); however, a distinction should be noted because of a novel characterization of the mRNA. is used to indicate when the binding site of the mRNA is available for binding, and is used when the binding site is occupied by a ribosome . In other words, represents the number of mRNA with an available binding site (i.g. unbound at the first 10 codons), and is the number of mRNA whose binding site is currently occupied.

Note that the definition of is not exactly indicative of the starting region that includes the Shine-Dagalno sequence, but rather whether the mRNA is available for another ribosomes to bind to it or not. Li and colleagues argue that the binding region contains up to the first 10 codons of the mRNA (or the approximate occupancy of one ribosome on an mRNA). This classification allows the model to be structurally simple, while capturing key biological features of the process.

General framework

The interactions between different components of the system can be characterized as follows. Unbound ribosomes can bind to free mRNA to become initiating ribosomes. We let represent this binding action. In principle, the more abundant and , the higher should be, so f is a monotone function with respect to both variables. Furthermore, . For mathematical simplicity, f is also assumed to be differentiable everywhere except at . This interaction is represented below, where is the functional form of the binding process. We note that the function form of may contain parameters that relates to the binding rate of to . Such parameters maybe under the influence of nutrient limitations.

Initiating ribosomes can abort translation, which revert to unbound ribosomes and free mRNA. Otherwise, the initiating ribosomes proceed to the 11th codon to become working ribosomes, which also releases the binding site, or increases the number of mRNA with available binding site. It is worth pointing out that the classification of mRNA allows the model to naturally incorporate the possibility that multiple ribosomes can translate a single strand of mRNA at the same time (e.g. polysomes). An alternative method to denote this multiple binding of ribosomes to mRNA is to use a chemical binding formulation, where the mass action approximation would be , where n is the average number ribosomes translating one mRNA strand at a given time and α is a first-order binding rate. However, this process cannot account for the fact that ribosomes can only bind to mRNA one at a time. Another method is to use the TASEP framework as mentioned in the introduction; however, this method is perhaps more suitable if we want to accurately capture the microscopic nature of the translation process.

During the working phase, ribosomes do not abort translation and continue translating at a constant rate until the stop codon. We note that previous studies suggest that ribosomes can abort translation after the initiation phase [21], [22]. However, we elected to ignore this possibility in our model formulation and will discuss the implications of this choice in the discussion section. Upon arriving at the last coding codon, the working ribosomes complete the translation and unbind the mRNA to become free ribosomes (or ribosomal sub-units). In Li et al., the elongation rate is assumed to be independent of the progression rate from initiating to working ribosomes and the translation abortion rate. However, this creates several potential problems such as parameter unidentifiability and model interpretation. Instead, we note that the model readily provides a platform to connect these different transition rates.

We (2) that elongation rate is a function of the nutrient composition of the medium. Since the composition in each experiment is fixed, is taken to be a constant for each experiment. accounts for how many amino acids (aa) are added to the polypeptide chain per second (sec), so it has the unit (). Thus, the amount of initiating ribosomes becoming working ribosomes per second is given by , where is the length of mRNA corresponding to the initiating stage (e.g. is taken to be the first 10 codons in Li et al.). Recall that this is possible because the section of the mRNA that determines its binding availability is not restricted to just the binding site. Similarly, is the fraction of the coding mRNA region () translated by a working ribosome in one second, where is the full length of the mRNA (i.g. from the 5' cap to the 3' end).

Regarding the translation aborting events during the initiating phase, we introduce a function which is the probability of aborting a translation. To simplify the biological mechanisms regarding , such as the regulatory actions of guanosine pentaphosphate (p)ppGpp under nutrient limitation [23], we (3) to be dependent on the nutrient composition of the medium, so it also stays relatively constant within one experiment. Furthermore for consistency with the modeling framework, the aborting event is (4) to be decided at the last coding codon of mRNA corresponding to . It is worth noting that if the initiating region includes up to the first 10 codons (as in the original study), then the aborting event can only happen when the ribosome completely occupies the initiating region (prior to entering the working region) because a single ribosome occupies about 10 codons. Thus, the rate at which becomes is given by , which serves the original purpose of . On the other hand, the rate at which aborts translation is given by , which serves the original purpose of . Together, these processes can be represented as: The differential equation system for this translation model takes the following form.

Note that the model does not contain any degradation or production terms, which is due mainly to the (5) of time scale separation. The processes within translation happen on a much faster time scale (seconds) than the production and degradation of mRNAs and ribosomes (hours). Hence, degradation and production terms can be omitted. Furthermore, if we assume the total amount of ribosomes and the total amount of mRNA are constant for each experiment (perhaps using a quasi-steady state (6) of the transcription process), then and with lead to a simplification of the full system to a two dimensional model. It is worth pointing out that the dependence of and on nutrient composition only accounts for the nutrients specifically reserved or allocated for the translation process. Furthermore, without incorporating the transcription process, we assume the machinery needed for translation is at an optimal constant (). Regarding this observation, Li et al. attempted to study the effect of nutrient up-shift with an extended model; however, without explicitly incorporating the transcription process to account for temporal variation in the total number of ribosomes and mRNA, it may be misleading to extend their results to study the transient and asymptotic dynamics of cell growth after nutrient up-shift.

Li et al. model formulation

In Li et al., the authors (7) takes the form of mass action. This assumption is justified as the number of ribosomes is much larger than the number of mRNA at almost any given time, especially at steady state growth. Furthermore, the effective rate at which ribosomes proceed to the 11th codon to become working ribosomes (), the rate at which initiating ribosomes abort translation (), and the elongation rate () are (8) to be independent. We note a small inaccuracy regarding the term denoting the translation from the last coding codon to the stop codon, which should be the quotient between the elongation rate and the length of mRNA corresponding to working ribosomes (i.g., just the coding sequence on mRNA instead of the entire mRNA). However, since the length corresponding to the working ribosomes is nearly all of mRNA, we will use the original formulation, which does not affect the dynamical properties of the model.

For each nutrient limitation condition, is treated as constant (along with other parameters). Thus, under these assumptions, the model in Li et al. takes the following form.

Using the conservation properties of the system, we can reduce equations (14)–(16) to two dimensional ODEs. A diagram of the system is presented in Fig. 1 and summaries of the model parameters and experimental values are presented in Tables 1 and 2, respectively. Recall the growth rate expression that we introduced in Equation (3). To obtain a more useful form, we first let and be the fraction of rRNA in RNA and the mass of rRNA in ribosomes, respectively. Then the total number of ribosome, , is given by , where ξ is the mass ratio of total RNA to total protein. Define to be the fraction of working ribosomes (i.g. ), then the new growth equation becomes:

Figure 1

The schematics for ribosomal dynamics for Equations (14)–(16) adapted from Fig. 2b in the study by Li and colleagues [12]. Note that the classification for mRNA depends on whether the binding site is available for binding (M) or not (M).

Table 1

Parameter ranges and definitions for E. coli. The ad hoc ranges are based on the estimated values in Li et al. [12]. The ranges indicate the min and max range for each parameter. During the parameter estimations, a more refined range is used for each scenario to obtain the best convergence. The model dynamic variables are listed at the end. The experimental variables and parameters are given in Table 2. The parameter indicated with ⁎ represents similar quantities with the parameters in the parentheses from Loladze and Elser [31].

Symbol	Definition	Ranges	Ref
kf	effective rate constant for unbound ribosomes and free mRNA to initiate translation (1/sec)	[1.0e − 5,2.4e − 4]	ad hoc
kr	rate constant for initiating ribosomes to abort translation (1/sec)	[5e − 3,5e − 2]	ad hoc
kel⁎(σ_a)	elongation rate (aa/sec)	[6,20]	ad hoc
kp	rate constant for initiating ribosomes to transition into working ribosomes (1/sec)	[0.3,1.4]	ad hoc
fm	fractional mass of mRNA	0.02	[24]
Vc	cell volume (m3)	10−18	[25]
Cp	concentration of proteins (g/m3)	2.4 × 105	[26]
mr	mass of rRNA component in a ribosome (g)	2.8 × 10−18	[27]
mnuc⁎(m_r)	average mass of a nucleotide in RNA (g)	5.4 × 10−22	[28]
maa⁎(m_a)	average mass of an amino acid in protein (g)	1.8 × 10−22	[29]
Maa	average length of the initiating region (aa)	10	[30]
Naa	average length of mRNA (aa)	300	[30]
na	N content averaged over 20 amino acids	17%	[31]
nr	N content averaged over 4 nucleotides	15%	[31]
pr	P content averaged over 4 nucleotides	9%	[31]
la	length of RNA polymerase (aa)	3400	[31]
lr	length of rRNA (nt)	4560	[31]
σr	transcription elongation rate (nt﹨sec)	71	[31]
ϕr	fraction of total protein that is RNA polymerase actively transcribing rRNA	0.2%	[31]
Pm	total protein mass in a cell (g)	CpVc	[12]
Rt	total number of ribosomes per cell	P_mξf_rm_r	[12]
Mt	total number of mRNAs per cell	P_mξf_mN_aam_nuc3	[12]
Ru	Number of Unbound ribosomes (count)
Ri	Number of initiating ribosomes (count)
Rw	Number of working ribosomes (count)
Mf	Number of mRNA with available binding site (count)
Mb	Number of mRNA with bound binding site (count)

Table 2

Summary of the mean of the processed experimental values for wild type E. coli from Li et al. [12]. The values of f are used as input and the values of k are used to cross check with our independent estimation of k. The parameter indicated with ⁎ represents similar quantities with the parameters in the parentheses from Loladze and Elser [31].

Symbol	Definition	C-limit	N-limit	P-limit
ξ	Mass ratio of total RNA to total protein	0.16	0.16	0.08
μ	steady-state growth rate (h−1)	0.09	0.09	0.09
Φ_R_u	Fraction of ribosomes not bound to mRNA	0.69	0.49	0.35
Φ_R_i	Fraction of ribosomes located at the first 10 codons of mRNA that may abort translation prematurely	0.03	0.03	0.03
ΦR_w⁎(ϕ_a)	Fraction of ribosomes located after the first 10 codons of mRNA, which contributes to the protein production	0.29	0.47	0.62
fr	Fractional mass of rRNA among total RNA	0.64	0.59	0.51
kel	Elongation rate (aa/sec)	12.5	7.5	12.5

Li and colleagues carefully parametrize the models with global constants and their experimental measurements [12]. In addition to the estimable parameters, the values of and are unknown. Instead, a new parameter called the saturation parameter (S) was introduced. As S is the only unknown, the value of S can be approximated. Finally, S was used to characterize the difference of cell growth under different nutrient limitation conditions.

Regarding cell growth under various nutrient limitations, there are three main findings by Li and colleagues corresponding with three different mechanisms under C-, N-, or P-limitation. Their experimental results provide estimates of the fractions of free, initiating, and working ribosomes, i.g. respectively, the elongation rates, and the abundance of ribosomes () calculated via the mass ratio of total RNA to total protein.

Under C-limitation, there is a low fraction of working ribosomes, which results in lower translation rate. However, the pool of ribosomes and the elongation rate are both high.

Under N-limitation, the elongation rate is slowed (due to stalling at glutamine codons) leading to a slow translation. However, the fraction of working ribosomes is moderate and a large pool of ribosomes compensates for the slow elongation.

Under P-limitation, there is a low abundance of ribosomes, which slows down translation. However, the lack of ribosomes is made up for by the high fraction of working ribosomes and fast elongation rate.

Parameter estimation

The main objective of our parameter estimation is to estimate possible values of and using the experimental values for , (at steady state), and ξ in each condition (C-limited, N-limited, and P-limited). Table 2 contains the experimental values from Li et al. [12]. Table 1 contains the ranges and fixed values of all parameters. Since one of our goals is to quantify/verify the difference in elongation rate under different nutrient-limited conditions, we will also estimate , which means we will need to re-estimate (since is calculated at steady-state based on ).

Data: Li and colleagues performed various experiments using batch and continuous (chemostat) cultures, where the experiments using batch cultures were mainly for supporting findings from the continuous cultures. For our purpose, we will use the three sets of data for the E. coli strain NCM3722 at growth rate of 0.09 under three different nutrient limitations. In each experiment involving the wild type E. coli, Li and colleagues determine the mass ratio of total RNA to total protein and quantify rRNA and ribosome fractions. The elongation rate was also measured using lacZ induction; since the values for are from Fig. 6.a in the supplementary materials of Li et al., we approximate their values from that figure. The means of the processed experimental data are summarized in Table 2. Note that, in our parameter estimation, we will try to estimate the elongation rate directly and independently, so we will use their experimentally measured value to cross check. The values for are used as input for the model. For specific details of the experiments, please refer to the method section in Li et al. [12].

Our approach is similar to that of an Approximate Bayesian computation scheme, for additional details see [32]. Our goal is not to obtain a precise estimation of the posterior distribution, but rather to find admissible estimations that satisfy all biological measurements at the same time. Let and . To determine the range of in step 1, we started with the baseline estimates from Li et al. whenever possible. If the range of certain parameters was not available, we were able to obtain an initial estimate using the upper and lower bounds of other parameters in combination with experimental estimated values of x at steady state. We slowly expand the range, about 10–20% each time, within biologically reasonable values until the fitting found no acceptable values near the lower and upper bounds.

First, we created samples from some uniform distributions .

Then, we simulated data set corresponding to .

Finally, we compared with the experimental data from Li et al. [12]. In this case, we wanted each component of to be within 10% error of the corresponding element in . Thus, we accepted the proposed parameter sets if where representing , respectively, and is a measure of the absolute distance between its arguments.

To test for the possibility that a completely different set of parameters might also satisfy the acceptable criteria, we also expanded the lower and upper bounds by many-fold (10 to 100) with the exception of parameters that were well-constrained biologically (e.g., should not be too close to 0 or much larger than 20, etc.). While more systematic approaches exist in literature, the stable convergence in our results demonstrates the sufficiency of our method in estimating the parameter distributions.

Global sensitivity analysis

To study the influence of parameters under each nutrient profile and in general, we calculate the partial rank correlation coefficients (PRCC) using 2000 samples obtained from the Latin hypercube Sampling. The details of the scheme can be found in the studies by [33], [34]. We vary the parameters within 1% range assuming uniform distribution. As a separate scenario, we also vary the parameters within a 25% range to study the impact of more uncertainty in the parameter values. The sensitivity and uncertainty analyses are carried out for each nutrient limitation using the estimated and fixed parameters.

Results

Mathematical analysis

Recall that in the study by Li and colleagues, the focus was to establish a mathematical model that can explain the experimental observations. Thus, certain assumptions were made during the model formulation; however, whether the dynamical behaviors of the model satisfy those assumptions was not considered. Furthermore, the existence of a unique globally asymptotically stable positive steady state was assumed and used as part of several key algebraic manipulations. For these reasons, we aimed to verify the basic biological assumptions and the asymptotic behaviors of the model.

Theorem 1, Theorem 2 establish that each variable is always positive and bounded by the biological bounds ( and ). Theorem 3 shows that the system always exhibits a unique biologically relevant positive steady state. Finally, Theorem 4 shows the local asymptotic stability of the positive steady state, which is indeed globally asymptotically stable by Theorem 5. Fig. 2 shows an example that demonstrates the stability of the system. These stability results provide mathematical justification for the algebraic manipulations that were done in Li et al. [12]. The proofs are given in subsection A.1 of the Appendix A.

Figure 2

An example to demonstrate that the unique positive steady state is globally asymptotically stable. (a) Phase-plane with trajectory – note that the positive steady state exists within the biological bounds of R. (b) Effects of varying model parameters on the asymptotic behaviors of the model. While the stability does not change, the relative ratios between and can be altered by changing k or k. This is consistent with the steady state expression from Theorem 3, .

Positivity

The system is positive given positive initials.

Boundedness

Each compartment of ribosome is bounded by and each compartment of mRNA is bounded by .

We note that the positivity and boundedness of the general model also hold similarly.

Positive equilibrium

With the biological bounds of and , then there is a unique positive equilibrium at , where and with , and

Local stability

The biological equilibrium at is locally asymptotically stable.

Global stability

The biological equilibrium at is globally asymptotically stable in the region bounded by

Computational results

Recall that our parameter estimation process focuses on finding sets of parameters that allow the model to fully capture the experimental values within an acceptable error range. Table 3 gives the estimated range of the parameters (, , , ). Important simulation results are shown in Figs. 3 and 4. Multiple rounds of simulations confirmed that there were no significant differences in the convergence of the estimated parameters, see Figure 7, Figure 8 in Appendix subsection A.2

Table 3

Parameter estimations for the purpose to provide consistency between the model prediction and experimental results. The tables show the mean and standard deviation of the parameter values. The median and 25–75% quartiles are shown in Fig. 4 with the corresponding fitting in Fig. 3.

Params	C-lim	C-sd	N-lim	N-sd	P-lim	P-sd
kf	1.76 × 10−5	4.14 × 10−6	2.51 × 10−5	5.15 × 10−6	1.54 × 10−4	3.39 × 10−5
kr	2.77 × 10−2	1.29 × 10−2	2.77 × 10−2	1.29 × 10−2	2.79 × 10−4	1.30 × 10−2
kel	1.42 × 101	1.97 × 100	9.28 × 100	1.20 × 100	1.62 × 101	1.90 × 100
kp	4.59 × 10−1	6.43 × 10−2	4.94 × 10−1	6.23 × 10−2	1.12 × 100	1.40 × 10−1

Figure 3

Model vs. experiments – the red star is the estimated value from the experiments. The red line is the median of the parameter estimation. The box represents the 25%–75% quantiles. The variables used for fitting are for: (a) μ, (b) , (c) , (d) , (e) ξ. The stars staying within the box means the simulation results agree well with the experimental results. The resulting parameter estimates are given in Table 3 and Fig. 4.

Figure 4

Parameter estimation from data fitting in Fig. 3 – the red star is the estimated value from the experiments (for the elongation rate only). The red line is the median of the parameter estimation. The box represents the 25%–75% quantiles. The red crosses are the outliers. The estimated parameters are: (a) k, (b) k, (c) k, (d) k.

Figure 7

Testing the simulation convergence in Fig. 3 under four different scenarios: (a) R = 0.9R,R = 0.05R,R = 0.05R; (b) R = 0.8R,R = 0.10R,R = 0.10R; (c) R = 0,R = R,R = 0.0; (d) R = 0,R = 0,R = R. The convergence results are stable and agree with the result in Fig. 3 in all cases.

Figure 8

Testing the simulation convergence in Fig. 4 under four different scenarios: (a) R = 0.9R,R = 0.05R,R = 0.05R; (b) R = 0.8R,R = 0.10R,R = 0.10R; (c) R = 0,R = R,R = 0.0; (d) R = 0,R = 0,R = R. The convergence results are stable and agree with the result in Fig. 4 in all cases.

Using the estimated parameters in each case, we simulate the dynamics of ribosomes under C-, N-, and P-limitation, see Fig. 5. For the simulation, we assume that initially all ribosomes are unbounded, meaning . We observe that in all three scenarios, the dynamics of and are monotonically decreasing and increasing, respectively. reaches an early peak within the first few hours (1.71 hr under C-limitation, 5.14 hr under N-limitation, and 2.13 hr under P-limitation) before settling into a steady state. The steady state of is reached (within 10%) at 24.75 hr under C-limitation, 41.61 hr under N-limitation, and 22.17 hr under P-limitation. Additionally, under C-limitation, the ribosomal dynamics take 24.75 hr and 37.59 hr to reach within 10% of the equilibrium for and , respectively. Under N-limitation, the ribosomal dynamics take 41.61 hr and 42.91 hr to reach within 10% of the equilibrium for and , respectively. And under P-limitation, the ribosomal dynamics take 21.56 hr and 17.40 hr to reach within 10% of the equilibrium for and , respectively.

Figure 5

Ribosomes dynamics under different nutrient limitation conditions – simulated using the mean values from the parameter estimations presented in Table 3. The simulations are for: (a) C-limitation, (b) N-limitation, (c) P-limitation.

We note that, while it is non-optimal to use the same set of initial conditions for the comparisons in Fig. 5, our main goal was to compare the rate at which the ribosome profile settles to the steady states under each nutrient limitation. Thus, by starting with the same initial values, the combined effect of all parameters dictates the reaction rate, allowing us to make relative comparisons across treatment. We also tested the ribosomal dynamics under four different scenarios, while there are qualitative differences, the ribosomal dynamics under C- and P-limitation is still faster than under N-limitation, see Fig. 9 in Appendix subsection A.2.

Figure 9

Testing the ribosomal dynamics, as compared to Fig. 5, with four different initial conditions: (a) R = 0.9R,R = 0.05R,R = 0.05R; (b) R = 0.8R,R = 0.10R,R = 0.10R; (c) R = 0,R = R,R = 0.0; (d) R = 0,R = 0,R = R. While there are qualitative differences in the dynamical behaviors (especially in the third and fourth scenarios), under C- and P-limitation, the convergence rate is faster compared to the situation under N-limitation.

Fig. 6 shows the global sensitivity results with respect to μ for large uncertainty in the parameters (10–200%). The result is with respect to the growth rate under C limitation; however, similar results are obtained for N- and P-limited conditions. We also carried out sensitivity analysis when the uncertainty is small (1%). Under both assumptions about the uncertainty level, we observe that μ is sensitive to all parameters with the exception of (all results are statistically significant except for ). Moreover, we note that seems to be the most impactful parameter. The sensitivity analysis is consistent with what we observe in our parameter estimation, where the values and ranges of are similar across all three nutrient limitations, see Fig. 4. This makes sense from our generalized derivations, where there is a hidden relationship between , and . This suggests that least one of the three parameters is not identifiable, a recurring problem in mathematical modeling of biological systems [35].

Figure 6

Global sensitivity with respect to the growth rate μ given large uncertainty in parameter values. The results indicate that k is the most sensitive parameter, while k is the least sensitive parameter. Additionally, the elongation rate k is among the most sensitive parameters, but fails behind k. We suspect this is an artifact of the model formulation where k is assumed to be independent of k.

Connection to the growth rate hypothesis

The growth rate hypothesis relates the abundance of P-rich ribosomal RNA to differences in maximal growth rate. But when the maximal growth rate is similar under different nutrient limitation conditions, it is of interest to study the N:P ratio of cell under each condition. Switching from considering Protein:RNA to N:P helps to identify the stoichiometric mass balance aspect of growth-limitation and gives a first-principle mechanistic understanding as to how organisms alter their elemental allocation to optimize growth. It also gives us simple predictions that can be tested empirically, so we can work iteratively to better understand nutritional controls on organismal biology.

Using the framework established by Loladze and Elser [31], we shall estimate the N:P ratio of the E.coli in the three nutrient limitation conditions in the study by Li and colleagues. Let θ denote the N:P mass ratio of the system. Note that the N content within a cell can be approximated by (N content average over 20 amino acids) × (mass of all protein) + (N content average over 4 nucleotides) × (mass of all rRNA). And the P content within a cell can be approximated by (P content average over 4 nucleotides) × (mass of all rRNA). Recall that is the mass of all proteins inside the cells. Let r be the mass of all rRNA; then the N:P mass ratio can be expressed as where and represent the N content average over 20 amino acids and 4 nucleotides, respectively, while represents the P content average over 4 nucleotides. Loladze and Elser estimated that is between 15% to 18%, is approximately 15%, and is about 9%, so we take , , and to be 17%, 15%, and 9%, respectively, for our calculations. By using a model that couples protein and rRNA synthesis, Loladze and Elser reformed the above expression into the equivalence: The definitions of the parameters are recorded in Tables 1 and 2. Here, we note that and are the same as and in the Li et al. model. Furthermore, the and also represent similar quantities in both models. The bridge between the two models is the unique, globally, asymptotically stable equilibrium that allows us to pass the estimations from one model to the other. Thus, we can first estimate the values of and with the Li et al. model, then use them to estimate the N:P mass ratio in each limitation condition, which then gives us the actual N:P ratio via a conversion based on their atomic numbers. The results are presented in Table 4.

Table 4

Estimation of N:P ratio under each nutrient limitation condition. The values are converted from θ using atomic ratio of N:P (14.007u:30.97376u). The first three estimates (without ⁎) used the experimentally reported values in Li et al. The second three estimates (with ⁎) used the numerically fitted values.

Condition	C-lim	N-lim	P-lim	C-lim⁎	N-lim⁎	P-lim⁎
N:P	9.76	9.68	12.57	10.16⁎	10.35⁎	13.80⁎

Discussion

The effects of nutrient conditions on cell growth have attracted significant modeling effort due to their wide implications, such as in the fields of ecology [36], [37], [38], [39], [40] and agriculture [41], [42], [43], [44], [45]. Furthermore, these effects introduce interesting dynamical properties extending the Lotka-Volterra framework [36], [37], [38], [46], [47]. While the overall aim is to reach a deeper understanding of biological systems by means of constructing a quantitative framework, the mathematical and biological properties of the model are often neglected. Without careful consideration, models may violate fundamental principles, such as mass conservation, or give unrealistic implications about the biological systems that they are meant to describe. In this study, we analyzed the model that Li et al. [12] formulated to describe their experimental data. Our main findings are summarized below.

A general framework of Mathematical models can be used to capture the essences of highly complex and nonlinear dynamic processes in a biological system. However, to do so, certain assumptions and generalization must be made. Thus, building upon the model formulation by Li and colleagues, we created a simple and generalizable mathematical framework for the translation process of cells such as E. coli. We stated all major assumptions related to the model formulation and pointed out possible limitations, especially with regards to the lack of the transcription process. Furthermore, when certain parameters are kept independent of one another, our theoretical and simulation results suggest a possible issue with parameter identification and interpretation.

Existence, uniqueness, and stability of the positive steady state. We analyzed the mathematical properties of the model developed by Li and colleagues to demonstrate that it satisfies the basic biological assumptions such as positive invariance and boundedness, specifically by the biological bounds established by the total number of ribosomes () and the total number of mRNA () from the model formulation. Furthermore, we confirm the existence of a unique positive steady state under all possible ranges of the parameter values, which is also globally asymptotically stable. This validates the algebraic manipulations based on the steady state values in Li et al. [12].

Model consistency with experimental results. Our results show that the model can produce simulation results that agree well with experimental results under all three scenarios of C-, N-, and P-limitations. This supports the usefulness of the framework in aiding our understanding of the nonlinear relationships between growth, biochemical and stoichiometric under differential nutrient limitation. Furthermore, we were able to estimate the values of the effective rate constant for unbound ribosomes and free mRNA to initiate translation () and the rate constant for initiating ribosomes to abort translation (), which are previously unknown due to grouping into the saturation parameter S. While stays relatively constant across three nutrient regimes (perhaps due to identifiability issue), we observe significant variation in the estimated values of between P-limited vs. C- and N-limited environments. Even when transcription and energetics are not explicitly incorporated, the model is able to capture the observed experimental differences under three different nutrient profiles.

Under P-limitation, the high fraction of working ribosomes and fast elongation rate are not sufficient to account for the lack of total ribosomes. Assuming the rate that unbound ribosomes and free mRNA initiate translation, , is similar across three nutrient conditions, we were unable to obtain simulations consistent with the experimental results in Li et al. Instead, by letting vary within an extended range (Table 1), we were able to simulate the experimental results. Specifically, we find that, under P-limitation, is several times larger than under C- or N-limitation, see Table 3. This implies that, other than having a higher fraction of working ribosomes and fast elongation rate, P-limited cells may also require a fast transition rate from initiating to working ribosomes to achieve the optimal growth rate. Additionally, our estimate of (transition rate from initiating ribosomes to working ribosomes) is also significantly higher under P-limitation, while remaining relatively similar under C- and N-limitation, see Fig. 4(c). These estimations suggest that, in order to overcome the lack of total ribosomes, multiple processes must work more rapidly. We suspect that this finding is mainly a consequence of lower value of ξ under P-limitation (see Table 2), which leads to a much lower amount of total ribosomes. For instance, a possible explanatory hypothesis is that suppose a similar amount of energy is reserved for translation across three different nutrient limitations. When there are less working components, it may appear that each component is working faster or more efficient due to a higher per capita energy allocation.

The transition rate from initiating to working ribosomes is the most important factor regardless of which nutrient is limiting. The results of our sensitivity analysis demonstrate that all parameters are sensitive with respect to the growth rate, with the most sensitive parameter being and the least sensitive parameter being . This partially explains why there is not much difference between the estimated values of across nutrient profiles. On the other hand, as pointed out in the model formulation, parameter unidentifiability maybe an issue for , and . Moreover, the difference in sensitivity becomes clearer under higher uncertainty, but the overall sensitivity order remains the same. This is surprising since the elongation rate, which directly impacts the growth rate, is intuitively thought to be the most important parameter that governs the process. A possible explanation is that the impact of is amplified artificially due to the assumption of independence between and . In our general formulation, since is proportional to , the sensitivity of should increase. Alternatively, if we consider translation as a transport problem (e.g., cars moving along a single highway), then represents the high traffic point (or the bottle neck) of the system. Hence, controls the flow of the system. This effect is further amplified by the assumption that there is no aborting events after initiation. A simple way to address this is to simply allow to be a function of . In other words, we can replace with a function of and a predetermined length of the initiating region on the mRNA strand. In addition, we may consider having a low but equal aborting rate after initiation. By doing so, we expect to be the most significant parameter in the general framework.

Ribosomal dynamics under C- and P-limitation occur on a faster time scale than under of N-limitation. Nitrogen contributes in many aspects of cell growth, such as in the production of proteins, DNA, and RNA. Using the values obtained from our estimations, we simulated the ribosomal dynamics under C-, N-, and P-limitation, see Fig. 5. We found that different types of nutrient limitation not only require cells to use different coping mechanisms, but also change the time scale of cell response. Under C- and P-limitation, we observed ribosomes reaching close to equilibrium (within 10%) significantly faster than under N-limitation for all three factions of the ribosomes. This may be due to the slow elongation rate under N-limitation, but other mechanisms may also play a role. We speculate that our observation is due to the involvement of N in forming macromolecules that are necessary in cell growth. Alternatively, since carbon is also present in these macromolecules in addition to energy pathways, and ATP concentrations in N-limited organisms are also high, our observation may also be an effect of energy limitation in the C- and P-limited treatments in that they reach equilibrium faster due to coupled energy and elemental limitation for ribosome synthesis. Unfortunately, we would need a higher resolution model, particularly with transcription, to study this hypothesis in more detail.

Connection to the growth rate hypothesis. We connected the Li et al. model with the framework by Loladze and Elser [31] to estimate the N:P ratio under the three nutrient limitation conditions. This is possible because both models share similar underlying processes and are intrinsically globally stable. Our results (Table 4) show consistent findings between experimental estimates and numerical estimates of the N:P ratios. We observed that the N:P ratios under P-limited are higher than under N-limited conditions at the same growth rate, which is expected. However, the Li et al. model does not contain the description of the transcription process, so we used the estimates from the study by Loladze and Elser for the transcription parameters (fraction of actively transcribing rRNA, , and the transcription rate, ). This may explain the similar N:P ratios under C-limited and N-limited conditions. Furthermore, maybe smaller under P-limited condition, which would lead to a higher N:P ratio.

Comparison of estimated parameters with literature

The translation abortion rate (): previous studies suggest that it is possible for ribosomes to dissociate at equal rates at every codon (even after the initiating phase). By using an exponential decay model, the drop-off rate was estimated to be on the order of to per elongation step in stress-free environments [21], [22]. Our estimated range for falls within this range and further shows that aborting events occur at a relatively constant rate for similar growth rates. However, it is difficult to accurately isolate the aborting rate using a simple exponential decay model because the exponential drop-off observed in the ribosome density along the mRNA is due to multiple factors. Thus, a single-parameter model cannot distinguish them. If data permit, a direct extension to the model would be to examine the possibility of aborting events taking place after initiation.

The translation elongation rate (): our estimated values for the translation elongation rates agree with the experimentally derived values from the original study by Li and colleagues. Furthermore, our estimation also falls within existing estimated range for elongation rate 6–20 aa/sec [20], [48], [49], [50]. When comparing to similar experimental settings, our estimates also agree with the estimates by Iyer et al. for C and N starvation E. coli of about 12.9 aa/sec [13]. Additionally, the hopping rate is often set around 10 aa/sec in TASEP models. [19].

The translation initiation rate (): previous studies using the TASEP modeling framework and experimental data for E. coli growth estimated the ribosome-dependent initiation rate to be around 0.09–0.15 per second [17], [19]. This value is equivalent to in our model, which works out to be approximately 0.11, 0.097, 0.19 per second under C-, N-, P-limited conditions, respectively, at equilibrium. This means our estimates for the initiation rates are within range of previous estimates.

The transition rate into working ribosomes (): is a unique parameter specific to our model, so it is difficult to find comparable estimates. On the other hand, we note that it is likely to depend on the length of the initiating of mRNA strand, the elongation rate, and the aborting rate. In future iteration, it is advisable to replace with a function of the elongation rate to eliminate this uncertainty.

There is significant value in the ability to predict cell growth based mainly on information about nutrient supply conditions. Being able to do so would allow us to naturally link supply to demand to biological function. For example, if we can make predictions about what the C:N:P ratios should be across different nutrient limitation treatments for a diverse group of organisms accurately, we can theoretically determine the growth rates off of elemental measurements of any animals out of the field. As there is currently no reliable way to carry out such task, it is difficult to extend these fundamental relationships governing growth on organismal scales to population and ecosystem-level dynamics. To this end, mathematical models are useful tools to facilitate this possibility. However, to do so, models need to be able to capture the expected qualitative behaviors and satisfy the standard biological constraints of the system. In this study, we provide a generalized framework for the translation process and study the mathematical and biological properties of the model. For simplicity, our modeling approach must make certain assumptions such as the independence between the translation and transcription process. Furthermore, transcription and translation processes are energy intensive, linked by the molecule Adenosine Triphosphate (ATP). Cells require energy to synthesize protein and DNA, which occur when energy is available via ATP reactions that themselves use phosphorus [51]. Thus, direct extension of the model may include explicit incorporation of the nutrient supply conditions and the coupling of transcription and translation processes with energy [52].

With regards to an extension to incorporate energy, the study by Scott et al. [4] inspires a simple possibility. Since cells must partition their energy usages into multiple processes, we can consider several compartments, each related to a major process such as transcription and translation elongation rates. For example, a compartment can be used to keep track of the amount of energy used for translation. Then, must influence the value of , so we can introduce a new variable for the translation elongation rate that is a function of and , perhaps . At equilibrium, we can argue that the allocation strategy is optimal, so an equilibrium must be reached, which should also be proven mathematically. This allows for a simplification of all the energy compartments into constants that represent their relative fraction of the total energy at equilibrium. Then the effect of the energy allocation strategy can be studied by estimating the values of these constants using energy data. If the optimal condition is loosened for generality, we would instead arrive at a much more complex problem of resource allocation, where cells are able to adjust their resource allocation continuously with respect to the differential changes in the environment for optimal growth. On the other hand, such system may be able to accurately reveal the different strategies that cells use in response to different growth condition.

Mathematically, certain interaction forms may provide interesting nonlinear dynamical behaviors for the model that contains biological insights [53]. For example, a more biologically realistic functional form of may result in multiple equilibria within the biologically bounded region, or none at all. This may have significant implications for the biological system and help deepen our qualitative and quantitative understanding of the growth process.

Declarations

Author contribution statement

Tin Phan: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Changhan He: Performed the experiments; Analyzed and interpreted the data; Wrote the paper.

Irakli Loladze: Conceived and designed the experiments; Analyzed and interpreted the data; Wrote the paper.

Clay Prater: Analyzed and interpreted the data; Wrote the paper.

Jim Elser, Yang Kuang: Conceived and designed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Funding statement

YK, CH, IL and TP were supported by and [DEB-1930728 and R01GM131405-02]. JJE was supported by [DEB-1930816] and CP was supported by [DEB-1930736].

Data availability statement

Data included in article/supplementary material/referenced in article.

Declaration of interests statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.