Component Characterization in a Growth-Dependent Physiological Context: Optimal Experimental Design

Abstract

Synthetic biology design challenges have driven the use of mathematical models to characterize genetic components and to explore complex design spaces. Traditional approaches to characterization have largely ignored the effect of strain and growth conditions on the dynamics of synthetic genetic circuits, and have thus confounded intrinsic features of the circuit components with cell-level context effects. We present a model that distinguishes an activated gene’s intrinsic kinetics from its physiological context. We then demonstrate an optimal experimental design approach to identify dynamic induction experiments for efficient estimation of the component’s intrinsic parameters. Maximally informative experiments are chosen by formulating the design as an optimal control problem; direct multiple-shooting is used to identify the optimum. Our numerical results suggest that the intrinsic parameters of a genetic component can be more accurately estimated using optimal experimental designs, and that the choice of growth rates, sampling schedule, and input profile each play an important role. The proposed approach to coupled component–host modelling can support gene circuit design across a range of physiological conditions.

1. Introduction

Proposed applications of synthetic biology demand complex synthetic constructs involving dynamic internal regulation. Novel analytic and experimental approaches will be needed to efficiently navigate the corresponding design space [1,2]. Model-based design approaches promise to (partially) replace costly experiments with computer simulations. A range of theoretical approaches to automated or computer-assisted design have been published [3,4,5,6,7,8,9,10]. When supported by reliable characterization of genetic regulatory components, automated design algorithms can greatly increase the efficiency of design. As an example, Nielsen et al. demonstrated efficient automated design of very large genetic logic circuits from a carefully designed and characterized regulatory library [11]. Systematic characterization of genetic components can include the generation of standardized data sheets [12], the use of standardized relative units [13], and the predictive characterization of circuit dynamics [14]. This type of characterization is demanding and resource-intensive; to date it has rarely been accomplished. The resulting knowledge deficit is a major bottleneck to the wider use of model-based design.

A drawback of standard approaches to characterization is that they fail to distinguish host physiology from the intrinsic properties of the construct [15]. The physiological state incorporates, among other features, (i) the DNA quantity and gene copy number, (ii) available RNA polymerases (RNAPs) and ribosomes, and (iii) the cell volume, all of which impact gene expression [16,17]. Calibration of model parameters without accounting for the cell’s physiology results in aggregate parameters that describe lumped effects from both host and component. Such a model cannot be trusted to extrapolate beyond the physiological state in which it was calibrated. Separating host state and component behaviour is a difficult and multivariate problem. Experimentalists cannot directly perturb or even measure many physiological properties. The cell’s physiological state can be modulated indirectly by external or internal perturbations including modulation of (i) nutrient sources [17,18], (ii) antibiotics [18], (iii) gene expression burden [18], and (iv) metabolic fluxes [19,20], as well as temperature, pH, and osmolarity. The aggregate effect of each of these perturbations influences growth rate, but predicting host properties from the perturbation or the growth rate is not generally possible. However, for nutrient-limited growth, it has been shown that the exponential growth rate acts as a summary statistic for the physiological state of an E. coli host cell [16,17]. Klumpp and Hwa demonstrated how nutrient-limited growth rates can then be used as an aggregate predictor of the physiological effects on gene expression [16].

In [16], Klumpp and Hwa stop short of proposing an explicit dynamic model for coupling physiology to gene expression. They focus primarily on steady-state behaviour. More recent works have developed coupled models of gene expression, host physiology, and growth rates [21,22,23]. Here, we focus specifically on the use of a coarse-grained model to empirically predict physiological properties from an observed exponential growth rate. Our model distinguishes intrinsic parameters of the genetic construct from extrinsic parameters of the host’s physiological state. This distinction allows for extrapolation across physiological conditions and for reuse of the estimated component parameters. We have aimed to keep the parameter set small to maximize both identifiability and interpretability.

Noise and nonlinearity make precise estimation of the parameters that characterize biomolecular systems a challenging task [24,25]. Optimal experimental design (OED) tools offer a means to improve the efficiency of data collection for model calibration [26,27]. Despite its potential to increase experimental efficiency and the precision of parameter estimates, OED has not seen widespread use in laboratory experiments within systems or synthetic biology. Two notable exceptions are reported by Bandara et al. and Ruess et al.; both groups implemented optimal experimental design for efficient calibration of dynamic biological models [28,29]. OED techniques are especially promising when coupled with optogenetic or microfluidic techniques, which allow for a broad range of dynamic perturbations in vivo [30]. Here, we employ optimal experimental design algorithms originally demonstrated for chemical and bioprocess engineering applications [31,32,33] for the characterization of a genetic component from simulated data.

We develop our physiologically aware model of gene expression for an exponentially growing E. coli population, because this is the only host system for which the necessary data has been collected. We use nutrient quality in exponential phase as a predictable controller of the cell’s growth rate and relevant physiological state. We use the model to optimally design dynamic induction experiments across a set of growth rates to simulate efficient estimation of the intrinsic parameters of the genetic component. We adopt an optimal design approach that expresses the experimental design as an optimal control problem, which can be efficiently solved using numerical optimal control methods, such as multiple shooting. We treat the induction profile, growth rate, and sampling schedule as experimental controls. Their selection is simultaneously optimized over multiple sub-experiments. This simultaneous optimization of sampling schedule and multiple experimental perturbations, including growth rates, is an improvement over previously published accounts of OED in synthetic biology (although see [34]). Previously published analyses either have assumed constant sampling rates or have chosen sampling schedules in a secondary optimization step [29,35], both of which are sub-optimal [36]. We use numerical simulations to demonstrate that the optimal experiments outperform intuitively designed experiments: the optimal designs improve both parameter estimation accuracy and out-of-sample prediction accuracy. We further demonstrate that the use of multiple growth rates is important for model identifiability over realistic parameter ranges.

2. Materials and Methods

In this work, we design optimal experiments for calibrating a dynamic model of expression of a genomically integrated gene, induced by an activating transcription factor (TF). For simplicity, we assume that the controlled induction input is the activating transcription factor copy number, u. (More realistically, the controlled input would be some signal that influences the TF abundance.) Our population-averaged model incorporates both gene-specific intrinsic parameters and parameters that characterize aspects of the host physiology. The physiological parameters of the model are dependent on the steady-state exponential growth rate, $\lambda$, controlled via nutrient quality. The model describes mRNA copy number, $X_{rna}$, and protein copy number, $X_{prot}$,

$$\begin{array}{cc}\hfill \frac{d}{dt}\frac{X_{rna}}{V}& =\alpha \frac{g}{V}\frac{\frac{P_a}{\eta G}{K}_r+\frac{P_a{K}_{rt}}{{\left(\eta G\right)}^2}u\left(t\right)}{1+\frac{P_a}{\eta G}{K}_r+\left(\frac{K_t}{\eta G}+\frac{P_a{K}_{rt}}{{\left(\eta G\right)}^2}\right)u\left(t\right)}-\delta \xi \frac{X_{rna}}{V}\hfill \\ \hfill \frac{d}{dt}\frac{X_{prot}}{V}& =\frac{\beta \frac{R_f}{V}}{K_M+\frac{R_f}{V}}\frac{X_{rna}}{V}-\lambda \frac{X_{prot}}{V}.\hfill \end{array}$$

(1)

The gene’s intrinsic characteristics are captured by the intrinsic parameters $\alpha$, $K_r$, $K_t$, $K_{rt}$, $\delta$, $\beta$, and $K_M$, while expression is also dependent on the physiological parameters: V, g, $P_a$, G, $R_f$, and $\lambda$, which are growth-rate-dependent, and $\eta$ and $\xi$, which are fixed. The intrinsic parameters characterize each gene’s induction behaviour. In contrast, the physiological parameters reflect the state of the host cell (Figure 1). These parameters are summarized in

