Optimal Resource Allocation in a Bacterial Growth Model Under Cold Stress and Temperature

1  Introduction

Microorganisms are constantly exposed to fluctuations in their environment, particularly in nutrient availability and temperature. To survive and thrive under such conditions, they have evolved sophisticated strategies that enable dynamic physiological adaptation [1]. This adaptation is primarily achieved through large-scale reorganization of gene expression, including regulation of genes responsible for nutrient uptake, energy production, protein synthesis, and stress responses [2,3]. Central to this process is the efficient allocation of cellular resources, particularly protein synthesis capacity, across competing functional demands [4–6]. Among environmental variables, temperature is especially influential, directly affecting reaction kinetics, enzyme activity, and the stability of cellular structures [7]. Under cold stress, biochemical reactions slow markedly, and energy availability becomes limited, intensifying the challenge of maintaining balanced resource allocation [1].

Exposure to low temperatures triggers extensive physiological reorganization, including the redistribution of limited resources such as ribosomes, metabolic enzymes, and energy toward essential survival functions. In these conditions, messenger RNA (mRNA) molecules tend to form stable secondary structures, such as hairpins, that impede ribosome binding and slow translation, thereby reducing protein synthesis efficiency [1]. To overcome these obstacles, Escherichia coli induces cold-adaptive proteins (CSPs), particularly CspA, which is rapidly upregulated under low temperatures and can constitute up to 10% of total cellular protein, reaching intracellular concentrations around 100 μM [8–10]. Acting as RNA chaperones, CSPs destabilize inhibitory secondary structures in nucleic acids, facilitating efficient transcription and translation during cold adaptation and thereby restoring cellular homeostasis under low-temperature conditions.

CspA functions primarily as an RNA chaperone, binding transiently to single-stranded regions of mRNA with low sequence specificity and moderate affinity [8,11]. It binds cooperatively to single-stranded RNA stretches longer than 74 nucleotides, with a minimal concentration of approximately 2.7×10−5 M required for cooperative binding, well below the levels observed in cells under cold conditions (1×10−4 M) [8]. Functionally, CspA destabilizes inhibitory secondary structures of RNA, preventing the formation of stable hairpins that impede translation and making mRNA more accessible to ribosomes and ribonucleases. Multiple CspA molecules can bind concurrently to a single mRNA transcript, remodeling stabilized secondary structures at low temperatures and facilitating efficient translation initiation. By improving ribosome accessibility and translational efficiency, CspA plays a central role in maintaining protein synthesis capacity when low temperatures slow cellular processes [12].

Experimental studies indicate that low temperatures reduce translation efficiency due to stabilization of mRNA secondary structures, limiting protein synthesis and growth. To cope with this challenge, E. coli induces cold-shock protein CspA, which alleviates inhibitory RNA structures [9,12]. However, the synthesis of CspA consumes significant cellular resources, creating a trade-off. Excessive CspA production competes with ribosome biogenesis and other essential processes, potentially limiting overall growth. This suggests the existence of an optimal CspA level that balances the benefits of enhanced translation against the costs of resource allocation, providing a natural framework for applying optimal control theory to quantify growth strategies under cold-shock conditions. Optimal control provides a theoretical benchmark for quantifying such growth strategies under cold-shock conditions [13,14], and although exact optimal strategies are unlikely to be implemented physiologically, they represent a gold standard that can be used to evaluate how closely real bacterial regulatory mechanisms approximate optimal resource allocation.

Microbial growth and resource allocation have been extensively investigated from both experimental and theoretical perspectives [7,13,14]. Early systems-biology studies established empirical growth laws linking ribosome content, proteome partitioning, and cellular growth rate, revealing fundamental trade-offs in bacterial resource economics [15–17]. Building on these biological insights, several theoretical frameworks have sought to formalize microbial growth optimization using mathematical and control-theoretic approaches. For instance, Giuliodori et al. [12] formulated microbial growth in fluctuating environments as an optimal control problem, while Yegorov et al. [13,14] extended this work to include protein degradation, recycling, and temperature-dependent kinetics. More recently, Mairet et al. examined proteome allocation under temperature stress, showing that extreme conditions redirect resources to chaperone-mediated repair processes [18]. Despite these advances, quantitative models that specifically capture low-temperature adaptation and its associated trade-offs in Escherichia coli remain limited. Addressing this gap, the present study develops a dynamical model to investigate how bacteria allocate limited translational resources to balance growth and stress adaptation under cold-shock conditions.

