Multiscale modelling of bioprocess dynamics and cellular growth

Abstract

Background

Fermentation processes are essential for the production of small molecules, heterologous proteins and other commercially important products. Traditional bioprocess optimisation relies on phenomenological models that focus on macroscale variables like biomass growth and protein yield. However, these models often fail to consider the crucial link between macroscale dynamics and the intracellular activities that drive metabolism and proteins synthesis.

Results

We introduce a multiscale model that not only captures batch bioreactor dynamics but also incorporates a coarse-grained approach to key intracellular processes such as gene expression, ribosome allocation and growth. Our model accurately fits biomass and substrate data across various Escherichia coli strains, effectively predicts acetate dynamics and evaluates the expression of heterologous proteins. By integrating construct-specific parameters like promoter strength and ribosomal binding sites, our model reveals several key interdependencies between gene expression parameters and outputs such as protein yield and acetate secretion.

Conclusions

This study presents a computational model that, with proper parameterisation, allows for the in silico analysis of genetic constructs. The focus is on combinations of ribosomal binding site (RBS) strength and promoters, assessing their impact on production. In this way, the ability to predict bioreactor dynamics from these genetic constructs can pave the way for more efficient design and optimisation of microbial fermentation processes, enhancing the production of valuable bioproducts.

Background

There exists a constant challenge in developing efficient biological systems that use microorganisms such as bacteria in fermentation processes to achieve higher yields [1]. Researchers have directed their efforts towards the engineering of synthetic strains that exhibit superior efficiency, accelerated growth and greater resistance to external changes. However, the creation of new strains often requires substantial resources and time, often resulting in unfavourable outcomes [2, 3].

In response to these challenges, mathematical models have been developed to predict cellular behaviour and maximise the productivity of target metabolites [2, 4]. These models range from simple methodologies, such as the Monod kinetic mass balance [5], to more sophisticated approaches. One example is metabolic flux analysis (FBA), which quantitatively assesses cellular metabolic networks [6]. These models allow for a thorough evaluation of cellular physiology and are used to optimise and design new strains that meet specific objectives [7,8,9].

Despite the abundance of available models, they have limitations. Most face accuracy issues due to the metabolic burden caused by genetic overexpression [10]. This complicates the design of synthetic strains aimed at maximising the production of desired metabolites and cell growth.

The task becomes even more challenging when considering competition for limited cellular resources [11, 12]. These resources are optimally distributed to promote cell growth [13]. However, in recombinant cells expressing heterologous proteins, the majority of resources are diverted towards the expression of recombinant genes. This increases the production of these proteins but reduces the resources available for the expression of other genes responsible for the cell’s metabolic activities. As a result, cellular growth is directly affected [14, 15].

In order to gain a deeper understanding of the competition for limited cellular resources, gene expression and their relationship in fermentation processes, we have developed a multiscale model that effectively integrates cell growth, substrate consumption and allocation of limited resources with the dynamics of a bioreactor. Our goal is to overcome the limitations of current models by providing a more detailed and comprehensive understanding of cellular physiology and its interaction with the operational conditions of a bioreactor.

For its development, we expanded the whole-cell mechanistic model proposed by Weisse et al. [11], specifically designed for Escherichia coli. This model links transcription and translation with the allocation of cellular energy, ribosomes, and the proteome. Due to its characteristics, it has been employed in various studies to predict the burden in genetic circuits [14, 16]. The model expansion involved incorporating variables related to acetate and recombinant protein production, followed by its integration with the differential equations describing the dynamics of a bioreactor. We considered the use of a batch bioreactor and validated the model using data reported in the literature for different Escherichia coli strains [17]. Subsequently, we conducted simulations to analyse the behaviour of the volumetric productivity of biomass and recombinant protein in response to modifications in the internal model parameters, such as the transcription rate of the recombinant protein and the strength of the ribosomal binding site (RBS). In this way, we identified design strategies that maximise the productivity of the recombinant protein by optimising the aforementioned parameters. Finally, we evaluated changes in acetate productivity in response to variations in the induction of the heterologous gene and the substrate transporter enzyme, finding that acetate production is directly related to the cellular stress response caused by the overexpression of heterologous proteins.

Our main objective is for the multiscale model to be a robust tool that paves the way towards a more detailed understanding of cellular functioning in fermentation processes. We anticipate that this advancement will foster the development of new strategies for the design and execution of bioprocesses.

Results

Mechanistic model of cell growth

For the development of the multiscale model, we initially adapted a mechanistic model of cell growth developed by Weisse et al. [11], which is based on the allocation of limited resources, gene expression and cell growth. The mechanistic model describes the presence of four essential proteins that facilitate the vital functions of the cell corresponding to substrate transport enzymes $e_t$, maintenance enzymes q, metabolic enzymes $e_m$ and the protein part of the ribosomes r.

Each protein incorporates two additional molecules corresponding to $\text {mRNA}{x}$, called $m_x$ and a complex formed by the binding of the ribosome and $\text {mRNA}{x}$, called $c_x$. Additionally, we assume that no protein degradation occurs and the specific growth rate $\lambda$ is a function of the energy levels a, the ribosome complex $c_x$ and the translational elongation rate $\gamma$. Figure 1A provides a more detailed understanding of the functioning of the molecules and the ordinary differential equations presented below determine the dynamic behaviour of the mentioned variables.

$$\begin{aligned} \frac{d r}{dt}&= v_r(c_r,a) + \sum _{x \in \{r,e_t,e_m,q\}}(v_x(c_x,a)-k_b r m_x +k_u c_x) -\lambda r, \end{aligned}$$

(1)

$$\begin{aligned} \frac{d y}{dt}&= v_y(y,a) - \lambda y, \quad y\in \{e_t, e_m, q \}, \end{aligned}$$

(2)

$$\begin{aligned} \frac{d m_x}{dt}&= w_x(a) + v_x(c_x,a)-k_b r m_x +k_u c_x -(\lambda + d_m)m_x, \end{aligned}$$

(3)