Table 1. Summary of Intrinsic and Physiological Parameters.

Parameter Label	Intrinsic Parameter	Nominal Value	Feasible Range
Intrinsic Transcription Parameters
Promoter Escape Rate	$\alpha$	20 min−1	[1–30]
RNAP-Promoter Binding	$K_r$	40	[10–40]
TF-Promoter Binding	$K_t$	$5\times {10}^5$	[$2\times {10}^3$–$1\times {10}^6$]
TF-RNAP Interaction	$K_{rt}$	$1.09\times {10}^9$	[$4.02\times {10}^5$–$5.93\times {10}^{10}$]
Intrinsic mRNA Decay Parameters
mRNA Decay Rate	$\delta$	$2.57\times {10}^{-4}$μm−3 min−1	[$7.7\times {10}^{-5}$–$7.7\times {10}^{-4}$]
Intrinsic Translation Parameters
Max. Initiation Rate	$\beta$	$4.0$ min−1	[1–10]
Half-saturating Constant	$K_M$	750 μm−3	[750–1500]
Property Label	Physiological Property	Value at μ = 0.6 db/h	Value at μ = 3 db/h
Physiological Properties of Transcription
Gene Copy Number	g	1.4	5.7
Available RNAP	$P_a$	1000	4000
Genome-lengths of DNA	G	1.3	4.3
Physiological Properties of mRNA Decay
RNase Concentration	$\xi$	900 μm−3 min−1	900 μm−3 min−1
Physiological Properties of Translation
Free Ribosomes	$R_f$	600	7000
General Physiological Properties
Cell Volume	V	0.4 μm−3	2.24 μm−3
Growth Rate	$\lambda$	$0.7\times {10}^{-2}$ min−1	$3.5\times {10}^{-2}$ min−1

, along with nominal values and relevant ranges. Justification for these parameter values is provided in the next section.

2.1. Derivation of the Physiological Gene Expression Model

We develop the model for the case of an exponentially growing E. coli population, where growth rate is controlled by nutrient limitation. While some features of the model may generalize to other cell types, it is only in the case of E. coli that sufficient data has been collected to provide reasonable estimates of functional relations and parameter values.

We employ two standard measures of growth rate. The doubling rate $\mu$ is the inverse of the doubling time ${\tau }_{\mu }$. The exponential growth rate $\lambda =\ln \left(2\right)/{\tau }_{\mu }$ is used to describe population growth as $P_0{e}^{\lambda t}$, and is thus suitable to use as the dilution rate in a differential equation model. We treat the exponential growth rate, $\lambda$, as an independent variable determined by experimental conditions. We also assume the growth rate is constant throughout each experiment, which implies that no nutrient shifts occur.

2.1.1. Cell Volume and Mass, DNA Content, and Protein Mass

Prior work in bacterial physiology has established expressions for the physiological parameters V, G, and g in terms of $\lambda$. (While each of the parameters is expected to undergo random fluctuations and to vary with the cell cycle, these expressions approximate average values over time).

Cell volume has been shown to scale with growth rate exponentially [37,38]:

$$V=V_0{e}^{(C+D)\lambda }.$$

(2)

Here the constant $V_0$ is the “initiation volume” measured to be $V_0=0.28 {\mathsf{\mu }\mathrm{m}}^3$ by Si et al. [37]. The parameters C and D represent the time periods required to replicate the chromosome and to septate, respectively. We take $C=40$ min and $D=20$ min [38].

An expression for the average number of genome equivalent lengths of DNA, G, is given in terms of C, D, and $\lambda$ by Cooper and Helmstetter [38]:

$$G=\frac{1}{\lambda C}(e^{(C+D)\lambda }-e^{D\lambda }).$$

(3)

Individual loci vary in copy number based both on growth rate and on their location relative to the origin of replication. (At faster growth rates there are more copies of genes near the origin because the cell must have multiple rounds of DNA replication underway.) We use the constant $l_{ori}$ to indicate the gene’s position relative to the origin: $l_{ori}=0$ at the origin; $l_{ori}=1$ at the terminus. Bremer and Churchwood [39] provide a simple derivation for the average gene copy number, g, of a specific locus as

$$g=e^{((C+D)-l_{ori}C)\lambda }.$$

(4)

2.1.2. Total RNA Polymerase (RNAP)

We first note that the buoyant density of E. coli has been observed as constant across growth rates: Kubitschek et al. found $\rho =1.09 \mathrm{pg} {\mathsf{\mu }\mathrm{m}}^{-3}$ [40] while Basan et al. report $\rho =0.215 \mathrm{pg} {\mathsf{\mu }\mathrm{m}}^{-3}$ [41] (average of reported values). We use the dry weight data from [17] with Si et al.’s description of volume [37] to select an intermediate density of $\rho =0.55 \mathrm{pg} {\mathsf{\mu }\mathrm{m}}^{-3}$.

The average cell mass for a given growth rate is then

$$M_{Tot}=\rho V=\rho V_0{e}^{(C+D)\lambda }.$$

(5)

This total cell mass $M_{Tot}$ can be partitioned into fractions of protein and other constituents. We fit the protein fraction of the mass, denoted ${\mathsf{\Phi }}_{pr}$, to data from [17] with a linear function (Figure S1):

$${\mathsf{\Phi }}_{pr}={\kappa }_{pr}\lambda +{\mathsf{\Phi }}_{pr0}$$

(6)

with ${\kappa }_{pr}=-6.47$ min and ${\mathsf{\Phi }}_{pr0}=0.65$.

The total RNAP fraction of the overall protein mass, which we denote ${\mathsf{\Phi }}_p$, exhibits an approximately linear dependence on growth rate, which we fit to data from [17] (Figure S2) as

$${\mathsf{\Phi }}_p={\kappa }_p\lambda +{\mathsf{\Phi }}_{p0}.$$

(7)

with ${\kappa }_p=0.30$ min and ${\mathsf{\Phi }}_{p0}=0.0074$. We can then express the total mass of all RNAP protein, $M_{RNAP}$, as

$$M_{RNAP}=\underset{\mathrm{RNAP} \% \mathrm{of} \mathrm{Prot}}{\underbrace{\left({\kappa }_p\lambda +{\mathsf{\Phi }}_{p0}\right)}} \stackrel{\mathrm{Protein} \mathrm{Mass}}{\overbrace{\left({\kappa }_{pr}\lambda +{\mathsf{\Phi }}_{pr0}\right)\underset{\mathrm{Cell} \mathrm{Mass}}{\underbrace{\rho V_0{e}^{(C+D)\lambda }}}}}.$$

(8)

The total number of RNAPs per cell can be determined by dividing $M_{RNAP}$ by the protein mass per RNAP core enzyme, which we denote $m_{rnap}$ and is estimated as $6.3\times {10}^{-7}$ pg [42].

2.1.3. Available RNAP

Each RNAP can be classified by its state: freely diffusing, weakly DNA bound at a non-specific site, actively transcribing other genes, paused or non-functioning during transcription, or immature [43,44,45]. Below, we describe promoter binding in terms of a thermodynamic model that takes into account the observed rapid equilibrium between non-specifically bound and freely diffusing RNAPs [43]. We therefore define the combined pool of free and non-specifically bound RNAPs as the available RNAP pool, with molecular population size $P_a$. Denoting the fraction of RNAP in this available pool by ${\mathsf{\Phi }}_a$, we estimate, from Bakshi et al. [43] and Stracy et al. [46] (details in the Supplementary Materials Section 1.3):

$${\mathsf{\Phi }}_a={\kappa }_a\lambda +{\mathsf{\Phi }}_{a0}$$

(9)

with ${\kappa }_a=-9.3$ min and ${\mathsf{\Phi }}_{a0}=0.59$. Using this relation we have an expression for the available RNAP:

$$P_a=\frac{\rho V_0}{m_{rnap}}\left({\kappa }_a\lambda +{\mathsf{\Phi }}_{a0}\right)\left({\kappa }_p\lambda +{\mathsf{\Phi }}_{p0}\right)\left({\kappa }_{pr}\lambda +{\mathsf{\Phi }}_{pr0}\right)e^{(C+D)\lambda }.$$

(10)

2.1.4. Transcription Rate

We describe the initiation of transcription via a thermodynamic equilibrium model of promoter occupancy involving available RNAPs, transcription factors (TFs), promoter copies, and non-specific binding sites along the genomic DNA [47,48]. At a given growth rate, we assume there are $P_a$ available RNAP copies and $T_a$ active transcription factor copies diffusing along the genomic DNA, and that the DNA contains $N_s$ non-specific binding sites to which these DNA binding proteins may weakly attach and g copies of the regulated promoter of interest. Further, we assume that $N_s\gg P_a,T_a,g$ and that each binding of an RNAP or a TF to a non-specific site or a promoter is characterized by an associated binding energy; ${ϵ}_{rn}$ and ${ϵ}_{rg}$ for RNAP binding to the non-specific sites and promoters, respectively, and ${ϵ}_{tn}$ and ${ϵ}_{tg}$ for transcription factor binding to non-specific sites and promoters, respectively (all $ϵ$ are negative [49]).

We use these species and site counts to enumerate the possible arrangements of RNAP and TF across the genome, and we use the binding energies to derive Boltzmann weights for each arrangement [47]. We denote the differences between the energy involved in binding the promoter and the background non-specific binding as $\Delta {ϵ}_t={ϵ}_{tg}-{ϵ}_{tn}$ and $\Delta {ϵ}_r={ϵ}_{rg}-{ϵ}_{rn}$. (Note, ${ϵ}_{tn}>{ϵ}_{tg}$ and ${ϵ}_{rn}>{ϵ}_{rg}$ so that $\Delta {ϵ}_t$ and $\Delta {ϵ}_r$ are both negative [49].) We denote the Boltzmann weights as $K_r=e^{-\frac{\Delta {ϵ}_r}{k_{B}T}}$, $K_t=e^{-\frac{\Delta {ϵ}_t}{k_{B}T}}$, and $K_{rt}=e^{-\frac{(\Delta {ϵ}_r+\Delta {ϵ}_t+{ϵ}_{pt})}{k_{B}T}}$ (with $k_B$, Boltzmann’s constant, and T, temperature in degrees Kelvin) where ${ϵ}_{rt}$ is the binding energy between RNA polymerase and transcription factor when both are bound to the same promoter. Then, the equilibrium probability of a single promoter being occupied by an RNAP is (further details in the Supplementary Materials Section 1.4)

$$p_{bound}=\frac{\frac{P_a}{N_s}{K}_r+\frac{P_a{T}_a}{N_s^2}{K}_{rt}}{1+\frac{P_a}{N_s}{K}_r+\frac{T_a}{N_s}{K}_t+\frac{P_a{T}_a}{N_s^2}{K}_{rt}}.$$

(11)

With g promoter copies, and presuming that open complex and promoter escape occurs at a fixed rate $\alpha$ [50] we have the initiation rate as

$$\mathrm{Initiation} \mathrm{rate}=\alpha g\frac{\frac{P_a}{N_s}{K}_r+\frac{P_a{T}_a}{N_s^2}{K}_{rt}}{1+\frac{P_a}{N_s}{K}_r+\frac{T_a}{N_s}{K}_t+\frac{P_a{T}_a}{N_s^2}{K}_{rt}}.$$

(12)

We next establish estimates of the parameter values for transcription. To begin, multiplying G by the genomic density of non-specific binding sites, $\eta =5\times {10}^6$ sites/genome, yields an estimate of the total number of non-specific binding sites $N_s$ as a function of growth rate [47].

From Heyduk et al., we have the rate of open complex formation $\alpha$ for the phage lambda $P_R$ promoter as 19.2 min⁻¹ [50]. Assuming $P_R$ exemplifies a strong promoter, we estimate a reasonable range based on analysis of constitutive promoters and mutants of the $P_R$ sequences as $1-30$ min⁻¹ [50,51,52,53].

The constants $K_r$ and $K_t$ can be interpreted as ratios of dissociation constants for the DNA-binding species (RNAP and TF) binding to non-specific DNA versus the promoter sequence [47]. This provides a convenient method for constraining their feasible values from reported dissociation rates. Expressing them as such yields; $K_r=K_{rnap}^{ns}/K_{rnap}^{prom}$, $K_t=K_{tf}^{ns}/K_{tf}^{prom}$ and $K_{rt}=K_r{K}_t\exp ({ϵ}_{rt}/k_{B}T)$. Here $K_{rnap}^{ns}$ and $K_{tf}^{ns}$ are dissociation constants of the TF and RNAP with respect to non-specific DNA binding and $K_{rnap}^{prom}$ and $K_{tf}^{prom}$ are dissociation constants for promoter binding. The non-specific dissociation constant for RNAP, $K_{rnap}^{ns}$, has been observed to be approximately $\mathrm{10,000}$ nM [47]. The promoter-specific dissociation constant, $K_{rnap}^{prom}$, varies from promoter to promoter. For lacP1, it is approximately 550 nM and for T7 it is approximately 3 nM [47]. To represent a relatively weak constitutive leak for an inducible promoter, we expect values nearer to the lacP1 dissociation constant would be reasonable, and so presume $K_{rnap}^{prom}$ could range from 250 to 1000 nM. Using the above value for $K_{rnap}^{ns}$ and the range for $K_{rnap}^{prom}$, we estimate a feasible range for $K_r$ to be between 10 and 40 (unitless). We set $K_{rnap}^{prom}$ to be 250 nM, so our nominal value for $K_r$ is 40, suggesting a slightly stronger leak than found in uninduced lacP1. (Over the growth rates we consider, we found $K_r$ to be practically unidentifiable, with high variability in its estimated value. We therefore fixed it to the nominal value for our experimental design and fitting.)

To estimate $K_t$, Stormo suggests a reasonable range for the promoter-specific dissociation constant, $K_{tf}^{prom}$, of 0.01–1000 nM [54]. We assume we are working with a relatively strong activator and that $K_{tf}^{prom}$ lies in the range 0.01–5 nM. Stormo also suggests that non-specific binding dissociation constants are between three and six orders of magnitude less than the specific binding constants. For simplicity, we follow Bintu et al. and let $K_{tf}^{ns}=$ 10,000 nM [47]. This yields a range for $K_t$ between 2000 and ${10}^6$. We chose $K_{tf}^{prom}=0.02$ and $K_t=5\times {10}^5$ as nominal values.

As noted, the parameter $K_{rt}$ can be expressed as $K_{rt}=K_r{K}_t\exp (-{ϵ}_{rt}/k_{B}T)$ with ${ϵ}_{rt}$ representing the energy involved in the RNAP-TF interaction at the promoter. Few estimates exist for such values. Bintu et al. use demonstrative values of ${ϵ}_{rt}/k_{B}T$ ranging from $-3.5$ to $-4.5$. We allow a wider range from $-3$ to $-5$. This allows $K_{rt}$ to range from $4.02\times {10}^5$ to $1.61\times {10}^{10}$. We take a nominal value of $K_{rt}=1.09\times {10}^9$.

2.1.5. mRNA Degradation