Building on these principles, we formulate a dynamical model of E. coli under low temperature using ordinary differential equations (ODEs) to describe the interplay between ribosomes, metabolism, and CspA. The RNA-chaperone activity of CspA is represented by a modulation factor applied to the effective protein-synthesis rate, reflecting its role in alleviating inhibitory mRNA secondary structures. To focus on the primary trade-off between translation efficiency and ribosomal availability, we assume that other regulatory factors and stress responses not directly related to ribosome or CspA allocation are negligible. For the purpose of this model, we assume that CspA-mediated translation is the main mechanism driving adaptation under low temperatures. Under these assumptions, ribosome and CspA synthesis are dynamically regulated to maximize cumulative growth, defined as the total increase in cellular components (ribosomes, CspA, and metabolic intermediates), over a defined time horizon. Resource allocation is formulated as a dynamic optimal control problem in which the objective is to maximize cumulative cellular growth over a defined time horizon. In this framework, the control variables represent the fractional distribution of translational resources between ribosome synthesis and CspA production. Applying Pontryagin’s Maximum Principle yields the necessary conditions for optimal resource allocation and enables the identification of candidate strategies that balance growth with stress adaptation. This formulation supports the analysis of both open-loop and feedback control policies, offering insights into how temporal dynamics and environmental conditions shape optimal bacterial growth trajectories under cold-shock stress.

Overall, this study introduces a quantitative framework that integrates molecular cold-shock mechanisms with optimal resource allocation principles. By explicitly linking RNA-level adaptation to proteome allocation and growth optimization, the model offers new insight into how E. coli balances translational efficiency and biomass production at suboptimal temperatures. The results generate testable predictions for cold-shock adaptation and provide a foundation for future experimental and theoretical investigations into microbial growth strategies under low-temperature stress. The paper is organized as follows. Section 2 outlines the mathematical framework governing the dynamical model. Section 3 analyzes temperature-dependent steady states and stability. Section 4 introduces the optimal control problem and characterizes optimal allocation strategies. Section 5 presents the key conclusions and highlights potential directions for future research.

2  Formulation

Consider a population of self-replicating prokaryotic cells, such as Escherichia coli, in a continuous stirred-tank reactor (CSTR) of constant volume. Let P, M, R, and C [g] denote the total masses of precursor metabolites (amino acids), metabolic machinery (enzymes involved in the uptake and conversion of nutrients into precursors), gene expression machinery (polymerases and ribosomes), and the cold-shock protein CspA, respectively. The metabolic machinery converts external substrates into precursors, while the gene-expression machinery and CspA utilize these precursors to synthesize the macromolecules required for cellular growth and maintenance. A schematic diagram (Fig. 1) illustrates the mass fluxes and catalytic interactions within the system. Based on the mass balances of these components, the system dynamics are described by a controlled set of ordinary differential equations (ODEs) that capture the time evolution of precursors, metabolic enzymes, ribosomes, and CspA (Eq. (1)). Unlike classical microbial growth models in which translational efficiency is prescribed as a fixed or externally scaled parameter, the present formulation treats translation as an endogenous, state-dependent process modulated by the intracellular concentration of the cold-shock protein CspA.

dPdt=VM−f(C)VR,dMdt=(1−α)(1−β)f(C)VR−γMM(t),dRdt=α(1−β)f(C)VR−γRR(t),dCdt=βf(C)VR−γCC(t),(1)

Figure 1: Flow diagram represents the flow or distribution of nutrients between different parts of the cellular machinery.

with following initial conditions,

P(0)=P0≥0, M(0)=M0≥0,R(0)=R0≥0, C(0)=C0≥0,(10a)

where f(C) denotes the CspA-dependent modulation of translational efficiency, VM and VR are the metabolic and translational fluxes, respectively, t∈[0,T] (hours) is the time variable, and T is the time horizon. A detailed description of each parameter and function involved in the proposed model (Eq. (1)) can be found as follows:

•   VM=VM(t) [g⋅h−1] is the metabolic flux, i.e., the rate at which precursor metabolites P are synthesized by the metabolic machinery M from external nutrients.

•   VR=VR(t) [g⋅h−1], modulated by f(C), denotes the effective translational flux, i.e., the rate at which ribosomes R consume precursor metabolites to synthesize cellular proteins. In contrast to classical growth models where translational efficiency is prescribed as a fixed or externally scaled parameter, VR here is an endogenous, state-dependent quantity regulated by the intracellular concentration of the cold-shock protein CspA. The modulation function f(C) accounts for the RNA-chaperone activity of CspA under cold-shock conditions, whereby destabilization of inhibitory mRNA secondary structures enhances ribosome accessibility and protein synthesis efficiency.