$$\begin{aligned} \frac{d c_x}{dt}&= k_b r m_x - k_u c_x v_x(c_x,a)-\lambda c_x, \quad x\in \{e_t, e_m, q, r\}, \end{aligned}$$

(4)

$$\begin{aligned} \lambda&= \frac{\gamma (a)}{M} \sum c_x, \ \ \ \text {with} \quad \gamma =\frac{\gamma _{max} a}{K_{\gamma }+a}. \end{aligned}$$

(5)

where $d_m$ represents the mRNA-degradation rate, $k_b$ the mRNA-ribosome binding rate, $k_u$ the mRNA-ribosome unbinding rate, $\gamma _{max}$ the maximum translational elongation rate, M the total cell mass and $K_{\gamma }$ the translation elongation threshold. The translation rate $v_x$ and transcription rate $w_x$ follow the following expression:

$$\begin{aligned} v_x&= \frac{c_x \gamma _{max} a}{n_x (a+K_{\gamma })}, \end{aligned}$$

(6)

$$\begin{aligned} w_x&= \frac{w_{x_{max}} a}{\theta _x + a}, \end{aligned}$$

(7)

where $n_x$ corresponds to the amino acid length, $w_{x_{max}}$ is the maximum transcription rate and $\theta _x$ is the transcription threshold. The transcription rate applies to all enzymes in the cell, except for the maintenance enzymes q. These maintenance enzyme are self-regulated, allowing them to maintain stable levels despite varying cell growth rates [11] and complies with the following equation:

$$\begin{aligned} w_q=\frac{w_q a}{\theta _x + a} \frac{1}{1+(q/K_q)^{h_q}}, \end{aligned}$$

(8)

where $K_q$ corresponds to the autoinhibition threshold constant and $h_q$ corresponds to the Hill coefficient of autoinhibition.

Expansion of the mechanistic cell growth model

Based on this mechanistic framework, we have extended the original model by incorporating new variables and parameters that specifically address the production of a recombinant protein as a product of interest and the synthesis of acetate as an unwanted byproduct.

To expand the mechanistic model of cell growth, we incorporated a distinction in the metabolic enzymes $e_{m}$ produced by the cell, classifying them into two types: the fermentation metabolic enzymes $e_{m_f}$, which participate in the fermentative pathway and are responsible for the direct production of acetate and the respiration metabolic enzymes $e_{m_r}$, in charge of the metabolic activities associated with cellular respiration. This differentiation allows for a more detailed and functional description of the cellular processes involved in energy production [18, 19]. Thus, in our model, the metabolic enzymes are defined as follows:

$$\begin{aligned} e_m = e_{m_f} + e_{m_r}. \end{aligned}$$

(9)

As shown in Fig. 1A, the external substrate is imported by the cell and transformed into internal substrate, which is then metabolised into energy through the action of metabolic enzymes. According to the principles of the mechanistic model, the amount of energy available to the cell directly correlates with the quality of the nutrient, denoted as $n_s$ [19,20,21,22]. Where the metabolic enzymes involved in fermentation and respiration processes generate different amounts of energy. For this reason, we consider that the nutrient quality is distributed differently for the fermentation and respiration metabolic pathways. Thus, the nutrient quality will be comprised of the nutrient quality used in the fermentative pathway expressed as $n_f$ and the nutrient quality used for the respiratory pathway expressed as $n_r$ [18, 19] (see Supplementary Information, Appendix A for more details). Additionally, Fig. 1B shows the dynamics of cellular energy, providing a clear perspective on its impact on cellular metabolism.

Fig. 1

General overview of the multiscale model that considers cellular mechanisms in the dynamics of a bioreactor. A The substrate S is introduced into the cell via transporter enzyme $e_t$, transforming it into internal substrate $S_i$. This internal substrate is subsequently converted into energy a through the action of fermentation metabolic enzyme $e_{m_f}$ and respiratory metabolic enzyme $e_{m_r}$. This energy is primarily used in transcription processes, generation mRNA molecules $m_x$, which then bind to free ribosomes r to form the ribosome–mRNA complex $c_x$. This complex participates in translation, achieving production of the six proteins that make up the model. The acetate Ac is secreted to the outside of the bioreactor. However, the acetate is assimilated by the cell and transformed in internal acetate $Ac_i$ when the levels of the main substrate S drop to near-zero for subsequent transformation into energy a. B Cell energy levels depend mainly on the quality of the supplied substrate $n_s$ and respiratory $e_{m_r}$ and fermentation $e_{m_f}$ metabolic enzymes and to a lesser extent, on acetate consumed by the cell and its quality $n_{s_{Ac}}$

The dynamics of the internal substrate are directly influenced by the substrate import rate to the bioreactor, $v_{imp}$, as well as by the catalytic rate of the fermentation metabolic enzyme, $v_{cat_f}$ and catalytic rate of the respiration metabolic enzyme, $v_{cat_r}$. These factors culminate in the derivation of the following mathematical expression:

$$\begin{aligned} \frac{dS_i}{dt}=v_{imp}(e_t,S)-v_{cat_{f}}(e_{m_f},S_i)-v_{cat_{r}}(e_{m_r},S_i)-\lambda S_i, \end{aligned}$$

(10)

where S represents the substrate concentration in the bioreactor and $S_i$ represents the internal substrate concentration. All kinetic rates adopt the Michaelis–Menten form. The parameters $v_{cat_f}$ and $v_{cat_r}$ are detailed in Fig. 1B, while the substrate import rate, $v_{imp}$, is defined as follows:

$$\begin{aligned} v_{imp}=\frac{e_t K_{cat} S}{K_m + S}, \end{aligned}$$

(11)

where $K_{cat}$ and $K_m$ correspond to the catalytic and saturation constants, respectively.

The formation of acetate is postulated as a physiological response to metabolic overflow [20]. In this context, fermentation metabolic enzymes directly participate in the synthesis and utilisation of acetate as a new carbon source once the main substrate has been consumed [17, 23, 24].

The consumption of acetate by the cell begins when the concentration of the main substrate drops below $2\%$ of its initial level [10]. Acetate is considered to energetically contribute to the cell due to its nutritional value, which is quantified as $n_{s_{AC}}$ (see Supplementary Information, Appendix B). The acetate intracellular level is described by the following mathematical expression:

$$\begin{aligned} \frac{dAc_i}{dt}=v_{con_{Ac}}(e_{m_f},Ac)-v_{cat_{Ac}}(e_{m_f},Ac_i)-\lambda Ac_i , \end{aligned}$$

(12)

where $Ac_i$ represents the concentration of acetate within the cell, Ac represents the concentration of acetate in the bioreactor, $v_{con_{Ac}}$ the acetate consumption rate and $v_{cat_{Ac}}$ the catalytic rate related to acetate. All kinetic rates involved in acetate consumption, acetate production and energy generation from acetate follow Michaelis–Menten kinetics, according to the following expressions:

$$\begin{aligned} v_{con_{Ac}}&= \frac{e_{m_f} K_{cat_{con_f}} Ac}{K_{m_{con_{f}}} + Ac}, \end{aligned}$$

(13)

$$\begin{aligned} v_{cat_{Ac}}&= \frac{e_{m_f} K_{cat_{Ac_i}} Ac_i}{ K_{m_{Ac_{i}}} + Ac_i}, \end{aligned}$$

(14)

where $K_{m_{con_f}}$ and $K_{m_{con_f}}$ correspond to the saturation constants for acetate consumption and internal acetate consumption, respectively. Meanwhile, $K_{cat_{con_f}}$ and $K_{cat_{Ac_i}}$ correspond to the catalytic constants for acetate consumption and the conversion of internal acetate into energy, respectively.

Finally, to introduce the production of a recombinant protein of interest as part of the mechanistic model, we consider a heterologous protein denoted as p which expresses this protein of interest, consuming cellular resources and directly affecting cell growth [11]. In this way, the equations related to the production of a recombinant protein are analogous to Eqs. 2, 3 and 4, where this time $y \in \{e_{m_f}, e_{m_r}, e_t, q, p \}$ and $x \in \{ r, e_{m_f}, e_{m_r}, e_t, q, p \}$.

Based on the above, it is established that the total amount of proteins in the cell is limited, with the proteomic fraction summing to a total of 1. This is defined as follows for a recombinant strain:

$$\begin{aligned} 1=\phi _R + \phi _p + \phi _{e_{m_f}} + \phi _{e_{m_r}} + \phi _q + \phi _t. \end{aligned}$$

(15)

This way, the reaction list used by the model, taking into account the new variables, is presented in the Supplementary Information, Appendix G, Figure G1.

Dynamic model for a batch fermenter

We developed a batch fermentation model based on a mathematical analysis of the mass balance in this type of bioreactor. We assume that the volume of the culture broth remains constant ($\frac{dV}{dt}=0$) and the specific growth rate determines the change in the number of cells over the course of the fermentation time. The dynamics of cell growth are described by the following expression:

$$\begin{aligned} \frac{dN}{dt}=\lambda N. \end{aligned}$$

(16)

The substrate within the batch bioreactor is completely mixed in the fermenter and undergoes a import rate by the cell denoted as $v_{imp}$, which directly depends on the abundance of transport enzymes. Therefore, the dynamics of the substrate concentration follow the following differential equation:

$$\begin{aligned} \frac{dS}{dt}=-v_{imp}(e_t,S) \frac{N}{Y_\frac{S}{N}}, \end{aligned}$$

(17)

where $Y_\frac{S}{N}$ corresponds to the efficiency of substrate consumption by the cell.

The acetate is secreted by the cell and its dynamics in the bioreactor is captured from the following mathematical expression:

$$\begin{aligned} \frac{dAc}{dt}=v_{prod_{Ac}(e_{m_f},S_i)} N-v_{con_{Ac}(e_{m_f},Ac)} N, \end{aligned}$$

(18)

where $v_{prod_{Ac}}$ represents the acetate production rate and is directly dependent on the internal substrate $S_i$, fermentation metabolic enzymes $e_{m_f}$, catalytic constant for acetate production $K_{cat_{prod_f}}$ and saturation constants for acetate production $K_{m_{prod_f}}$, following the Michaelis–Menten equation as described by the following expression:

$$\begin{aligned} v_{prod_{Ac}}=\frac{e_{m_f} K_{cat_{prod_f}} S_i}{ K_{m_{prod_f}} + S_i}. \end{aligned}$$

(19)

We define that intracellular enzymes are not secreted by the cell. Consequently, we characterise the production rate of a recombinant protein as the number of molecules of the heterologous protein p produced per cell during the generation of cell growth [10, 11]. Thus, the total quantity of heterologous protein molecules will increase in parallel with the growth of the cell population. Accordingly, the concentration of the protein in the bioreactor can be expressed by the following mathematical equation:

$$\begin{aligned} \frac{dP}{dt}=\lambda p N. \end{aligned}$$

(20)

Finally, the mechanistic nature of the model allows for the direct incorporation of tunable parameters commonly used in circuit design, such as gene induction strength or ribosomal binding site strength (see “Methods” section for further details). Research suggests that increases in gene induction rate and RBS strength lead to a reduction in cellular growth due to the limited availability of free ribosomes [14]. For this reason, through a series of simulations, we will establish how variations in these parameters affect the volumetric productivity of a batch bioreactor. Consequently, the volumetric productivity was calculated using the following expression:

$$\begin{aligned} P_{V_X}=\frac{X}{t_f}, \end{aligned}$$

(21)

where X represents the concentration in grams per liter of biomass, recombinant protein, and acetate. The final fermentation time, denoted as $t_f$, is reached when the biomass attains the stationary phase The biomass concentration is calculated based on the relationship between cell dry weight and specific growth rate, as discussed in the Supplementary Information, Appendix C and D, and detailed in Equation D5.

Fitting and validation of the multiscale model