The data in [17] indicate that the mRNA decay rate is linear in mRNA copy number, with a growth-rate-independent rate constant. In [16], the authors hypothesize that this is due to the maintenance of a constant concentration of RNase E, the primary ribonuclease involved in mRNA decay initiation in E. coli [55]. (RNase E exhibits auto-regulation which appears to keep its concentration constant [56,57,58]. Moreover, RNase E appears to be expressed in excess, resulting in insensitivity to small changes in its abundance [58,59].) Under this assumption and using mass action, we have a simple model for mRNA decay as follows

$$\mathrm{mRNA} \mathrm{Decay} \mathrm{Rate} (\mathrm{copy} \# \mathrm{per} \min )=\delta \frac{E}{V}{X}_{rna}.$$

(13)

Here E is the copy number of RNase E. Assuming $\frac{E}{V}$ is constant across growth rates (and therefore $E\propto V$), we can express the mRNA degradation rate as

$$\mathrm{mRNA} \mathrm{Decay} \mathrm{Rate} (\mathrm{copy} \# \mathrm{per} \min )=\delta \xi X_{rna}.$$

(14)

Here $\xi$ is the constant concentration of relevant degrading enzymes in the host.

The parameter $\delta$ represents the susceptibility of the the mRNA transcript to RNase degradation, and acts as a mass action constant. The half-life of mRNAs can range from 1–10 min; Chen et al. suggest the mean RNA half life is near to $2.5$ min [60,61,62]. Based on a constant concentration of RNase E of $\xi =900$ μm⁻³ [63] and taking a nominal half-life of 3 min, we suggest $\delta \approx 2.57\times {10}^{-4}$ μm⁻³ min⁻¹.

2.1.6. Total and Free Ribosome Populations

A linear relation for the fraction, ${\mathsf{\Phi }}_r$, of protein mass that is composed of ribosomal protein was derived in [18]:

$${\mathsf{\Phi }}_r={\kappa }_r\lambda +{\mathsf{\Phi }}_{r0}.$$

(15)

Fitting this model to data from Bremer [17] yields estimates of ${\kappa }_r=5.48$ min and ${\mathsf{\Phi }}_{r0}=0.030$, (see Figure S4). The total ribosomal protein mass $M_{Rib}$ is then

$$M_{Rib}=\left({\kappa }_r\lambda +{\mathsf{\Phi }}_{r0}\right)\left({\kappa }_{pr}\lambda +{\mathsf{\Phi }}_{pr0}\right)\rho V_0{e}^{(C+D)\lambda }.$$

(16)

Each ribosome has a mass of $2.7$ MDa of which 35% is protein [64]. The individual ribosomal protein mass is therefore $m_{rib}=1.57\times {10}^{-6}$ pg. We then have the number of ribosomes per cell, $R_{Tot}$, as

$$R_{Tot}=\left({\kappa }_r\lambda +{\mathsf{\Phi }}_{r0}\right)\left({\kappa }_{pr}\lambda +{\mathsf{\Phi }}_{pr0}\right)\frac{\rho V_0{e}^{(C+D)\lambda }}{m_{rib}}.$$

(17)

Based on results in Dai et al. [65], we have assumed that about 10% of the ribosomes are inactive. Assuming that these are free, we denote this fraction as ${\mathsf{\Phi }}_f=0.1$ and have

$$R_f={\mathsf{\Phi }}_f{R}_{Tot}=\frac{{\mathsf{\Phi }}_f\rho V_0}{m_{rib}}\left({\kappa }_r\lambda +{\mathsf{\Phi }}_{r0}\right)\left({\kappa }_{pr}\lambda +{\mathsf{\Phi }}_{pr0}\right)e^{(C+D)\lambda }.$$

(18)

2.1.7. Translation Rate

We employ the Michaelis–Menten model of translation initiation proposed by Borkowski et al. [66]. In terms of the free ribosome concentration $R_f$,

$$\mathrm{Translation} \mathrm{Rate} (\mathrm{copy} \# \mathrm{per} \min )=\frac{\beta \frac{R_f}{V}}{K_M+\frac{R_f}{V}}{X}_{rna}.$$

(19)

In this model, the mRNA’s ribosome binding site (RBS) is characterized by two constants: $\beta$, the maximal translation initiation rate per mRNA, and $K_M$, a half-saturating constant specific to the given RBS. A detailed justification for this model is provided in the Supplementary Materials Section 1.6.

To estimate $\beta$, we note that translation initiation rates on the order of $\beta =4 {\min }^{-1}$ have been observed for lacA [67]. Additional studies of RBS activity suggest a wide variation in expression levels induced by RBS alterations, ranging over orders of magnitude [68,69]. We presume that a reasonable range for the maximal translation initiation rate is from 1 min⁻¹ to 10 min⁻¹.

The values of $K_M$ reported by Borkowski et al. [66] are reported in arbitrary units. They observe a range of saturation levels from near linear behaviour to near constant saturation [66] from which we can propose a range of $K_M$ values based on the observed free ribosome concentration, which in our model is between 1000 μm⁻³ and 3000 μm⁻³. We let $K_M$ range between 750 μm⁻³ and 1500 μm⁻³, with a nominal value of $K_M=750$ μm⁻³, which is approaching saturation (and hence the near constant translation efficiency observed by Klumpp and Bremer [16,17]).

2.2. Optimal Experimental Design

To illustrate calibration of the intrinsic parameters of the model, we consider a set of dynamic induction experiments conducted over multiple growth rates, where the experimenter can select (i) the (constant) growth rate, (ii) the time-varying induction profile, and (iii) the sampling schedule. We define an experiment as a set of $N=3$ sub-experiments, each with a constant growth rate $\lambda$ (possibly repeated), organized into vector $\mathit{\lambda }=[{\lambda }^{\left(1\right)},{\lambda }^{\left(2\right)},{\lambda }^{\left(3\right)}]$. (The exponential growth rates, ${\lambda }^{\left(i\right)}$, are reported as doubling rates ${\mu }^{\left(i\right)}$ (db/hr) in the results for interpretability). In each sub-experiment, the population begins at rest and responds to a dynamic induction signal in the form of a time varying transcription factor copy number, $\mathit{u}\left(t\right)=[u^{\left(1\right)}\left(t\right),u^{\left(2\right)}\left(t\right),u^{\left(3\right)}\left(t\right)]$. (Such an input $u^{\left(i\right)}\left(t\right)$ could be implemented by, e.g., a calibrated induction system [70] or closed loop control of fluorescently tagged TFs [71].) For computational efficiency, we have restricted each $u^{\left(i\right)}\left(t\right)$ to be a piecewise constant, with 6 constant control values delivered for 100 min each, for a total duration of $t_f=600$ min. We have constrained the maximum value of each $u^{\left(i\right)}\left(t\right)$ to $u_{max}=1000$ and the minimum to $u_{min}=0$. This wide range was selected to ensure it spans the unsaturated range of the promoter well in excess. For each sub-experiment, we also select sampling schedules for both mRNA and protein species. A single sampling schedule for a specific species is defined as a list of points, $s_{\left(species\right)}=\{p_1,\dots p_l,\dots ,p_L\}$. Each $p_l=({\tau }_l,c_l)$ consists of a time ${\tau }_l\in [0,600]$ and a positive integer count $c_l$ of the number of samples taken at that time. We have restricted the design such that $c_l\in \{1,2,3\}$. The total number of samples for each species, $c_{Tot}$ is constrained to be no more than 12. The sampling schedule for each sub-experiment, $S^{\left(i\right)}$, consists of a schedule for each species $S^{\left(i\right)}=\{s_{\left(rna\right)}^{\left(i\right)},s_{\left(prot\right)}^{\left(i\right)}\}$. These are collected into an overall schedule: $\mathit{S}=\{S^{\left(1\right)},S^{\left(2\right)},S^{\left(3\right)}\}$.

Optimization over the sampling schedules as defined above involves integer programming, which poses significant computational challenges. We avoid these by employing a relaxation approach [33,72,73]. We treat the sampling schedule as a continuous sampling density such that $s_{\left(species\right)}^{\left(i\right)}$ corresponds to a density $w_{\left(species\right)}^{\left(i\right)}\left(t\right)$. Each sub-experiment has sampling density $W^{\left(i\right)}\left(t\right)=\{w_{\left(rna\right)}^{\left(i\right)}\left(t\right),w_{\left(prot\right)}^{\left(i\right)}\left(t\right)\}$ and the overall sampling schedule is $\mathit{W}\left(t\right)=[W^{\left(1\right)},W^{\left(2\right)},W^{\left(3\right)}]$. These sampling densities are restricted to be a piecewise constant over 48 equal intervals (each of length $12.5$ min) during each sub-experiment. Each sampling density $w_{\left(species\right)}^{\left(i\right)}\left(t\right)$ is upper bounded by $w_{max}=\frac{3}{12.5}$, and the integral of the sampling distribution, ${\widehat{w}}_n^{\left(species\right)}$ is bounded by the maximal number of samples, $c_{Tot}$:

$${\widehat{w}}_{\left(species\right)}^{\left(i\right)}={\int }_0^{t_f}{w}_{\left(species\right)}^{\left(i\right)}\left(t\right)dt\le c_{Tot}.$$

(20)

The relaxed schedule must be discretized to arrive at an implementable sampling schedule [33]. In practice, the optimal densities are often bang-bang [73] and so can be easily discretized by the sampling interval length (12.5 min) to recover integer sample counts. To account for non-integer values, we applied the Sum-Up-Rounding strategy, a common heuristic in integer programming, to recover integer solutions that approximate the optimal discrete schedule [72]. This discretization strategy effectively rounds the continuous density while respecting the sampling constraints [72]. After generating the sampling counts, $c_l$, from the sampling densities, we fixed their associated times ${\tau }_l$ to the centre of each of the 12.5 min intervals.

To characterize a genetic construct with the induction experiments as defined above, we seek an experimental design that can accurately estimate the intrinsic parameter set $\theta$. As stated in Section 2.1.4, the parameter $K_r$ proved to be practically unidentifiable under the model specifications described above. Consequently, we fixed it to its nominal value, and took the intrinsic parameter set to be estimated as $\mathit{\theta }=[\alpha \beta K_t{K}_{rt}\delta ]$. These parameters characterize the regulated promoter, its regulating transcription factor and the downstream gene sequence, independent of growth context. We denote an estimate of the true parameters as $\widehat{\mathit{\theta }}$. For our objective, we seek to minimize the determinant of the covariance matrix $\det (\mathrm{cov}(\widehat{\mathit{\theta }}))$, in what is known as a D-optimal design [74]. The determinant of the covariance matrix is also known as the generalized variance [75]; minimizing it is equivalent to minimizing the volume of the confidence ellipsoid defined by the covariance matrix [76]. Because the model is nonlinear and the true parameter vector is, by definition, unknown, it is not possible to estimate the parameter covariance matrix a priori. Instead of estimating $\mathrm{cov}(\widehat{\mathit{\theta }})$ directly, we use the Fisher information computed at an initial guess ${\mathit{\theta }}_{\mathit{o}}$ as a proxy. In linear models (or in the limit of large sample sizes or low signal to noise ratio), the inverse of the FIM is asymptotically equivalent to $\mathrm{cov}(\widehat{\mathit{\theta }})$ [77]:

$$\mathcal{I}{(\mathit{\theta },t_f)}^{-1}\approx \mathrm{cov}(\widehat{\mathit{\theta }}).$$

(21)

In this case minimizing $\det (\mathrm{cov}(\widehat{\mathit{\theta }}))$ is equivalent to maximizing the determinant of the Fisher information matrix, $\det \left(\mathcal{I}\right)$. In the non-linear finite sample regime in which our experiments are conducted, this relation is admittedly tenuous; the use of an initial guess ${\mathit{\theta }}_{\mathit{o}}$ also introduces potential errors. However, the process of iterating the local approximation and successively updating ${\mathit{\theta }}_{\mathit{o}}$ after each experiment has been numerically demonstrated to yield convergence to the true parameter set in some cases [26,27]. We thus define our objective as ${\mathsf{\Theta }}_D\left(\mathcal{I}\right)=\det \left(\mathcal{I}\right)$, which is generally known as the D-optimality score [76]. (Other options such as A, E, and E-modified optimality minimize other properties of the covariance [74,76,78].)

To determine the optimal set of experiments, we follow an optimal control-based procedure described in [31,32,33]. (This approach to OED was demonstrated by application to chemical and bioprocess models. It has seen limited use in systems biology [28] and has not previously been applied to component characterization in a synthetic biology context.) To formulate a control problem we must state the objective, the D-optimal score ${\mathsf{\Theta }}_D\left(\mathcal{I}\right)=\det \left(\mathcal{I}\right)$, in terms of the model dynamics. We restate the model in the following generic form:

$$\frac{d\mathit{X}}{dt}=\mathit{F}(\mathit{X},\mathit{\theta },{\lambda }^{\left(i\right)},u^{\left(i\right)}\left(t\right)).$$

(22)

Here $\mathit{X}=[X_{rna},X_{prot}]$ and $\mathit{F}$ is the right hand side of the model expressed in Equation (1).

We determine the local sensitivities of $X_{rna}$ and $X_{prot}$ with respect to each parameter ${\theta }_i$, by solving the following system of sensitivity equations

$$\frac{d}{dt}\frac{\partial \mathit{X}}{\partial {\theta }_j}=\frac{\partial \mathit{F}}{\partial {\theta }_j}+\frac{\partial \mathit{F}}{\partial \mathit{X}}\frac{\partial \mathit{X}}{\partial {\theta }_j}.$$

(23)

Because the parameters are dimensioned, their scalings can vary widely, leading to poor conditioning of the Fisher information matrix and associated computations [31]. To rectify this, we scale the sensitivities by the parameter values (i.e., we use logarithmic sensitivities)

$${\overline{\mathit{X}}}_{{\theta }_j}=\frac{\partial \mathit{X}}{\partial \log \left({\theta }_j\right)}={\theta }_i\frac{\partial \mathit{X}}{\partial {\theta }_j}.$$

(24)

Applying this change of variables yields

$$\frac{d{\overline{\mathit{X}}}_{{\theta }_j}}{dt}={\theta }_j\frac{\partial \mathit{F}}{\partial {\theta }_j}+\frac{\partial \mathit{F}}{\partial \mathit{X}}{\overline{\mathit{X}}}_{{\theta }_j}.$$

(25)

The initial conditions for the state variables and their sensitivities are constrained to be at steady state with respect to a zero induction input ($u=0$):

$$\mathit{X}(t=0)={\mathit{X}}_{SS}(\mathit{\theta },\lambda ,u=0)$$

(26)

$$\begin{array}{cc}\hfill {\overline{\mathit{X}}}_{{\theta }_j}(t=0)& ={\overline{\mathit{X}}}_{{\theta }_i,SS}(\mathit{\theta },\lambda ,u=0).\hfill \end{array}$$

(27)

Because the steady state depends on both the growth rate and the initialized parameter values, these initial conditions are not constants. They therefore appear as nonlinear constraints in the optimization problem.

We denote the sampling variance for observations of species by ${\sigma }_{species}^2$. We assume normally distributed errors with variances equal to 5% of the species quantity:

$${\sigma }_{rna}^2=\left(0.05\right)X_{rna}$$

(28)

$${\sigma }_{prot}^2=\left(0.05\right)X_{prot}.$$

(29)

We further assume no covariance between species measurements. The (scaled) Fisher information matrix $\mathcal{I}(\mathit{\theta },t)$ can then be written for the relaxed problem as a differential equation, as per [33]

$$\frac{d}{dt}{\mathcal{I}}_{jk}(\mathit{\theta },t)={\overline{X}}_{{\theta }_j}\mathsf{\Omega }\left(t\right){\mathsf{\Sigma }}^{-1}\left(t\right){\overline{X}}_{{\theta }_k}{\mathcal{I}}_{jk}(\mathit{\theta },0)=0$$

(30)

where the matrices $\mathsf{\Omega }\left(t\right)$ and $\mathsf{\Sigma }\left(t\right)$ are defined as

$$\mathsf{\Omega }\left(t\right)=\left[\begin{array}{cc}{w}_n^{\left(rna\right)}\left(t\right)& 0\\ 0& w_n^{\left(prot\right)}\left(t\right)\end{array}\right],\mathsf{\Sigma }\left(t\right)=\left[\begin{array}{cc}{\sigma }_{rna}^2{X}_{rna}\left(t\right)& 0\\ 0& {\sigma }_{prot}^2{X}_{prot}\left(t\right)\end{array}\right].$$

(31)

Here $w_n^{\left(species\right)}\left(t\right)$ are sampling densities, as defined above. Further details on the integration of the Fisher information over continuous sampling densities can be found in [33,73]. Because samples are assumed to be independent across time points as well as species, the Fisher information for the experiment is additive across both. (Here we use the Fisher information for homoskedastic error, despite the heteroskedasticity of the model. We found this simpler formulation performed equivalently and was more tractable.) Therefore, we can express the cumulative Fisher information for a given sub-experiment from time 0 to $t_f$ as

$${\mathcal{I}}_{jk}^{\left(i\right)}(\mathit{\theta },t_f)={\int }_0^{t_f}{\overline{X}}_{{\theta }_i}^{\left(j\right)}{\mathsf{\Omega }}^{\left(i\right)}\left(t\right){\mathsf{\Sigma }}^{-1}\left(t\right){\overline{X}}_{{\theta }_k}^{\left(i\right)}dt.$$

(32)

The Fisher information for the total experiment, ${\mathcal{I}}_{Tot}$, is then the sum over all sub-experiments

$${\mathcal{I}}_{Tot}=\sum _{i=1}^N{\mathcal{I}}^{\left(i\right)}(\mathit{\theta },t_f).$$

(33)

The overall objective (D-optimality score) is then

$${\mathsf{\Theta }}_D\left({\mathcal{I}}_{Tot}\right)=\det \left(\sum _{i=1}^N{\mathcal{I}}^{\left(i\right)}\right).$$

(34)

The matrix is symmetric, so only the diagonal and lower triangular elements need to be computed.

The value of ${\mathsf{\Theta }}_D\left({\mathcal{I}}_{Tot}\right)$ can vary over many orders of magnitudes for different experiments, and so to improve numerical accuracy we take the logarithm of the determinant as the objective. Finally, because most optimization packages are minimizers, we invert the sign:

$$\begin{array}{c}\hfill \underset{\mathit{u},\mathit{\lambda },\mathit{W}}{\min }-\ln \left({\mathsf{\Theta }}_D\left({\mathcal{I}}_{Tot}\right)\right).\end{array}$$

(35)

For numerical stability, we compute the determinant using QR factorization (details in Section 2 of the Supplementary Materials).

In summary, we formulate our OED optimal control problem as

$$\begin{array}{c}\hfill \begin{array}{c}\hfill Objective:\\ \hfill \underset{\mathit{u},\mathit{\lambda },\mathit{W}}{\min }-\ln \left({\mathsf{\Theta }}_D\left({\mathcal{I}}_{Tot}\right)\right)& =-\ln \left(\det \left(\sum _{i=1}^N{\mathcal{I}}^{\left(i\right)}(\mathit{\theta },t_f)\right)\right)\hfill \\ \hfill Subjectto(\forall i\in \{1,\dots ,N\left\}\right):\\ \hfill \frac{d{\mathit{X}}^{\left(i\right)}}{dt}& =\mathit{F}({\mathit{X}}^{\left(i\right)},\mathit{\theta },{\lambda }^{\left(i\right)},u^{\left(i\right)}\left(t\right))\hfill \\ \hfill \frac{d{\overline{\mathit{X}}}_{{\theta }_j}^{\left(i\right)}}{dt}& ={\theta }_j\frac{\partial \mathit{F}({\mathit{X}}^{\left(i\right)},\mathit{\theta },{\lambda }^{\left(i\right)},u^{\left(i\right)}\left(t\right))}{\partial {\theta }_j}+\frac{\partial \mathit{F}({\mathit{X}}^{\left(i\right)},\mathit{\theta },{\lambda }^{\left(i\right)},u^{\left(i\right)}\left(t\right))}{\partial \mathit{X}}{\overline{\mathit{X}}}_{{\theta }_j}^{\left(i\right)},\forall {\theta }_i\in \mathit{\theta }\hfill \\ \hfill \frac{d{\widehat{w}}_{sp}^{\left(i\right)}}{dt}& =w_{sp}^{\left(i\right)}\left(t\right),\forall sp\in \{rna,prot\}\hfill \\ \hfill \frac{d{\mathcal{I}}_{jk}^{\left(i\right)}}{dt}& ={\overline{X}}_{{\theta }_j}^{\left(i\right)}{\mathsf{\Omega }}^{\left(i\right)}\left(t\right){\left({\mathsf{\Sigma }}^{\left(i\right)}\left(t\right)\right)}^{-1}{\overline{X}}_{{\theta }_k}^{\left(i\right)},\forall {\theta }_j,{\theta }_k\in \mathit{\theta }\hfill \\ \hfill With Constraints:\\ \hfill {\mathit{X}}^{\left(i\right)}(t=0)& ={\mathit{X}}_{SS}^{\left(i\right)}(\mathit{\theta },{\lambda }^{\left(i\right)},u=0)\hfill \\ \hfill {\overline{\mathit{X}}}_{{\theta }_j}^{\left(i\right)}(t=0)& ={\overline{\mathit{X}}}_{{\theta }_j,SS}^{\left(i\right)}(\mathit{\theta },{\lambda }^{\left(i\right)},u=0),\forall {\theta }_j\in \mathit{\theta }\hfill \\ \hfill {\mathcal{I}}_{jk}^{\left(i\right)}(\mathit{\theta },0)& =0\forall {\theta }_j,{\theta }_k\in \mathit{\theta }\hfill \\ \hfill {\widehat{w}}_{sp}^{\left(i\right)}\left(0\right)& =0,sp\in \{rna,prot\}\hfill \\ \hfill {\widehat{w}}_{sp}^{\left(i\right)}\left(t_f\right)& <c_{Tot},sp\in \{rna,prot\}\hfill \\ \hfill 0& <w_{sp}^{\left(i\right)}\left(t\right)<w_{max},sp\in \{rna,prot\}\hfill \\ \hfill 0& \le u^{\left(i\right)}\left(t\right)\le u_{max}\hfill \\ \hfill {\lambda }_{min}& <{\lambda }^{\left(i\right)}<{\lambda }_{max}.\hfill \end{array}\end{array}$$

(36)

This problem formulation matches that used by Telen et al. and Hoang et al. [31,33]. We used a multiple shooting algorithm to solve the optimal control problem [32]. We implemented this algorithm in CasADi, a rapid prototyping optimal control toolbox, using its MATLAB interface [79]. CasADi uses algorithmic differentiation to compute first and second derivatives with respect to both the objective and the constraints. This higher-order information is used for optimization in an interior-point barrier method implemented in the nonlinear programming package IPOPT [80]. Further details of the multiple shooting algorithm and optimization settings can be found in Section 2 of the Supplementary Materials.

3. Results

3.1. Comparing Lumped and Physiologically Aware Models

Models of gene expression typically use lumped parameters that confound physiological effects with intrinsic features. Such models are usually calibrated against data collected from cells grown in a particular medium (and so at a single growth rate). To illustrate the consequences of neglecting growth effects, we consider a lumped version of our growth dependent model (Equation (1))

$$\begin{array}{cc}\hfill \frac{d\left[X_{rna}\right]}{dt}& =A\frac{B+Cu}{1+B+\left(D+C\right)u}-E\left[X_{rna}\right]\hfill \\ \hfill \frac{d\left[X_{prot}\right]}{dt}& =F\left[X_{rna}\right]-\lambda \left[X_{prot}\right]\hfill \end{array}$$

(37)

where parameters $A=\frac{\alpha g}{V}$, $B=\frac{P_g{K}_r}{\eta G}$, $C=\frac{P_a{K}_{rt}}{{\left(\eta G\right)}^2}$, $D=\frac{K_t}{\eta G}$, $E=\delta \xi$, and $F=\frac{\beta R_f}{V(K_{M}V+R_f)}$ are treated as growth-rate-independent constants. (The one exception is the protein decay rate, which is equal to the exponential growth rate.) Of course, this parameter lumping does not impact model behaviour for the growth rate at which the model is calibrated, but accurate extrapolation to other growth rates cannot be assured. This loss of accuracy is illustrated in Figure 2. Panels A and B show predictions of the full growth-dependent model and the lumped model, both calibrated to fast growth rates (3 db/h). As the growth rate drops, significant deviation in the predicted steady state copy number of both protein and mRNA is observed. In Figure 2C,D, we see similar deviation for the steady state concentrations of these species. (We note that the relative deviations in concentration are comparatively small. Protein concentration is important for, e.g., metabolic enzymes. In contrast, copy number is more important for DNA binding proteins, which tend to disperse along the DNA rather than the total cell volume [43,71].) A complementary measure of inaccuracy due to lumping appears in Figure 2E, which shows how the lumped parameter values vary with growth rate when the full model is used to determine their values (in comparison to their values calibrated at the reference fast growth rate of 3 db/h). The divergence in these parameters sets has a significant effect on the model’s dynamic behaviour. Figure 2F shows the root mean squared error (RMSE) between the two models’ response to the input profile depicted in Figure 2G. The difference in dynamic behaviour between the lumped and full model under this dynamic induction, at a slow growth rate of 0.5 db/h, is shown in Figure 2G.

3.2. Null and Optimal Experimental Designs

We next explore how experimental design can improve the estimation of the intrinsic parameters of the full growth-dependent model (Equation (1)). We have selected a null experiment, shown in

Table 2. Null (non-optimal) experimental designs.

Experiment	Growth Rates (db/h)	Induction Pattern (log10)	Sampling Schedule	Optimality
Null Experiment	{0.6,1.8,3}			−63.8
Growth Variant	{2,2.5,3}			−61.5
Sampling Variant	{0.6,1.8,3}			−62.6
Induction Variant	{0.6,1.8,3}			−60.5

, to represent a reasonable non-optimized design. It consists of a set of logarithmically (base 10) distributed pulses, delivered at evenly spaced growth rates, with evenly spaced samples taken every half-hour (mRNA and protein samples are taken at the same time). We propose this as a sensible first experiment for fitting a dynamic model. We designed the null experiments assuming limited information on the dynamics of the specific gene expression system. They were selected to involve a reasonably comprehensive set of perturbations and measurements. For comparison, we examine three variants of the null experiment, also shown in Table 2, each with a perturbation to either the growth rates, induction profile, or sampling rate of the null. The growth variant is identical to the null case, except it is performed over a narrower range of growth rates. The sampling variant differs from the null by redistributing the samples so that their rate is halved but the sampling number at each point is doubled (with the same total number of samples). The induction variant uses linearly (rather than logarithmically) spaced strengths of induction pulses. Shown with each design is its predicted optimality score (the negative log determinant of the scaled Fisher information matrix computed at the true parameter set: Equation (35)). As expected, each variant provides less information than the null (as measured by D-optimality).

In contrast, Figure 3 depicts an optimal experiment for our model. This design was generated by using the true nominal parameter vector as the algorithm’s initialization parameters. This is an artificial scenario (the true parameters are not known initially, by definition), but it serves to demonstrate the difference between intuitive designs and optimal designs selected by our method. In Figure 3, each column describes a sub-experiment, labeled with the corresponding optimal growth rate. The top row, (A–C), depicts the optimal input profiles; the middle row, (D–F), shows the system response in both mRNA and protein copy number. The last row, (G–I), shows the sampling densities (in the shaded regions) as well as the results of the Sum-Up-Rounding discretization scheme (depicted by the stem plots, where the height is an integer representing the sample count). The optimality score of this design is −69.6. (In this case, the optimal design selected the maximum growth rate twice. We suspect this is due to higher sensitivities at faster growth. More complex models typically demand observations over a wider range of growth rates).

Comparing the null and optimal designs we see that, while the null experiments have an intuitive appeal because of their wide, even distribution of experimental choices, none achieves a score similar to the optimized design. Further, while null variants exhibit differences in optimality scores, when compared to the optimized experiment these differences are small. This suggests that manually tuning aggregate design measures is less effective than the holistic optimization provided by the OED approach.

3.3. Utility of Optimal Designs for Parameter Identification and Prediction

The difference in optimality scores between the null and optimized experiments suggests significant improvements in parameter estimation using the optimized design. However, the statistical theory used to derive the optimality score can only be guaranteed to hold a priori for linear models in the limit of small observation variances and large samples [77]. To validate that our optimality scores correspond to improved parameter estimation accuracy, we used simulated experiments to assess the correlation between theoretical and observed variability. To generate simulated data, we simulated the model using the nominal true parameters and added normally distributed observation error with a variance equal to 5% of the corresponding species count. We used a multiple shooting approach for parameter estimation; obvious outliers were removed (details in Section 2 of the Supplementary Materials). Figure 4A shows, on the vertical axis, the logarithm of the generalized variance (determinant of the observed covariance matrix) for the collection of parameter estimates (each corresponding to an experimental design). This measure of fit variability is computed from the collection of parameter fits to independent simulations of the given experiment. From this set of estimates we computed the observed covariance matrix and the generalized variance. The optimality score of the design is shown on the horizontal axis (lower is better), computed at the true parameter value. Figure 4A shows that the optimality score (objective of the OED approach) and the observed parameter covariance (sampled measure of design ’quality’) correlate well. In particular, we note that the optimal design achieves both a better optimal score and a smaller generalized parameter variance when compared to the non-optimal designs. This suggests that the objective function (using the homeostatic FIM) is a useful measure of design quality, despite the nonlinearity and heteroskedasticity of the model. We also note that adopting a uniform sampling schedule or reducing the range of growth rates is expected to result in considerably worse performance, as evidenced by comparison between null variants.

In addition to accuracy of parameter estimates, accurate predictions of the system behaviour for out-of-sample conditions is also important for model-based design. For linear regression models, D-optimal designs that minimize the generalized variance of the parameter estimate are equivalent to designs that minimize the upper bound on prediction variance (General Equivalence Theorem) [81]. This guarantee does not generalize to our dynamic, nonlinear, and heteroskedastic model. To verify if prediction accuracy indeed correlates with D-optimality for our model, we ran a second simulation study. For each of the parameter estimates used in calculating the fit variability for the experimental designs, we simulated the model response for a new dynamic experiment. We chose out-of-sample conditions: growth rates and induction levels not used in fitting. We simulated the model with the true nominal parameter values, and computed the integrated squared error (ISE) between the true and fitted values in each case (thus including error at all time points). Figure 4B shows the generalized parameter variance along the horizontal and the log ISE on the out-of-sample data set on the vertical axis. The plots show that the generalized parameter variance correlates reasonably well with prediction accuracy: the optimal design performs better than any of the null designs. This reflects the linear case, in which the D-optimal design sets an upper-bound on the expected prediction variance [81].

So far, we have compared designs in the idealized (and artificial) scenario in which OED is applied to the model parametrized by its true values. We relaxed that assumption by first generating five intrinsic parameter sets from a uniform distribution spanning the intervals specified in the feasible ranges from Table 1. We then ran the OED algorithm using models parametrized by these ‘perturbed’ parameter sets. The optimality scores of these designed experiments, when evaluated against simulated data generated by the true parameter set, are plotted in Figure 4A,B, depicted as X’s.

4. Discussion

We have proposed a physiologically aware model of gene expression in E. coli that accounts for the effects of nutrient limitation on the host physiology. The model is more complex than standard models of gene expression, but this comes with several benefits. The model naturally suggests a partitioning of the parameter set into (1) intrinsic parameters that characterize the genetic component and (2) a set of empirically derived parameters characterizing the host’s physiological state as a function of nutrient-mediated growth rate. The intrinsic parameters can thus be reused across a range of growth conditions. In particular, because the intrinsic parameters can be linked to sequence properties of the component (e.g., promoter affinities, RBS strength, and mRNA stability), they could provide insight as to which aspect of a component’s DNA sequence could be altered to achieve a desired effect across a variety of contexts. Future work could attempt to link such intrinsic parameters to sequence properties directly. In contrast, the extrinsic parameters can be used to predict the host’s context based on the observed growth rate in a range of nutrient conditions.

Accurate estimation of the intrinsic parameter set is more difficult than estimation in context-naive models. We have demonstrated that optimal experimental design can mitigate this challenge by identifying maximally informative experiments. We applied a comprehensive OED platform, optimizing over constant growth rate, time-dependent induction profile, and sampling schedule. Many current experimental design approaches in systems biology use general purpose optimization algorithms that scale poorly to such multivariate and non-linear designs [30,74]. This often leads to non-optimal selection of certain design variables [36]. To achieve computational efficiency in our optimization tasks, we re-interpreted the problem via optimal control, with growth rates, sampling schedules, and induction levels as control variables. This optimal control framing of OED problems has so far been underutilized in systems biology; it allows access to the rich tool-set developed by control and process engineers, including direct multiple-shooting and collocation methods [31,32,33]. These methods are ideally suited for use with emerging experimental tools such as optogenetic induction systems and automated culture and microfluidic devices [30]. Our numerical results suggest that the multiple shooting algorithm provides an efficient method to optimize experiments, improving both parameter estimation accuracy and prediction accuracy over unoptimized designs.

The optimal experimental design algorithm presented here could be improved in a number of ways. We used the simplest and most tractable form of the Fisher information. More accurate formulations account for the sample variance’s parametric dependence [29] and for the constraints imposed by initial conditions [82]. Additionally, our approach provides a design that may be only locally optimal. Current iterative designs are based on the assumption that the iterated algorithm does not become trapped in some local minimum. Past numerical studies of iterative applications of OED have shown consistent convergence to the true parameter set (global optimum), but there is scope for more rigorous investigation [26,27]. Improved computational efficiency may allow users to address multiple scenarios and larger models, or even implement online experimental design algorithms that can update the design in real-time [83,84]. Our algorithm was implemented with CasADi, which provides rapid-prototyping capabilities for control problems, algorithmic differentiation and interfaces to powerful non-linear programming solvers like IPOPT [79,80]. This tool facilitates rapid implementation, but requires some mathematical expertise. Further packaging of OED tools for use by experimentalists would no doubt increase adoption in systems and synthetic biology.

The experiments proposed in this work are demanding, requiring sampling of multiple species and control of inputs over extended time periods. In practice, it is not currently possible to precisely control the transcription factor count in vivo. It would suffice to map the experimental input (e.g., light, chemical inducer) to the expected average copy number of TF per cell. This could be implemented through precise calibration experiments, achievable with current techniques [70]. Future experiments may be able to implement closed-loop control using real-time measurements to ensure the robust tracking of desired inputs. It also may be possible to modify the OED algorithm to account for the effects of input variability (an errors-in model). In this work, we have assumed time-series measurements of both mRNA and protein quantities. Emerging experimental tools, combined with assays like RNA-seq [85], will likely enable these complex dynamic experiments in the future [30].

For the model considered here, observing only one species results in a structurally unidentifiable parameter set. Any parameter lumping to alleviate unidentifiability leads to a loss of modularity in characterization. Even with observations of both species, the parameter $K_r$, which characterizes promoter leak, had to be excluded from the analysis because it was practically unidentifiable over reasonable parameter ranges. In practice, it may be possible to ignore promoter leakiness in many cases. However, for strong leaks, it is important to identify the parameter $K_r$; fixing it, as we have, may introduce significant bias in other parameter estimates. In addition to these specific identifiability concerns, the re-usability of parameter estimates is limited by the use of relative units, which are specific to the measurement instrument they are calibrated on. The methods used in this work should ideally be implemented using absolute units, which will allow for comparison between parameter values calibrated on different instruments and in different labs [86,87,88].

Our description of physiological state was restricted to a single case: exponentially growing E. coli host cells, with growth rate determined by nutrient quality. Our model formulation was made possible by significant experimental efforts into characterizing precisely this physiological response [17]. Beyond nutrient limitation, several other relevant growth perturbations are of interest to synthetic biologists, including translation-inhibiting antibiotics, gene expression burden, and metabolic knockdowns. Phenomenological growth laws and coarse-grained proteome models have been extended to some of these conditions [18,19,41]. Recent results suggest that certain metabolic perturbations may effect proteome partitioning (and possibly RNAP and ribosomal fractions) in a manner similar to nutrient limitation [19,20]. Certain physiological properties, such as the ribosome-to-protein ratio, can also be predictably linked to the growth rate controlled by expression burden or antibiotic dosage [18,41]. The linear ribosome-to-protein ratios may even generalize across species (Figure S1 of [18]). The details necessary to link these coarse-grained models to the specifics of gene expression, as originally proposed in Klumpp and Hwa [16] (and further elaborated on here) are still lacking. While the nutrient-limited growth theories have taken decades to build, modern techniques can likely accelerate this host characterization process, allowing generalization across strains, growth inhibition conditions, and potentially even microbial species.

Extensions of phenomenological growth theories to heterologous protein expression burden and modification of metabolic fluxes could have significant impact on metabolic engineering. These effects are especially relevant to the emerging sub-discipline of dynamic metabolic engineering, where synthetic gene expression circuits dynamically modulate both fluxes and heterologous enzyme expression [89,90,91]. Design of such systems is challenging; the engineered regulatory components modulate expression burden and enzymatic activity, which in turn affects the growth rate, broader cell physiology, and the regulatory component behaviour itself [92,93,94]. An extended growth theory combined with the optimal experimental design algorithm and gene expression models presented here could be valuable for guiding future work in this area.

We have proposed the coupling of physiologically aware modelling and OED techniques to address limitations in the accuracy and generalizability of component characterization. Current gene expression models often fail to account for the host or environmental state, and are calibrated with poorly constrained parameter estimates using ad hoc experimental designs. These shortcomings contribute to the limited use of model-based design in synthetic biology. Poor parameter estimates result in lackluster predictive accuracy, and context-naive model predictions cannot generalize, resulting in recalibration for each new design and circumstance, and thus little advantage of model-based design over trial-and-error approaches. Our model has yet to be validated in vivo, but wider use of such coarse-grained, context-aware models and optimal experimental designs will ideally maximize researchers’ return on experimental investment and aid efforts to rationally design gene expression circuits.