•   γM>0 [h−1], γR>0 [h−1] and γC>0 [h−1] represent the degradation rates of the metabolic machinery M, gene expression machinery R and cold shock protein C, respectively.

•   β=β(t) and α=α(t) are dimensionless resource allocation functions (control strategies), such that at any time t∈[0,T], β(t) represents the fraction of precursor mass allocated to CspA, while the remaining fraction 1−β(t) is divided between ribosomes and metabolic machinery according to α(t) and 1−α(t), respectively.

•   α∈[0,1] and β∈[0,0.2] denote the admissible range of resource allocation fractions.

•   Cells in the population are assumed to share a constant cytoplasmic density, with δ>0 [L⋅g−1] representing its inverse.

•   The structural volume of the cell population is defined as

V(t)=β⋅(M(t)+R(t)+C(t))[L](2)

representing the volume occupied by the macromolecules of the metabolic and gene expression machinery and cspA, excluding monomeric precursors.

•   The concentrations of the cellular components are defined as

p(t)=P(t)V(t)[g.L−1],r(t)=R(t)V(t)[g.L−1],c(t)=C(t)V(t)[g.L−1],m(t)=M(t)V(t)=1δ−r(t)−c(t)[g.L−1]vR(t)=VR(t)V(t)[g.L−1.h−1],vM(t)=VM(t)V(t)[g.L−1.h−1].

Using the relation in Eq. (2), we can rewrite it as:

m(t)=1δ(1−r(t)−c(t)).(3)

•   The growth rate of the self-replicating system is

μ(t)=1V(t)dV(t)dt=δ[vRf(c)+(γM−γR)r(t)+(γM−γC)c(t)−γMδ].

Now, after the modification as mentioned above, the proposed model (Eq. (1)) can be written as follows:

dpdt=vM−f(c)vR−pδ[vRf(c)+(γM−γR)r(t)+(γM−γC)c(t)−γMδ],drdt=α(1−β)f(c)vR−γRr−rδ[vRf(c)+(γM−γR)r(t)+(γM−γC)c(t)−γMδ],dcdt=βf(c)vR−γCc−cδ[vRf(c)+(γM−γR)r(t)+(γM−γC)c(t)−γMδ],(4)

where the dynamics of m is omitted, as it can be calculated directly from r and c according to relation (3).

Following the work of [13], the metabolic flux can be described as

vM(t)=eMm(t)=eM(1δ−r(t)−c(t)),(5)

where eM>0 [h−1] represents the effective input from the environment, capturing the concentration and quality of external nutrients. It can be expressed in Michaelis-Menten form as

eM(t):=kMs(t)KM+s(t),

with kM>0 [h−1] and KM>0 [g⋅L−1] denoting the rate and half-saturation constants of metabolism and s(t) the nutrient concentration in the medium. This approximation holds when s(t)≈const or s(t)≫KM. The translational flux is governed by the gene expression machinery and follows a Michaelis-Menten relation:

vR(t)=kRr(t)p(t)KR+p(t),

with kR>0 [h−1] the catalytic rate constant and KR>0 [g⋅L−1] the half-saturation constant. The modulation function gives the effective translation adjusted for the RNA-chaperone activity of CspA,

f(c)=1+ηcnkcn+cn,(6)

where c(t) is the intracellular CspA concentration, kc [g⋅L−1] the activation constant, η a folding-efficiency parameter, and n the Hill coefficient. Here n=2 reflects CspA’s experimentally observed cooperative binding to single-stranded RNA under cold-shock conditions [8]. The effective translational flux incorporating CspA activity is then f(c)vR(t).

Now, let us rescale the model to simplify its mathematical analysis. We introduce the following dimensionless variables and constants:

t^=kRt,p^(t^)=δp(t),r^(t^)=δr(t),m^(t^)=δm(t)c^(t^)=δc(t)EM=eMkR,K=δKR,Kc=δkc,ΓM=γMkR,ΓR=γRkRΓC=γCkR.

Substituting the nondimensional variables into Eq. (2) yields

m^(t^)+r^(t^)+c^(t^)=1.(7)

The system expressed in terms of concentrations (Eq. (4)) can be transformed into a dimensionless form using the rescaled variables and constants:

dp^dt^=EM(1−r^(t)−c^(t))−(p^(t)r^(t)K+p^(t))(1+ηc^2(t)Kc2+c^2(t))−p^(t)[p^(t)r^(t)K+p^(t)(1+ηc^2(t)Kc2+c^2(t))+(ΓM−ΓR)r^(t)+(ΓM−ΓC)c^(t)−ΓM],dr^dt^=α(1−β)(p^(t)r^(t)K+p^(t))(1+ηc^2(t)Kc2+c2(t))−ΓRr^(t)−r^(t)[p^(t)r^(t)K+p^(t)(1+ηc^2(t)Kc2+c^2(t))+(ΓM−ΓR)r^(t)+(ΓM−ΓC)c^(t)−ΓM],dc^dt^=β(p^(t)r^(t)K+p^(t))(1+ηc^2(t)Kc2+c^2(t))−ΓCc^(t)−c^(t)[p^(t)r^(t)K+p^(t)(1+ηc^2(t)Kc2+c^2(t))+(ΓM−ΓR)r^(t)+(ΓM−ΓC)c^(t)−ΓM],(8)

with the initial conditions:

p^(0)=p^0,r^(0)=r^0,c^(0)=c^0,(9)

where

m^(t^)=1−r^(t^)−c^(t^),(10)

and

α(t^)∈[αmin,αmax]=[0,1],β(t^)∈[βmin,βmax]=[0,0.2],∀t^∈[0,T^].

In the scaled system (Eq. (8)), the dynamical expression of m^ is given by Eq. (10). It is also evident that the concentrations m^, r^, and c^ are constrained to the interval [0,1] due to physical limitations imposed by the relation in Eq. (7).

2.1 Parameter Specification

To model bacterial growth and resource allocation under cold-shock conditions, we used kinetic parameters previously measured and validated at 37°C [13]. Because direct experimental data at 15°C are limited, most parameters were adjusted using the Q10 temperature coefficient, a standard approximation for accounting for temperature-dependent reaction rates in microbial physiology:

k(T)=k(Tref)⋅Q10T−Tref10,Q10=2.5,Tref=37∘C, T=15∘C.

This Q10-based scaling is widely used to predict temperature effects on microbial and biological process rates [19]. Parameters derived in this manner preserve the relative ratios among reaction rates while accounting for the overall slowing of metabolic and translational processes under cold-shock conditions. Parameter values associated with CspA-related processes, such as translation efficiency and protein half-life, were taken directly from experimental measurements performed at 15°C in previously published studies. The main parameter values used in the model are summarized in Table 1.

3  Effects of Temperature

Since almost all the parameters of the proposed model (Eq. (4)) depend on the temperature of the environment. Therefore, we analyze the impact of different levels of temperature on the dynamics of the proposed model in this section.

3.1 Steady State Analysis

Following the resource allocation framework of Giordano et al. (2016), we determine the steady state

(p^opt∗,r^opt∗,c^opt∗)∈G

of the system (Eq. (8)), at which the growth rate μ^ attains its maximum. In contrast to approaches where analytical parameterizations of the steady state are available, here the nonlinearities in f(c^) and vR(p^,r^) preclude closed-form solutions. To solve this system numerically, we employ a two-stage procedure. First, a coarse grid search is performed to identify promising regions in the parameter space (α,β). This is followed by a refined optimization using the global, derivative-free differential_evolution algorithm [20] from SciPy. Such numerical strategies are increasingly used in systems biology to overcome convergence issues associated with steady-state constraints [21].

The resulting optimal allocation at the reference temperature 15°C is

(p^opt∗,r^opt∗,c^opt∗,αopt,βopt)=(0.0121,0.0560,0.0056,0.0559,0.0088),

with an associated growth rate

μ^∗=0.122h−1,corresponding to a doubling time Td=ln⁡2μ^∗≈5.69h.

At this optimum, the proteome fractions are distributed as

r^∗:c^∗:m^∗=5.6%:0.6%:93.8%,

where m^∗=1−r^∗−c^∗ denotes the fraction of metabolic proteins. This indicates that only a small investment in CspA is sufficient to sustain growth, while the vast majority of resources remain allocated to metabolic functions under cold-shock conditions.

The convergence of the dynamical system towards the steady state shows that p^(t), r^(t), c^(t), and μ^(t) approach their optimal values for the given initial conditions (Fig. 2). These simulations confirm that the numerically computed optimum corresponds to a stable equilibrium of the system and provide a quantitative baseline for the temperature and control analyses presented in the following sections.

Figure 2: Simulation of the dynamical system showing convergence of p^, r^, c^ and μ^ to their steady-state values at 15∘C. The solid lines show the dynamical behavior of the proposed model, and the dotted lines show the respective steady states.