We fitted the model to four distinct strains of Escherichia coli, namely W3110, VAL22, VAL23 and VAL24, as reported by Lara et al. [17]. The VAL22 strain is a unique mutant with an inactivated poxB (EC 1.2.5.1) gene encoding pyruvate oxidase (EC 1.2.3.3). VAL23 is a double mutant incapable of producing lactate and formate due to deletions in the ldhA (EC 1.1.1.27) and pflB (EC 2.3.1.54) genes, which encode lactate dehydrogenase (EC 1.1.1.27) and pyruvate-formate lyase (EC 2.3.1.54), respectively. VAL24 is a triple mutant with deleted ldhA (EC 1.1.1.27) and pflB (EC 2.3.1.54) genes, along with an inactivated poxB (EC 1.2.5.1) gene. W3110 corresponds to the parental strain. In all cases, glucose was considered as the substrate and a green fluorescent protein GFP was chosen as the desired product. Table 1 presents the operational conditions of the bioreactor used experimentally for each strain [17] and all the parameter values estimated and used by the model are found in the Table F3 and Table F4 respectively present in the Supplementary Information.

Table 1 Initial conditions of the batch bioreactor

The model was adjusted using biomass and substrate concentrations reported by Lara et al. [17] for each strain (see Fig. 2A). The parameter fitting was based a objective function F:

$$\begin{aligned} F=\sum ^{n}_{i=1} \left[\left(N_{i}-{\hat{N}}_{i}\right)^2 +\left(S_{i}-{\hat{S}}_{i}\right)^2\right], \end{aligned}$$

(22)

where n corresponds to the number of observations, $N_i$ and $S_i$ denote the observed variables for biomass and substrate, respectively and $\hat{N_i}$ and $\hat{S_i}$ represent the model predictions for biomass and substrate, respectively.

By adjusting the multiscale model to focus on biomass and substrate concentrations, we were able to predict the dynamics of GFP and acetate, as shown in Fig. 2B. This predictive capability was validated by comparing it with published experimental data, highlighting the model’s adaptability to various cultivation conditions. However, it is important to recognise that while the model captures general trends, discrepancies may exist in the dynamics of acetate and GFP, which can be attributed to metabolic complexities that are not explicitly represented in the model. Despite these limitations, the model allowed us to adjust parameters and develop optimal strategies to maximise the production of the desired metabolite.

Figure 2C presents a comparison between the multiscale model and an unstructured model that applies Monod kinetics [5] (see “Methods” section). We used the objective function F, specified in Eq. 22, to fit the biomass and substrate concentrations for the four Escherichia coli strains previously examined. A low value of F indicates a better model fit. In this way, it is evident that the unstructured model achieves greater accuracy in fitting for the four strains, showing more significant variability in substrate concentration compared to biomass, except in the VAL24 strain where the opposite occurs.

The comparison revealed that the multiscale model does not outperform the unstructured model in terms of fitting the objective function. However, it was observed that the W3110 strain exhibits very similar results in both models. Thus, due to the similarity in the fitting of the objective function, we decided to consider this strain as the selected one to be used in the multiscale model for the future simulation analyses that will be developed in this research.

On the other hand, unlike the unstructured model, the selection of the multiscale model offers clear advantages in terms of applicability, as it allows for the incorporation of a wide range of variables and parameters that describe cellular behaviour in detail. These characteristics make the multiscale model an essential tool for advancing future studies and developments in the field of biotechnology.

The fitting and subsequent validation of the model entail the estimation of a series of parameters that are part of the multiscale model. For this reason, we performed a sensitivity analysis for each estimated parameter [25] and represented it in a spider plot shown in Fig. 2C (details in “Methods” section). The results reveal parameters with low sensitivity, such as $k_{\gamma }$, which is consistent with the sensitivity analyses conducted by Weisse et al. [11]. In contrast, the transcription rate of recombinant proteins $w_p$ shows greater sensitivity due to its direct influence on the allocation of limited cellular resources. Specifically, a high value of $w_p$ leads to a greater allocation of limited resources dedicated to the expression of the heterologous protein, thereby reducing essential metabolic activities for the cell that affect cell growth [14].

Fig. 2

Model validation based on experimental data reported by Lara et al. [17]. A This figure shows the model’s ability to capture the growth dynamics of the strains by adjusting biomass and substrate concentrations. Both parameters were used to calculate the objective function F, which allowed for the model fit. Circles represent experimental data, while lines show the fitted results from the multiscale model. B Concentrations of GFP and acetate obtained from the model’s predictive capability for each analysed strain. C Comparison of the objective function for the four strains, considering an unstructured model that follows Monod kinetics versus the multiscale model. In both cases, the objective function only considered biomass and substrate concentrations. A small F value indicates a close fit between the model and the reported experimental data. The spider plot represents the sensitivity analysis conducted using the extended Fourier amplitude sensitivity test (eFAST) [25], for the six estimated parameters: the translation rate of the recombinant protein $w_p$, transcription rate of the enzymes responsible for fermentation $w_{e_{mf}}$, transcription rate of the enzymes responsible for respiration $w_{e_{mr}}$, transcription rate of the substrate transporter enzymes $w_t$, quality of the substrate used in the fermentation pathway $n_f$ and translation elongation threshold $K_\gamma$

Finally, based on the obtained results, we propose the following steps for other proteins and expression hosts for which temporal data is available: (1) identify the strain, culture type, conditions, and recombinant protein of interest, (2) collect experimental data on the dynamics of biomass and substrate, with a particular focus on substrate selection, as the kinetic parameters employed in Eq. 11 vary depending on the type of substrate used, (3) adjust the model to the experimental data of biomass and substrate using Eq. 22, which will allow for the estimation of the six parameters presented in the spider plot in Fig. 2C (if the fit is inadequate, additional parameters for estimation may need to be incorporated, as detailed in Supplementary Information Appendix E and Figure G3), and (4) simulate the multiscale model using the new parameter estimates to predict recombinant protein and acetate production.

Analysing productivity in response to heterologous gene induction and ribosomal binding strength

The mechanistic nature of the model allows for the direct incorporation of adjustable parameters, such as the induction strength of a heterologous gene $w_p$ or the ribosomal binding site strength, which are widely used in the design of synthetic circuits [10, 14, 15]. This integration enables a direct linkage of cellular behaviour to productivity analysis within a bioreactor.