In addition to the reference case at 15∘C, we computed the steady states for other temperatures (13∘C, 17∘C, and 19∘C) using the same numerical procedure to validate our adopted technique. The resulting temperature-dependent steady-state values of (p^∗,r^∗,c^∗,α∗,β∗,μ^∗) are summarized in Table 2 against each temperature. For each temperature, the associated dynamical simulations (see Figs. 3–5) show that p^(t), r^(t), c^(t), and μ^(t) converges smoothly to their respective steady-state values listed in the Table 2. This behavior confirms that, at every temperature considered, the model admits a unique equilibrium toward which trajectories consistently converge.

Figure 3: Simulation of the dynamical system showing convergence of p^, r^, c^ and μ^ to their steady-state values at 13∘C. The solid lines show the dynamical behavior of the proposed model, and the dotted lines show the respective steady states.

Figure 4: Simulation of the dynamical system showing convergence of p^, r^, c^ and μ^ to their steady-state values at 17∘C. The solid lines show the dynamical behavior of the proposed model, and the dotted lines show the respective steady states.

Figure 5: Simulation of the dynamical system showing convergence of p^, r^, c^ and μ^ to their steady-state values at 19∘C. The solid lines show the dynamical behavior of the proposed model, and the dotted lines show the respective steady states.

At steady state across the analyzed temperature range (13°C–19°C), the model predicts relatively small CspA fractions. This is consistent with growth-optimized proteome allocation, in which CspA acts primarily as a regulatory component rather than a dominant protein species. While CspA can reach high abundance during acute cold shock, such levels reflect transient stress responses following rapid temperature downshift and are not representative of adapted steady-state growth. Under steady conditions, excessive CspA expression imposes translational and metabolic costs, leading cells to preferentially allocate resources toward ribosomal and metabolic proteins, resulting in a minimal-sufficient CspA level.

3.2 Stability Analysis

The local stability of a steady state of a dynamical system can be assessed using the Jacobian matrix approach [22].

Theorem 1: A steady state is locally asymptotically stable if all the eigenvalues of the Jacobian matrix evaluated at the steady state have negative real parts.

For system (8), the Jacobian matrix J is obtained by differentiating the right-hand side with respect to the state variables (p^,r^,c^). The resulting Jacobian matrix is shown below

J=(−Af(c^)(1+(p^2+2p^K))+B1−em−p^(1+p^)K+p^f(c^)−p^(ΓM−ΓR)−em−A1B(1+p^)−p^(ΓM−ΓC)f(c^)A(α(1−β)−r^)f(c^)p^K+p(α(1−β)−2r^)−2r^(ΓM−ΓR)−ΓRA1B(α(1−β)−r^)−(ΓM−ΓC)r^f(c^)A(β−c^)p^K+p^f(c^)(β−c^)−c(ΓM−ΓC)βA1B−ΓC−[c^B+f(c^)]A1−B1−(ΓM−ΓC)).

The Jacobian is evaluated at the steady state (p^∗,r^∗,c^∗). where

A=r^K(K+p^)2,A1=p^r^K+p^,B=ηKc(Kc+c^)2(ΓM−ΓR)r^+(ΓM−ΓC)c^−ΓM.

Because of the strong nonlinearities, we compute the eigenvalues of the Jacobian matrix J numerically at each temperature’s steady state. The results are reported in Table 3.

Across all temperatures examined (13°C–19°C), the Jacobian eigenvalues remain strictly negative, confirming that the steady state at each temperature is locally asymptotically stable. Together with the dynamical simulations (Figs. 2–5), which show that trajectories consistently converge to these steady states, the results demonstrate that the system reliably settles into its temperature-dependent equilibrium. This indicates that the qualitative stability properties of the model are preserved throughout the entire physiological temperature range analyzed.

3.3 Steady-State Growth

We consider E. coli growth processes in the temperature range T0≤T≤T1, where T0=13∘C (286.15 K) and T1=19∘C (292.15 K). One can see how the optimal dimensionless steady-state growth rate μ^opt∗ varies with temperature (Fig. 6). This temperature interval was selected to represent mild cold-acclimation conditions under which E. coli remains metabolically active and cold-shock regulation is functionally relevant. Within this range, temperature effects on metabolic and translational processes can be approximated using continuous scaling relationships such as Q10. The applicability of model predictions outside this interval is limited: as temperature decreases below this range, growth becomes increasingly dominated by physical and survival constraints not captured by the present formulation, while as temperature increases beyond this range, cold-shock regulation progressively loses relevance and standard growth physiology dominates. For each temperature, the steady-state values of the allocation parameters and concentrations (αopt∗,βopt∗,popt∗,ropt∗,copt∗) were determined numerically as those that maximize the growth rate. Because the effective growth time scale changes with temperature, growth rates computed at different temperatures are compared using the same dimensionless time variable t^, chosen here to correspond to 15∘C. All parameter values used in these simulations are taken from Table 1, and temperature variation is applied according to the scaling based on Q10 described in the parameter specification section.

Figure 6: Temperature-dependent behavior of the optimal dimensionless steady-state growth rate μ^∗ for parameters set in Table 1.

Next, we analyzed how bacterial growth depends on the allocation parameter β, which represents the fraction of cellular resources directed toward the synthesis of the cold-shock protein CspA. For each temperature, and for each fixed value of β, the remaining allocation parameter α and the concentrations p, r, and c were numerically optimized at steady state yielding (αopt∗,popt∗,ropt∗,copt∗) to maximize the steady state growth. At low temperatures, moderate allocation to CspA enhances growth by improving translational efficiency through alleviation of inhibitory mRNA secondary structures. However, excessive allocation to CspA diverts translational resources away from ribosome synthesis and metabolic protein production, ultimately reducing the overall growth rate. Across the temperature range considered, the optimal value of β remains the same. Values of β higher than this optimum lead to reduced growth, highlighting the trade-off between stress protection and resource allocation for growth, and indicating that additional allocation to CspA offers no further benefit at higher temperatures and instead diverts resources away from ribosome and metabolic protein synthesis. The growth rate is maximized when β is set to its steady-state optimal value, βopt∗ (Fig. 7).

Figure 7: Effect of the control parameter β on the bacterial growth rate over the temperature range T0=15∘C to T1=18∘C. For each value of β and at each temperature, the allocation parameter α and the concentrations p^, r^, and c^ were numerically optimized at steady state. The growth curve reaches its overall maximum when β takes the steady state optimal value, βopt∗. Increasing β beyond this optimum reduces the growth rate across all temperatures, indicating that CspA contributes positively to growth only up to a certain level.

4  Optimal Control Problem

Optimal control theory provides a principled framework for quantifying how cells should allocate limited resources over time to maximize growth under environmental constraints [13,23–25]. In our setting, the control variables α(t) and β(t) represent the fractions of translational capacity devoted to ribosome synthesis and CspA production. Applying Pontryagin’s Maximum Principle yields the necessary optimality conditions and identifies candidate allocation strategies that balance biomass production with the costs of cold-shock adaptation. We formulate the optimal resource allocation problem for system (Eq. (8)) as the maximization of cumulative biomass production over the finite horizon [0,T^]. Following previous work [13], we assume that E. coli regulates its proteome to maximize structural biomass under prevailing environmental conditions. This leads to the objective functional

J(α(t^),β(t^))=∫0T^μ(t)dt=∫0T^μ^(t^)dt  ⟶  maxα(t^),β(t^),

where

μ^(t^)=μkR=p^(t)r^(t)KR+p^(t^)(1+ηc^2Kc2+c^2)+(ΓM−ΓR)r^(t^)+(ΓM−ΓC)c^(t^)−ΓM,

which represents the cumulative biomass production over the time horizon [0,T^]. The optimization is carried out over all measurable open-loop control functions α(t^) and β(t^).

We assume that environmental variations, such as nutrient upshifts or downshifts, occur as instantaneous shifts [13]. For example, Escherichia coli transitions rapidly between habitats such as the mammalian gut, water bodies, soils, and sediments, where intermediate fluctuations are negligible. Each shift specifies a new value of the environmental input eM, defining a distinct optimal allocation problem. Between two consecutive shifts, eM is treated as constant, representing bacterial adaptation over the finite interval between environmental changes.

Earlier studies formulated growth maximization as an infinite-horizon problem under the overtaking optimality criterion, implicitly assuming that trajectories reach and remain at the steady state corresponding to maximum growth rate. However, the formal existence of such overtaking optimal strategies was not proven, and infinite-horizon problems cannot be solved numerically. Finite-horizon simulations with free terminal states revealed small deviations from steady state near the end of the interval, regarded as artifacts under turnpike theory [26].

Following the approach of Yegorov et al. (2018) [13], we therefore consider two finite-horizon formulations: (i) without a terminal condition, where the final state is free, and (ii) with a fixed terminal condition, constraining the state variables to reach the steady state corresponding to maximum growth:

(p(T^),r(T^),c(T^))=(p^opt∗,r^opt∗,c^opt∗),T^∈(0,+∞).(11)

4.1 Optimality Conditions