Figure 3A presents a comparison between two strains: a wild-type strain and a recombinant strain expressing proteins of interest. The wild-type strain, in the absence of a heterologous gene, allocates all its limited resources to cellular growth, thereby increasing the proteomic fraction of ribosomes $\phi _R$, which is directly linked to an increase in the specific growth rate [13] (details in “Methods” section). On the other hand, the expression of a protein of interest by a recombinant strain does not provide any benefit to the cell as the expression of this protein merely consumes cellular resources, occupying a significant portion of the cellular proteome and reducing the availability of other proteins that are essential for cellular growth and maintenance [14].

Fig. 3

Productivity as a function of recombinant protein induction and changes in ribosomal binding strength. A Proteomic fraction assignment is performed for a wild-type strain and a recombinant strain. The expression of a heterologous gene does not provide any contribution to the cell and only consumes limited cellular resources, occupying part of the proteome. The proteomic fraction of the rest of the proteins undergoes a significant decrease, directly affecting cellular growth. Here, $\phi _R$, $\phi _q$, $\phi _{et}$, $\phi _{e_{m_f}}$, $\phi _{e_{m_r}}$ and $\phi _p$ correspond to the proteomic fraction of ribosomes, maintenance enzymes, transport enzymes, fermentation enzymes, respiration enzymes and recombinant protein, respectively. B An increase in the expression of a heterologous gene leads to a decrease in cell growth, biomass and acetate productivity. However, at very high induction levels, the cell cannot sustain the production of recombinant proteins due to the lack of free ribosomes, depicted in a pie chart for the entire analysed induction range. C There is a delicate balance between strength RBS and the induction rate of the heterologous gene, ensuring elevated levels of productivity and shorter cultivation times. D Productivity strategies were considered by optimising biomass and recombinant protein production, with both variables considered as part of the objective function J. A weighting variable $\alpha$ was introduced with a range between 0 and 1 to determine maximum biomass productivity ($\alpha$=1) and maximum GFP productivity ($\alpha$=0)

Therefore, it is crucial to understand the cellular response to variations in the induction of a heterologous gene in order to extrapolate the results to a bioreactor through volumetric productivity analysis [10]. Figure 3B illustrates the trends in specific growth rate, biomass productivity, GFP production and acetate production at different levels of induction of the heterologous gene.

An increase in the induction rate $w_p$ promotes a more pronounced distribution of limited cellular resources towards GFP production, reducing growth rate and biomass and acetate productivity. Figure 3B shows how biomass yield is affected by an increase in GFP productivity, a consequence of an elevation in the induction rate, reaching its peak in a proximity range of induction around 1500 [molecules/min/cell]. However, surpassing this value, GFP productivity begins to decline due to the scarcity of available ribosomes, which play a crucial role in translational activity. Initially, free ribosomes are primarily engaged in metabolic activities compromising cellular growth [13]. Nevertheless, as we increase the induction rate, the count of free ribosomes decreases sharply, reaching levels close to zero, halting cellular growth and essential metabolic functions.

The ribosomal binding site strength plays an essential role in the design of synthetic circuits, as it regulates the allocation of ribosomes to mRNA and thereby facilitates translation [26,27,28]. Due to its relevance, we have investigated its impact on the volumetric productivity of a batch bioreactor, examining the results under various induction rates of a recombinant protein. The data presented in Fig. 3C show the observed trends in specific cell growth rate, biomass and GFP productivity, as well as cultivation times. These findings reveal a trade-off between these parameters, which is crucial for optimising process productivity.

A high RBS strength increases competition for free ribosomes in the cell, capturing a greater proportion of these ribosomes and thereby reducing the availability of free ribosomes for other cellular activities [14, 15, 29, 30]. This can adversely affect general cellular functions and metabolic balance, which is manifested by a significant decrease in cell growth, volumetric productivity and high cell culture times.

Figure 3C illustrates the balance between ribosomal binding site (RBS) strength and induction $w_p$. It reveals that high levels of induction of a heterologous gene require a reduction in RBS strength, which in turn increases the availability of free ribosomes for metabolic activities essential for cellular growth, thereby enhancing biomass productivity. In contrast, the productivity of a heterologous protein is highly sensitive to variations in the induction rate. Low induction levels result in significantly reduced productivity due to insufficient activation, necessitating an increase in RBS strength to improve the affinity between the ribosome and the mRNA of the heterologous protein, thereby boosting its production. On the other hand, excessively high induction rates also lead to decreased productivity due to a shortage of free ribosomes, as detailed in Fig. 3B. Therefore, under high induction conditions, it is necessary to reduce RBS strength to allow for rapid dissociation between the ribosome and the mRNA of the heterologous protein. This facilitates the translation of other mRNA molecules, increasing the cell’s capacity to synthesise other essential proteins for its functioning.

Additionally, Fig. 3C shows that the highest values of volumetric productivity for recombinant protein, biomass, and cell growth are achieved when RBS strengths below 1 are considered. This highlights the importance of an optimal balance between heterologous gene induction and RBS strength to maximise productivity. For this reason, we have implemented an optimisation procedure to determine the optimal RBS strength and induction $w_p$ to maximise the productivity of recombinant protein and biomass. The optimisation involved creating a quadratic error objective function J, where biomass productivity and GFP production served as optimisation targets. This approach facilitates the development of a production strategy that simultaneously emphasises both biomass generation and GFP production, addressing various objectives within the bioreactor process.

$$\begin{aligned} \begin{array}{cc} \underset{k_b,w_p}{\text{minimise}} &\quad J(x,k_b,w_p), \\ \text {subject to:} &\quad dx/dt = h(x), \\ & \quad k_b,w_p \ge 0. \end{array} \end{aligned}$$

(23)

The objective function is of the form:

$$\begin{aligned} J=\alpha (N_{obj}-N_{max})^2+ (1-\alpha )(P_{obj}-P_{max})^2, \end{aligned}$$

(24)

where $N_{obj}$ corresponds to the maximum desired biomass productivity, $N_{max}$ corresponds to the maximum biomass productivity achieved by the model, $P_{obj}$ corresponds to the maximum desired recombinant protein productivity, $P_{max}$ corresponds to the maximum protein productivity achieved by the model and $\alpha$ corresponds to the weight of the objective function. The weight $\alpha$ can be used to control maximisation in recombinant protein or biomass production, the weight $\alpha$ plays a crucial role in our optimisation process, where a value of $\alpha =1$ indicates an exclusive focus on maximising biomass productivity, while a value of $\alpha =0$ signals an absolute prioritisation of GFP productivity.

Figure 3D illustrates a trade-off between RBS strength for values close to 1 and $w_p$, highlighting the direct relationship between the induction of the heterologous gene and the productivity of the recombinant protein. As induction values increase, a higher RBS strength is required to ensure optimal productivity, as this enhances the efficiency of the translation processes involved in the expression of the heterologous gene.

In this context, competition for limited resources and low induction generate a linear relationship within the induction range of 600 to 800 [molec/cell/min]. This behaviour suggests the need to increase the affinity between ribosomes and the mRNA of the heterologous gene to favour the production of heterologous proteins. In contrast, when the goal is to maximise biomass productivity in the absence of heterologous gene induction, cellular resources are directed toward growth and maintenance, eliminating the need for high RBS strength values. Since the translation of the mRNA from the heterologous gene is not active in this scenario, a lower RBS strength reduces the metabolic burden, allowing the cell to optimise resource usage for cellular growth.

The results suggest that prioritizing biomass productivity over recombinant protein productivity yields substantial outcomes. We observed an $86\%$ increase in the proteomic fraction associated with ribosomes and a potential $60\%$ improvement in their productivity. As explained in the previous paragraph, considering low induction rates allows for working with lower RBS strength values to avoid metabolic burden and promote resource allocation for cellular growth. Conversely, a preference for recombinant protein production intensifies cellular stress due to higher induction rates, requiring higher RBS strength. Notably, a strategy favouring recombinant protein production over biomass results in an extraordinary $93\%$ increase in GFP productivity, constituting $86\%$ of the proteomic fraction.

These insights underscore the inherent strategic flexibility in balancing biomass and recombinant protein productivity based on diverse objectives and operational considerations within the bioreactor process.

Analysis of acetate productivity against changes in the induction rate of the heterologous gene and the transporter enzyme

Previous studies have established a direct relationship between acetate production and substrate import rates, with an observed increase in acetate production as substrate consumption rates rise [23, 31]. Additionally, research has shown that it is feasible to regulate acetate production through the control of the induction of transport enzymes [32,33,34]. Given that our model can predict both acetate and recombinant protein production, we decided to directly explore the effects of varying the induction rate of a heterologous gene, $w_p$, as well as the induction of the transport enzyme, $w_t$, on productivity and specific growth rate.

In the first phase of the simulation, we analysed the behaviour of maximum acetate productivity in response to changes in the induction rate of the recombinant protein, considering different specific cell growth rates under low, medium, and high induction conditions. Figure 4A shows a direct and proportional relationship between protein induction and acetate productivity. The maximum levels of acetate productivity are reached when the induction levels of the heterologous gene are low or absent (wild-type strain). Under these conditions, the cell allocates its resources to growth and metabolism, naturally producing acetate as a byproduct.

We also found that the threshold growth rate required to initiate acetate formation is lower in cells induced to produce a recombinant protein, as demonstrated in previous studies [19, 35]. Figure 4A further compares the maximum acetate productivity with that of the recombinant protein, showing that increases in induction negatively impact acetate production, primarily due to the limitation in the availability of free ribosomes.

Since acetate production is correlated with specific growth rate and substrate consumption rate [20, 36], we conducted a second simulation to explore how specific growth rate and acetate productivity are affected by changes in the induction of the recombinant protein and the transport enzyme. In this context, Fig. 4B illustrates the increase in substrate import rate as an initial condition in response to an increase in $w_t$. This increase in the parameter elevates the concentration of transport enzymes within the cell, enhancing the consumption of the main substrate.

Fig. 4

Productivity as a function of recombinant protein induction and transporter enzyme induction. A Acetate productivity as a function of specific growth rate against variations in the induction rate of the recombinant protein. The blue line indicates the general trend, showing an increase in acetate productivity with increasing rates of recombinant protein induction. The colored circles represent different levels of induction: green for null induction corresponding to the wild-type strain, orange for medium induction, and yellow for high induction. The inset highlights the correlation between maximum GFP and acetate productivity as a function of recombinant protein induction, establishing the limitations generated by limited resources. B The effect of transporter enzyme induction on maximum acetate productivity and maximum specific growth rate. The lines represent different rates of recombinant protein induction, where the green line corresponds to the wild-type strain with null induction, the orange line represents medium induction, and the yellow line indicates high induction. The red points indicate acetate productivity at different levels of transporter enzyme induction, with the size of the circle proportional to the magnitude of induction. The inset shows the substrate import rate against variations in induction levels, highlighting an exponential increase in substrate consumption at higher levels of transporter enzyme induction. C A trade-off between recombinant protein induction and transporter enzyme induction is necessary to ensure maximum productivity and cellular growth. This balance ensures that both factors are maximised under the limitations of available cellular resources, optimising both biomass yield and recombinant protein production

Additionally, Fig. 4B examines the variation in acetate productivity as a function of changes in specific growth rate, limited by the induction of transport enzymes ($w_t$). A direct relationship was observed between these two parameters, where maximum acetate productivity is achieved at higher specific growth rates due to increased substrate consumption by the cell, allowing more carbon to enter the cell and be converted into energy for cellular activities. However, the overexpression of the recombinant protein through the induction of $w_p$ reduces the maximum acetate productivity due to the limitations imposed by the availability of free ribosomes, which are linked to the high production of GFP.

Finally, we simulated the behaviour of maximum biomass productivity, GFP, acetate, and specific growth rate in response to variations in the induction rate of the recombinant protein and the transport enzyme. Figure 4C shows a balance between both parameters that ensures higher values of productivity and cellular growth. High induction rates of the recombinant protein are counterbalanced by high induction rates of the transport enzyme, thereby enhancing substrate import and consequently increasing productivity and cellular growth.