In optimal control theory, Pontryagin’s Maximum Principle (PMP) [27,28] provides the first-order necessary conditions for open-loop control strategies. The PMP has been used extensively in microbial resource-allocation models [9,13,22] and is well suited for the present problem, where the controls α(t) and β(t) govern translational allocation to ribosomes and CspA. For our proposed model (Eq. (8)), the application of PMP yields necessary conditions consisting of the Hamiltonian, adjoint system, and optimality conditions, as summarized below. The Hamiltonian is defined as

H(p^,r^,c^,α,β,λ0,λp,λr,λc)= λp[EM(1−r^−c^)−f(c^)(p^r^K+p^)−p^(f(c^)(p^r^K+p^)+(ΓM−ΓR)r^+(ΓM−ΓC)c^−ΓM)]+λr[α(1−β)f(c^)(p^r^K+p^)−ΓRr^−r^(f(c^)(p^r^K+p^)+(ΓM−ΓR)r^+(ΓM−ΓC)c^−ΓM)]+λc[βf(c^)(p^r^K+p^)−ΓCc^−c^(f(c^)(p^r^K+p^)+(ΓM−ΓR)r^+(ΓM−ΓC)c^−ΓM)]+λ0[f(c^)(p^r^K+p^)+(ΓM−ΓR)r^+(ΓM−ΓC)c^−ΓM],(12)

where p^, r^, and c^ denote the non-dimensionalized state variables, λp, λr, and λc are the adjoint variables, and α, β are the time-dependent control variables. From Pontryagin’s Maximum Principle, the adjoint variables satisfy

dλjdt^=−∂H∂yj,j=1,2,3,(13)

which yields the following adjoint equations:

dλpdt^=Kr^K+p^f(c)(λp(1+p^+λr(r^−(1−β)α)+λc(c^−β)+λ0)+λp(p^r^K+p^f(c^)+(ΓM−ΓR)r^+(ΓM−ΓC)c^−ΓM),dλrdt^=λpEM+p^K+p^f(c^)(λp(1+p^+λr(2r^−(1−β)α)+λc(c^−β)+λ0)+(ΓM−ΓR)(λpp^+λr(2r^−1)+λcc^+λ0)+(ΓM−ΓC)λcc^,dλcdt^=λpEM+p^r^K+p^f′(c^)(λp(1+p^+λr(r^−(1−β)α)+λc(c^−β)+λ0)+(ΓM−ΓC)(λpp^+λrr^−+λc(2c^−1)+λ0)+λc(p^r^K+p^f(c^)+(ΓM−ΓC)r^)(14)

where

f(c^)=(1+ηc^2Kc2+c^2),f′(c^)=2ηc^Kc2(Kc2+c^2)2,

λ0≡0  or  λ0≡−1,

(λ0,λp(t^),λr(t^),λc(t^))≠(0,0,0,0),∀t^∈[0,T^].

The Hamiltonian maximum conditions are

Sα(t^):=∂H∂α=λr(t^)(1−β(t^))f(c^)p^(t^)r^(t^)K+p^(t^),Sβ(t^):=∂H∂β=f(c^)p^(t^)r^(t^)K+p^(t^)(λc(t^)−α(t^)λr(t^)).

where the optimal controls α(t) and β(t) are given by

α(t^)={αmin,ifλr(t^)<0,αmax,ifλr(t^)>0,αsing,ifλr(t^)=0,

β(t^)={βmin,ifλc(t^)−α(t^)λr(t^)<0,βmax,ifλc(t^)−α(t^)λr(t^)>0,βsing,ifλc(t^)−α(t^)λr(t^)=0,

which is a necessary condition for an optimal open-loop control. An admissible process

α(t^), β(t^), p^(t^), r^(t^), c^(t^), λ0, λp(t^), λr(t^), λc(t^)

that satisfies PMP conditions is called extremal. It is called normal if λ0<0 and abnormal if λ0=0. Since the controlled system is autonomous, that is, the right-hand sides of the state equations depend on time only through the state and control variables, the Hamiltonian is conserved along any extremal process.

If a switching function Sα(t^) or Sβ(t^) vanishes over some time subinterval, then the corresponding control α(t^) or β(t^) is singular, and the associated trajectory lies on the singular manifold of the system. We now examine the structure of singular arcs arising from these switching functions

Theorem 2 ([23]): Any singular arc us(t) of a normal extremal solution of (OCP) is of at least second order.

Proof: Assume that the process (φ,λ,u) satisfies the Pontryagin Maximum Principle (PMP) as a normal extremal, and set λ0=−1. To describe the singular arcs, suppose the switching functions vanish on a sub-interval τ=[t^1,t^2]⊂[0,T^]. The switching surfaces are defined by

Σα={(φ,λ)∈R2n∣Sα=0},Σβ={(φ,λ)∈R2n∣Sβ=0},

where n is the number of state variables.

On the singular arc for α, Sα=0 implies λr=0 (under the assumptions β≠1 and ϕ(p^,r^,c^)≠0). Similarly, for β, Sβ=0 implies λc=αλr. If both controls are singular, then λr=0 and λc=0. Here we define

ϕ(p^,r^,c^)=f(c^)p^r^K+p^,

as a convenient shorthand for the product term appearing in the switching functions. To compute the derivatives of the switching functions, we use the Poisson bracket operator

{f,g}=∑i=1n(∂f∂λi∂g∂φi−∂f∂φi∂g∂λi),(15)

where φi are the state variables and λi are the costate variables. For simplicity, we write the switching functions as

Sα(t^):=∂H∂α=ϕ(p^,r^,c^)(1−β)λr,Sβ(t^):=∂H∂β=ϕ(p^,r^,c^)(λc−αλr).

The first derivatives are

S˙α=λ˙r(1−β)ϕ(p^,r^,c^)+λr[(1−β)(ϕp^(p^,r^,c^)+ϕr^(p^,r^,c^)+ϕc^(p^,r^,c^))−ϕ(p^,r^,c^)β˙],S˙β=(λ˙c−αλ˙r−α˙λr)(1−β)ϕ(p^,r^,c^)+(λc−αλr)[(1−β)(ϕp^(p^,r^,c^)+ϕr^(p^,r^,c^)+ϕc^(p^,r^,c^))−ϕ(p^,r^,c^)β˙],

where

ϕp^(p^,r^,c^)=∂ϕ(p^,r^,c^)∂p^p^˙,ϕr^(p^,r^,c^)=∂ϕ(p^,r^,c^)∂r^r^˙,ϕc^(p^,r^,c^)=∂ϕ(p^,r^,c^)∂c^c^˙.

Along the singular arc, the first derivatives vanish identically, implying that λ˙r=0 and λ˙c=0. As the Poisson bracket method does not directly resolve this case with two singular controls, we examined the second derivatives S¨α and S¨β by evaluating them on (Sα,S˙α) and (Sβ,S˙β), respectively, using SageMath. The controls α and β do not appear in these second derivatives, confirming that the singular arc is of at least second order.

Now consider α as the sole control while β is fixed. The first derivative is

S˙α=λ˙r(1−β)ϕ(p^,r^,c^)+λr(1−β)(ϕp^(p^,r^,c^)+ϕr(p^,r^,c^)+ϕc(p^,r^,c^)),

which vanishes identically along the singular arc. To compute higher-order derivatives of Sα, we employ the Poisson bracket operator [29]. The Hamiltonian can be expressed as

H=H0+αH1,H1=Sα.

Applying this, we find

H˙1=∂H1∂φφ˙+∂H1∂λλ˙=∑i=1n(∂H∂λi∂H1∂φi−∂H∂φi∂H1∂λi)={H,H1}={H0+αH1,H1}=H01,

since {αH1,H1}=α{H1,H1}=0. The second derivative is then

H¨1=H˙01={H0,H01}+α{H1,H01}=H001+αH101.

To ensure the singular arc is at least of order two, we must show H101=0 along the switching surface Σα. Otherwise, a first-order singular control could be computed as u=−H001/H101. Because Hi depends only on λr,p^,r^,c^, the Poisson bracket becomes

H101={H1,H01}=∂H1∂p^∂H01∂λp−∂H1∂λp∂H01∂p^+∂H1∂r^∂H01∂λr−∂H1∂λr∂H01∂r^+∂H1∂c^∂H01∂λc−∂H1∂λc∂H01∂c^,=μ(p^,r^,c^)∂H01∂λr−λrμr(p^,r^,c^)∂H01∂r^−λrμp(p^,r^,c^)∂H01∂p^−λrμc(p^,r^,c^)∂H01∂c^,

where

∂H01∂r^=λ˙rϕr^(p^,r^,c^)+λr∂∂r^(ϕp^(p^,r^,c^)p^˙+ϕr^(p^,r^,c^)r^˙+ϕc^(p^,r^,c^)c^˙).

Substituting this expression into the previous equation yields

H101=ϕ(p^,r^,c^)[λ˙rϕr^(p^,r^,c^)+λr∂∂r(ϕp^(p^,r^,c^)p˙+ϕr^(p^,r^,c^)r˙+ϕc^(p^,r^,c^)c˙)]−λrϕr^(p^,r^,c^)∂H01∂r^−λrϕp^(p^,r^,c^)∂H01∂p^−λrϕc^(p^,r^,c^)∂H01∂c^,