Discussion

This study presents a mechanistic model capable of predicting the dynamics of recombinant protein and acetate evolution during a bioreactor culture process, integrating intracellular complexity and the universal cellular trade-offs related to proteome, energy, and ribosome limitations with the operational variables at bioreactor level. Unlike classical phenomenological models (e.g. Monod model), our approach enables the regulation of key synthetic biology parameters, such as heterologous gene induction and ribosome binding site strength, allowing for a more accurate representation of the metabolic processes involved in recombinant protein production. The multiscale model is grounded in a mechanistic representation of universal cellular trade-offs, including limited resources of energy, ribosomes, and proteins, which govern critical processes such as protein synthesis and gene expression. By reducing the demand for any of these resources, capacity is freed for the development of other intracellular processes.

Since our model integrates the internal behaviours of the cell with the operational parameters of a bioreactor, we generated various simulation scenarios. Initially, we compared the results for a wild-type strain and a recombinant strain. The comparison reveals that the expression of heterologous proteins can divert critical resources toward their production, negatively impacting cell growth by reducing the availability of ribosomes for essential metabolic functions. This phenomenon, known as “metabolic burden” and documented in the literature [13, 14, 19], suggests a necessary trade-off between protein production and cell growth, as evidenced by the decrease in biomass and acetate productivity as heterologous gene induction increases.

In this way, our results highlight the importance of finding an optimal balance between gene induction and RBS strength to maximise both biomass and recombinant protein production. An excessively strong RBS can exacerbate competition for ribosomes, thereby limiting the overall efficiency of the system. In this regard, optimising these parameters leads to an optimal productivity point, where ribosome availability is sufficient to support both cell growth and protein synthesis. This approach not only enables the prediction of productivity based on adjustable variables but also provides a valuable framework for the rational design of bioprocesses, contributing to improved efficiency in industrial protein production. The strategic flexibility in managing resources and operational parameters underscores the potential of this model to guide the optimisation of biotechnological processes.

Additionally, the model is capable of predicting acetate concentration in the bioreactor. For this reason, we studied the relationship between its production and the overexpression of heterologous proteins. The results confirm that acetate production is strongly influenced by heterologous gene induction and the availability of free ribosomes. Cells that do not express a recombinant protein tend to allocate more resources toward cell growth and acetate production as a byproduct. However, it is possible to observe that cells can reach high acetate productivity levels even at low cell growth rates when there is overexpression of a recombinant protein, as previously reported in the literature [19, 37]. In this context, the overexpression of heterologous proteins causes saturation of the cellular machinery, reducing the availability of proteins essential for central metabolism. As a result, the carbon flux is redirected towards fermentative pathways instead of respiratory ones, promoting metabolic overflow and increasing acetate production [35]. Essentially, acetate production is the metabolic response to the stress caused by heterologous protein overexpression, which explains the high productivity levels at low growth rates. This highlights the complex interaction between cellular processes and their impact on metabolite production, emphasizing the competition for limited cellular resources.

A key finding is the model’s ability to predict that acetate productivity can be mitigated or controlled through the regulation of transporter enzymes. An increase in the induction of these enzymes allows for greater substrate import, which, combined with moderate levels of recombinant protein induction, facilitates a better balance between acetate production and cell growth. This approach suggests a potential strategy for optimising biotechnological processes by preventing the accumulation of undesirable acetate while maximising the production of recombinant proteins and biomass.

Despite the positive results obtained with the multiscale model, some limitations arise, primarily related to the simplification in the representation of complex metabolic networks. Although the model provides accurate predictions regarding biomass, recombinant protein, and acetate production, it does not fully capture the network-level metabolic interactions that could influence the distribution of carbon and energy flows. However, we believe that the implementation of the model is highly useful for understanding how the overexpression of heterologous proteins leads to increased acetate production as a consequence of metabolic overflow. This limitation could be addressed in future studies by integrating metabolic flux analysis approaches, allowing for a more detailed interpretation and optimisation of the use of limited cellular resources.

Regardless of these limitations, we believe that the implementation of this type of model in bioprocesses represents a crucial first step toward designing more efficient strains. These models not only enable increased productivity but also provide a deeper understanding of cellular physiology in response to metabolic overloads. For example, Atkinson et al. [16] demonstrated, through a mechanistic model based on the work of Weisse et al. [11], that task allocation in microbial consortia for the degradation of complex substrates can significantly reduce the genetic burden on each cell, thereby improving both growth and overall system efficiency. Likewise, Faizi et al. [38] developed a model for allocating limited resources, identifying optimal strategies for proteome distribution, thus achieving greater productivity in phototrophic cultures under limiting light conditions in a chemostat. These examples illustrate how proteomic allocation models have already been successfully implemented, and therefore, we aspire for our research to serve as a foundation for developing new design and operational strategies. By integrating strain design and bioreactor operation, it will be possible to maximise both productivity and efficiency under complex industrial conditions.

Conclusions

This study introduces a refined multiscale model that builds upon the framework provided by Weisse et al. [11], incorporating significant new variables such as acetate and a the production of a recombinant protein. This advanced approach surpasses classical models, such as the Monod model, by integrating a broader spectrum of internal cell parameters that facilitate a detailed description and a profound understanding of cellular growth and production mechanisms.

The application of our model has demonstrated its ability to accurately predict how adjustments in the induction rate of proteins and the strength of the ribosome binding site affect the production rates of biomass, recombinant proteins and acetate within a bioreactor. This knowledge enables the optimisation of biotechnological processes, which is essential for the industrial and therapeutic production of recombinant proteins.

The results also show that careful regulation of heterologous gene induction and optimisation of RBS strength are crucial for maximising productivity, maintaining an optimal balance between biomass growth and protein expression. This underscores the value of our model as an essential tool for bioprocess engineering, proposing methods to enhance both the efficiency and sustainability of these processes.

In conclusion, the proposed multiscale model not only provides a comprehensive understanding of molecular interactions and bioreactor productivity but also establishes a solid framework for future improvements in biological engineering strategies. The advancements achieved represent a significant milestone toward the development of more effective and sustainable biotechnological processes, highlighting the importance of this study in advancing the field of biotechnology.

Methods

Implementation and configuration of simulation in the multiscale model All the differential equations in this research, corresponding to the multiscale model, were developed using the ODE15s differential equation solver available in the Matlab software. For the simulation, the initial conditions of the cellular variables were initially calculated, with the substrate held constant until reaching a steady phase. Once this state was achieved, the resulting values of these cellular variables were used as initial conditions for the multiscale model.

Parameter estimation The model was adjusted using experimental data proposed by Lara et al. [17]. In this adjustment, we worked with a total of four different strains, corresponding to a W3110 parental strain and three recombinant strains called VAL22, VAL23 and VAL24.

To perform the fitting of these four strains, we employed the “genetic algorithm” function and the “ODE15s” solver for stiff differential equations, available in the Matlab R2021b software.

For model calibration, six different parameters were estimated: the translation threshold $K_{\gamma }$, the nutrient quality used in the fermentation pathway $n_f$ and the transcription speeds of the transporter enzymes $w_t$, metabolic fermentation enzymes $w_{e_{mf}}$, metabolic respiration enzymes $w_{e_{mr}}$ and the recombinant protein $w_p$.

Only the concentration of biomass and substrate was adjusted, as shown in Fig. 2A, using the function F, defined in Eq. (22), through the genetic algorithm implemented as a function in Matlab. Figure 2B illustrates the predictive capability of the model, where, based on the performed adjustment, the dynamics in the concentration of GFP and acetate are effectively represented.

The estimation results for each strain are presented in Table F3, of the Supplementary Information.

Multiscale and Monod model comparison

Figure 2C presents a comparison of the biomass and recombinant protein adjustments using the objective function F, as defined in Eq. (22). The equations for the Monod model used for the fitting are as follows:

$$\begin{aligned} \frac{dN}{dt}&= \lambda N, \end{aligned}$$

(25)

$$\begin{aligned} \frac{dS}{dt}&= -\lambda \frac{N}{Y_{\frac{S}{N}}}, \end{aligned}$$

(26)

where the specific growth rate is given by $\lambda =\frac{\lambda _{\text {max}}S}{K_s + S}$, with $\lambda _{\text {max}}$ as the maximum specific growth rate and $K_s$ as the half-velocity constant [5].

Sensitivity analysis Figure 2C was constructed through a sensitivity analysis using the Fourier amplitude sensitivity test (eFAST) method, as developed by Marino et al. [25]. For each estimated parameter, the search curve was repeated four times, totaling 65 evaluations. Values closer to one in the results indicate higher sensitivity in the model.

Analysis and optimisation of bioreactor productivity by varying the induction rate $w_p$ and the RBS strength

To assess the variation in the productivity of a bioreactor in relation to the induction rate $w_p$, we explored a wide range of this parameter [14], ranging from $[10^0, \ 10^5]$. This range allowed us to analyse in greater detail the effects of induction of the heterologous gene and identify inflection points caused by the scarcity of free ribosomes, as evidenced in Fig. 3B.

Regarding the pie charts representing the availability of free ribosomes, we employed the expression derived from the work of Weisse et al. [11]:

$$\begin{aligned} R = r + R_t, \end{aligned}$$

(27)

where R represents the total number of ribosomes, r corresponds to the free ribosomes and $R_t$ indicates the ribosomes that are bound to the messenger RNA.

In Fig. 3C, a heat map is presented where cell growth, productivity and cultivation time vary based on $w_p$ and RBS strength. To conduct this simulation, we considered a range of values for $w_p$ identical to that analysed in Fig. 3B. However, concerning RBS strength, defined as the ratio between the association constant $k_b$ and dissociation $k_u$ ($\text {RBS}=\frac{k_b}{k_u}$), we varied only the association constant $k_u$ in the interval $[10^{-4} , \ 1]$, keeping $k_b$ constant at a value of 1 [14].

Regarding the optimisation presented in Fig. 3D, we considered the range of RBS and $w_p$ used in previous analyses. In the development of Eqs. 23 and 24, we applied a genetic algorithm implemented as a function in Matlab and simulated various scenarios given by the variation of a weighting parameter $\alpha$ with a range of values $[0 , \ 1]$.

The mass fraction of each protein in the productivity strategy of Fig. 3D was calculated following the expressions:

$$\begin{aligned} \phi _y&= \frac{n_y y}{M}, \end{aligned}$$

(28)

$$\begin{aligned} \phi _R&= \frac{n_r R}{M}, \end{aligned}$$

(29)

where $y \in \{ {e_{m_f}, e_{m_r}, e_{t}, p } \}$ and R corresponds to the total ribosomes in the cell, $n_y$ and $n_r$ denote the lengths of non-ribosomal proteins and ribosomes, respectively. The sum of all proteomic fractions must be equal to 1.

$$\begin{aligned} \phi _q= 1-\phi _R - \sum \phi _y. \end{aligned}$$

(30)

The values for cell growth and productivity correspond to the maximum output provided by the simulation for each scenario.

Productivity as a function of the induction of the transporter enzyme and the recombinant protein

In Fig. 4A, we allowed the specific growth rate to be a function of induction $w_p$ in the range of $[0 , \ 10^3]$, considering the maximum acetate productivity for each simulation.

On the other hand, Fig. 4B depicts acetate productivity and specific growth rate in relation to the induction of the transporter enzyme. Initially, we illustrate how the import rate of the main substrate increases with variations in $w_t$. We considered a range of $w_t$ from $[1 , \ 10]$, using a value of 1 to avoid model deceleration due to low consumption rates and adopting a maximum value of 10, as no significant changes were observed beyond this point. Additionally, to study the effect of overexpression of recombinant proteins on acetate productivity, we considered a smaller induction range of $w_p$ between $[0 , \ 2000]$ in order to avoid the cellular overload caused by the expression of the heterologous protein.

Finally, Fig. 4C presents a heatmap for cell growth and biomass, GFP and acetate productivity. For its development, we used an induction rate ranging from $[1 , \ 10^5]$ and a $w_t$ range from $[1 , \ 10]$.