Team:NUS Singapore/Model

NUS iGEM 2019

index

Overview Preliminary Simulations Curve-fitting Sensitivity Analysis and Predictions Learnings & Future Work References

Overview

Modelling techniques empower Synthetic Biology by offering a better understanding of biological systems and making useful predictions. Hence, Modelling served as an integral part of the Design aspect of our DBTL cycle. This led to a close interplay between Modelling and Experiments as the models were trained using experimental data and in turn made predictions to guide our experimental design.

Our approach to modelling our growth switch and growth knob control systems involved the following steps:

a) Preliminary Simulations
b) Curve-fitting
c) Sensitivity Analysis and Predictions

We have successfully modelled bacterial cell growth under the influence of our two control systems:

i) the HicA-HicB toxin-antitoxin system (growth switch) - which has the ability to switch the bacterial cells on and off.
ii) the SgrS glucose uptake system (growth knob) - which fine-tunes bacterial productivity.

Modelling these systems allowed us to predict their behaviour and understand how they could be optimised before incorporating them into the final design for demonstration.

PRELIMINARY SIMULATIONS

Our first step in modelling the control systems involved carrying out a literature review to understand the mechanisms at play. We developed ordinary differential equations that could mathematically describe cellular processes occurring in the systems. Using our preliminary models, we simulated the likely behaviour of the systems by writing a MATLAB script. These simulations provided answers to questions that the Wet Lab team had in the initial stages of planning their experiments.

Kinetics of HicA and HicB

The foundation of our growth switch lies in the HicA-HicB system which consists of the HicA toxin and HicB antitoxin that can interact in different ways to regulate growth. Our preliminary model for the HicA-HicB system was based on a mechanism proposed in literature (Winter et al., 2018) which included the synthesis, degradation and interaction of the toxin-antitoxin proteins (Fig. 1). Our control system was built to have inducible transcription of the toxin and antitoxin, whose kinetics were modelled by Hill equations. The subsequent processes of translation and interaction between the proteins were described by simple rate equations with the appropriate reaction orders. For example, it has been proposed that two HicA monomers bind with a HicB tetramer to form a complex, therefore the rate equation for complex formation will be second order with respect to HicA and first order with respect to HicB tetramer.

Growth regulation by the HicA-HicB system

The cell growth was modelled by a logistic equation as this model can capture the growth trends that we usually observe for E.coli. The growth rate tends to slow down with a growing population size which may be due to factors like lack of resources or accumulation of toxic substances. To account for the inhibitory effect of toxin on growth, a Hill equation term for repression was incorporated into the logistic equation.

In the absence of literature and experimental data on the HicA-HicB system, we made educated guesses about the parameter values based on a previously-established model of another toxin-antitoxin system known as HipA-HipB (Rotem et al., 2010).

Fig. 1: Proposed mechanism of the HicA-HicB toxin-antitoxin system.

Assumptions

1) All inducible and repressible systems were assumed to follow Hill equation characteristics.
2) Initial concentration of all species inside the cell was assumed to be 0.
3) The initial value of \(\mathit{OD}_{600}\) was assumed to be 0.2.
4) The production of HicA toxin was assumed to be unaffected by its effect on the global translation (Fig. 1).
5) The degradation of the HicA-HicB complex was considered negligible compared to its dissociation.

Equations of the Preliminary Model

Kinetic equations

\[\frac{\mathrm{d}[\mathit{HicAmRNA}]}{\mathrm{dt}} = \mathit{syn_{HicAmRNA}} \times \Bigg(\frac{[\mathit{IPTG}]}{k_{\mathit{IPTG}} + [\mathit{IPTG}]}\Bigg) - \mathit{deg_{HicAmRNA}} \times [\mathit{HicAmRNA}] \] \[\frac{\mathrm{d}[\mathit{HicBmRNA}]}{\mathrm{dt}} = \mathit{syn_{HicBmRNA}} \times \Bigg(\frac{[\mathit{ara}]}{k_{\mathit{ara}} + [\mathit{ara}]}\Bigg) - \mathit{deg_{HicBmRNA}} \times [\mathit{HicBmRNA}] \] \[\frac{\mathrm{d}[\mathit{HicA}]}{\mathrm{dt}} = \mathit{syn_{HicA}} \times \mathit{HicAmRNA} - \mathit{deg_{HicA}} \times [\mathit{HicA}] - 2k_a[\mathit{HicA}]^2[\mathit{HicB}] + 2k_d \times [\mathit{complex}] \] \[\frac{\mathrm{d}[\mathit{HicB}]}{\mathrm{dt}} = \frac{1}{4}\mathit{syn_{HicB}} \times \mathit{HicBmRNA} - \mathit{deg_{HicB}} \times [\mathit{HicB}] - k_a[\mathit{HicA}]^2[\mathit{HicB}] + k_d \times [\mathit{complex}] \] \[\frac{\mathrm{d}[\mathit{complex}]}{\mathrm{dt}} = k_a \times [\mathit{HicA}]^2 \times [\mathit{HicB}]^2 - k_d \times [\mathit{complex}] \]

Growth equations

\[\frac{\mathrm{d}\mathit{OD}_{600}}{\mathrm{dt}} = \mathit{OD}_{600} \times \mu_{\mathrm{max}} \times \Bigg(1 - \frac{\mathit{OD}_{600}}{\mathit{OD_{600_{\mathrm{max}}}} \times \Big(\frac{k_{\mathit{HicA}}}{k_{\mathit{HicA}} + [\mathit{HicA}]}\Big)} \Bigg) \]

Symbol	Description
\([\mathit{HicAmRNA}]\)	Concentration of HicAmRNA present in the cell
\(\mathit{syn_{HicAmRNA}}\)	Rate of transcription of the HicA gene, associated with the p_Lac promoter
\([\mathit{IPTG}]\)	Concentration of the inducer (IPTG) used to induce the transcription of HicAmRNA
\(k_{\mathit{IPTG}}\)	Concentration of IPTG required to achieve half of the maximum synthesis rate of HicAmRNA
\(\mathit{deg_{HicAmRNA}}\)	Rate of degradation associated with HicAmRNA
\([\mathit{HicBmRNA}]\)	Concentration of HicBmRNA present in the cell
\(\mathit{syn_{HicBmRNA}}\)	Rate of transcription of the HicB gene, associated with the p_BAD promoter
\([\mathit{Ara}]\)	Concentration of the inducer (Ara) used to induce the transcription HicBmRNA
\(k_{\mathit{Ara}}\)	Concentration of Arabinose required to achieve half of the maximum synthesis rate of HicBmRNA
\(\mathit{deg_{HicBmRNA}}\)	Rate of degradation associated with the HicBmRNA
\([\mathit{HicA}]\)	Concentration of HicA toxin protein present in the cell
\(\mathit{syn_{HicA}}\)	Rate of translation of HicAmRNA, associated with rbs34
\(\mathit{deg_{HicA}}\)	Rate of degradation associated with the HicA protein
\([\mathit{HicB}]\)	Concentration of HicB antitoxin protein present in the cell
\(\mathit{syn_{HicB}}\)	Rate of translation of HicBmRNA, associated with rbs34
\(\mathit{deg_{HicB}}\)	Rate of degradation associated with the HicAB protein
\([\mathit{complex}]\)	Concentration of HicA-HicB complex present in the cell
\(k_a\)	Association constant for toxin-antitoxin complex formation
\(k_d\)	Dissociation constant for toxin-antitoxin complex breakdown
\(\mu_{\mathrm{max}}\)	Maximum growth rate of the E.coli strain in the absence of induction
\(\mathit{OD}_{600_{\mathrm{max}}}\)	Final OD₆₀₀ value achieved by the E.coli strain in the absence of induction
\(k_{HicA}\)	Concentration of HicA toxin required to achieve half of the maximum repression of growth rate

Simulations of the preliminary model

Our preliminary models were used to simulate bacterial growth when different inducer concentrations were used to trigger the production of toxin and antitoxin (Figs. 2 and 3). A knowledge of these trends helped shape the initial experiments for characterising the HicA-HicB toxin-antitoxin system.

Fig. 2: Simulations of bacterial growth with varying inducer (IPTG) concentrations used to induce HicA toxin.

Fig. 3: Simulations of bacterial growth with varying inducer (Ara) concentrations used to induce HicB antitoxin.

Kinetics of SgrS and ptsG

Our growth knob focuses on the manipulation of cell glucose uptake by SgrS to achieve finely-tuned growth. In this control system, SgrS RNA binds and degrades the glucose transporter mRNA ptsG, hence regulating glucose consumption and growth of the cell (Fig. 4). To capture the effect of SgrS regulation of ptsG on cellular growth, our kinetic equations describing the endogenous transcription and degradation of SgrS and ptsG, and the association, dissociation and degradation of the SgrS-ptsG complex were constructed with close reference to a paper that modelled this endogenous mechanism (Fei et al., 2015).

Growth regulation by the SgrS system

The growth and glucose consumption of the cell was modelled by the Monod equation, which is usually used to describe microbial growth rates in an environment with limited nutrients. The effect of ptsG levels on glucose uptake of the cell is modelled by a Hill equation, linking the kinetics of the SgrS system to cellular growth.

Fig. 4: Illustration of the mechanism of action of SgrS.

Assumptions

1) The kinetic steps listed above are described with rate equations, while the inducible transcription of SgrS in our system is modelled by a Hill equation.
2) The limiting substrate for growth in our system has been assumed to be glucose.
3) The initial \(\mathit{OD}_{600}\) is assumed to be 0.1, the initial glucose concentration to be 0.2%, the initial SgrS concentration to be 0.1636 \(\mu\mathrm{M}\), the initial ptsG mRNA concentration to be 4.632 \(\mathrm{nM}\), and the initial SgrS-ptsG complex concentration to be 0.2525 \(\mathrm{nM}\).

Equations

Kinetic equations

\[ \mathit{syn_{SgrS_{\mathrm{ind}}}} = \mathit{syn_{SgrS_{\mathrm{max}}}} \times \Bigg(\frac{[\mathit{aTc}]}{[\mathit{aTc}] + k_{\mathit{aTc}}}\Bigg) \] \[ \frac{\mathrm{d}[\mathit{SgrS}]}{\mathrm{dt}} = \mathit{syn_{SgrS_{\mathrm{endo}}}} + \mathit{syn_{SgrS_{\mathrm{ind}}}} - \mathit{deg_{SgrS}} \times \mathit{SgrS} - k_{\mathrm{on}} \times \mathit{SgrS} \times {\mathit{ptsG}} + k_{\mathrm{off}} \times [\mathit{complex}] \] \[ \frac{\mathrm{d}[\mathit{ptsG}]}{\mathrm{dt}} = \mathit{syn_{ptsG}} - \mathit{deg_{ptsG}} \times \mathit{ptsG} - k_{\mathrm{on}} \times \mathit{SgrS} \times {\mathit{ptsG}} + k_{\mathrm{off}} \times [\mathit{complex}] \] \[ \frac{\mathrm{d}[\mathit{complex}]}{\mathrm{dt}} = k_{\mathrm{on}} \times \mathit{SgrS} \times \mathit{ptsG} - k_{\mathrm{off}} \times \mathit{SgrS_{ptsG}} - k_{\mathrm{deg}} \times [\mathit{complex}] \]

Growth equations

\[ k_G = k_{G_{\mathrm{max}}} \times \Big(1 + \frac{k_p}{\mathit{ptsG}}\Big) \] \[ \frac{\mathrm{d}\mathit{OD}_{600}}{\mathrm{dt}} = \mathit{OD}_{600} \times \mu_{max} \times \frac{[\mathit{glucose}]}{[\mathit{glucose}] + k_G} \] \[ \frac{\mathrm{d}[\mathit{glucose}]}{\mathrm{dt}} = -\frac{\mathrm{d}\mathit{OD}_{600}}{\mathrm{dt}} \div \mathit{Yield} \]

Symbol	Description
\([\mathit{SgrS}]\)	Concentration of SgrS RNA present in the cell
\([\mathit{ptsG}]\)	Concentration of ptsG mRNA present in the cell
\([\mathit{complex}]\)	Concentration of SgrS-ptsG complex present in the cell
\(\mathit{syn_{SgrS_{\mathrm{endo}}}}\)	Endogenous rate of transcription of the SgrS gene
\(\mathit{syn_{SgrS_{\mathrm{ind}}}}\)	Induced rate of transcription of the SgrS gene, associated with the p_tet promoter
\(\mathit{syn_{SgrS_{\mathrm{max}}}}\)	Maximum rate of transcription of the SgrS gene, associated with the p_tet promoter
\(k_{\mathit{aTc}}\)	Concentration of aTc required to achieve half of the maximum synthesis rate of Sgrs
\(\mathit{deg_{SgrS}}\)	Endogenous degradation rate of SgrS
\(\mathit{syn_{ptsG}}\)	Endogenous transcription rate of ptsG
\(\mathit{deg_{ptsG}}\)	Endogenous degradation rate of ptsG
\(k_{\mathrm{on}}\)	Association constant for SgrS-ptsG complex formation
\(k_{\mathrm{off}}\)	Dissociation constant for SgrS-ptsG complex breakdown
\(k_{\mathrm{deg}}\)	tRNase E-mediated degradation constant of SgrS-ptsG complex
\([\mathit{aTc}]\)	Concentration of the inducer (aTc) used to induce the transcription Sgrs RNA
\(\mathit{OD}_{600}\)	OD₆₀₀
\([\mathit{glucose}]\)	Concentration of glucose in the growth media
\(\mu_{\mathrm{max}}\)	Maximum rate of cell growth
\(k_G\)	Concentration of glucose required to achieve half of the maximum of cell growth rate
\(k_{G_{\mathrm{max}}}\)	Maximum of \(k_G\)
\(k_p\)	Concentration of ptsG required to achieve half of the repressive effect of ptsG on \(k_{G_{\mathrm{max}}}\)
\(\mathit{Yield}\)	Yield coefficient, biomass per mass of glucose utilized

Simulations of the preliminary model

In our preliminary model, we could simulate reduced bacteria growth within a range of inducer (aTc) concentrations (Fig. 5). This information provided the wet lab team with an initial direction in the experimental characterization of the SgrS glucose uptake system.

Fig. 5: Simulations of bacterial growth with varying inducer (aTc) concentrations used to induce SgrS.

CURVE-FITTING

After developing the preliminary models, our next step was curve-fitting which allowed us to verify if our models could reproduce the experimental cell growth well. At the same time, this would allow us to estimate the values for the parameters in the model. The curve-fitting of our models to the experimental data was also carried out by writing a MATLAB script. Through numerous iterations of the DBTL cycle, our preliminary model was continuously refined in order to better determine the behaviour of the system and make more accurate predictions. After multiple rounds of developing equations and fitting to experimental data, we established our final models which could best represent the behaviour of our systems. We prioritized goodness-of-fit and parsimony in choosing our final models. Based on the curve-fitting results, we then verified that our final models could emulate the trends observed in the experiments and we obtained the optimised values of the equation parameters.

Refining the preliminary model

Although our preliminary model described in the previous section gave a reasonable fit to the experimental data (Fig. 6), we recognised that it lacked reliability. This was not only because of insufficient data to validate such a complex model with multiple intermediate steps, but also due to the very few constraints on the parameter values. Furthermore, the results of the curve-fitting showed that the preliminary model could not represent the experimental data well enough. The model simulations were not able to reproduce the observed growth resumption after the induction of antitoxin. This implied that including an equation for the sequestration of toxin by the antitoxin through complex formation is not effective in describing the neutralizing effect of the antitoxin. The curve-fitting results also show that the model underestimates the delay in observing the growth arrest after toxin induction. Thus, we opted for a simpler model that could capture the key trends observed in the experimental data while being biologically reasonable.

Our final model for the HicA-HicB system assumes protein synthesis to be a single-step process that can be directly induced, thus neglecting the intermediate steps of transcription and translation. Instead of modelling the interaction between the toxin-antitoxin as a reaction, the antitoxin was considered to have a repressive effect on the toxin which has been modelled by a form of the Hill equation. The curve-fitting results of the final model showed that the new set of equations could effectively emulate key features such as the delay in growth arrest and resumption upon the induction of toxin and antitoxin respectively as well as the dose response of growth to different amounts of toxin (Fig. 7) and antitoxin (Fig. 8).

Fig. 6: Fitting of the preliminiary model to the HicA-HicB characterisation data, for varying inducer (ara) concentrations used to induce HicB antitoxin while a constant inducer (IPTG) concentration was used to induce HicA toxin.

Fig. 7: Fitting of the final model to the HicA-HicB characterisation data, for varying inducer (IPTG) concentrations used to induce HicA toxin.

Fig. 8: Fitting of the final model to the HicA-HicB characterisation data, for varying inducer (ara) concentrations used to induce HicB antitoxin while a constant inducer (IPTG) concentration was used to induce HicA toxin.

Assumptions

1) All inducible and repressible systems were assumed to follow Michaelis-Menten characteristics.
2) Initial concentration of all species inside the cell was assumed to be 0.
3) The initial value of \(\mathit{OD}_{600}\) was assumed to be 0.2
4) The production of HicA toxin was assumed to be unaffected by its effect on the global translation (Fig. 1).
5) The effect of HicB antitoxin on the HicA toxin has been assumed to be limited to a repression of the toxin's inhibtiory effect on growth.

Equations of the Final Model

\[\frac{\mathrm{d}[\mathit{HicA}]}{\mathrm{dt}} = \mathit{syn_{HicA}} \times \Bigg(\frac{[\mathit{IPTG}]}{[\mathit{IPTG}] + k_{\mathit{IPTG}}}\Bigg) - \mathit{deg_{HicA}} \times [\mathit{HicA}] \] \[\frac{\mathrm{d}[\mathit{HicB}]}{\mathrm{dt}} = \mathit{syn_{HicB}} \times \Bigg(\frac{[\mathit{ara}]}{[\mathit{ara}] + k_{\mathit{ara}}}\Bigg) - \mathit{deg_{HicB}} \times [\mathit{HicB}] \] \begin{align} & \frac{\mathrm{d}\mathit{OD}_{600}}{\mathrm{dt}} = (\mathit{OD}_{600} + \mathit{OD}_{600_{\mathrm{min}}}) \times \mu_{\mathrm{max}} \times \Bigg(1 - \Big(\frac{\mathit{OD}}{\mathit{OD}_{600_{\mathrm{max}}}}\Big)^l\Bigg) \times \Bigg(\frac{{k_{\mathit{HicA}}}^{nHic}}{{k_{\mathit{HicA}}}^{nHic} + {\mathit{HicA}_{\mathit{eff}}}^{nHic}}\Bigg) \\ & \text{where } \mathit{HicA_{eff}} = \mathit{HicA} \times \Big(\frac{k_{\mathit{HicB}}}{k_{\mathit{HicB}} + \mathit{HicB}}\Big) \end{align}

Parameter	Description	Optimized Value (3 s.f.)
\(\mathit{syn_{HicA}}\)	Maximum synthesis rate of HicA toxin	\(7.32 \times 10^{-8} \text{ } \mu\mathrm{M}\mathrm{min}^{-1}\)
\(\mathit{deg_{HicA}}\)	Degradation rate of HicA toxin	\(4.00 \times 10^ {-3} \text{ } \mathrm{min}^{-1}\)
\(k_{\mathit{IPTG}}\)	Concentration of IPTG required to achieve half of the maximum synthesis rate of HicA	\(0.0949 \text{ } \mathrm{mM}\)
\(\mathit{syn_{HicB}}\)	Maximum synthesis rate of HicB antitoxin	\(1.16 \times 10^{-7} \text{ } \mu\mathrm{M}\mathrm{min}^{-1}\)
\(\mathit{deg_{HicB}}\)	Degradation rate of HicB antitoxin	\(0.0130 \text{ } \mathrm{min}^{-1}\)
\(k_{\mathit{ara}}\)	Concentration of Arabinose required to achieve half of the maximum synthesis rate of HicB	\(0.0700\%\)
\(\mathit{OD}_{600_{\mathrm{min}}}\)	Parameter accounting for the inflection point in the growth curve	\(0.485\)
\(\mathit{OD}_{600_{\mathrm{max}}}\)	Final OD₆₀₀ value achieved by the E.coli strain in the absence of induction	\(1.35\)
\(\mu_{\mathit{max}}\)	Maximum growth rate of the E.coli strain in the absence of induction	\(0.0250 \text{ } \mathrm{min}^{-1}\)
\(l\)	Logistic coefficient accounting for gradient of the growth curve	\(0.165\)
\(k_{\mathit{HicA}}\)	Concentration of “effective” HicA toxin required to achieve half of the maximum repression of growth rate	\(1.35 \times 10^{-6} \text{ } \mu\mathrm{M}\)
\(n_{\mathit{Hic}}\)	Hill-coefficient for the repressive effect of “effective” HicA toxin on growth rate	\(6.06\)
\(k_{\mathit{HicB}}\)	Concentration of HicB antitoxin required to make the “effective” HicA toxin level half of the original HicA toxin concentration	\(1.34 \times 10^{-6} \text{ } \mu\mathrm{M}\)

Refining the preliminary model

Our preliminary model was not able to capture the key trends observed in the experimental data (Fig. 9). In the preliminary model, the repressible effect of SgrS was shown to be saturated beyond an inducer concentration of 10nM. However, further growth reduction was observed in the experimental data, when inducer concentrations were increased from 10nM to 1500nM. Additionally, we were unable to carry out experiments to validate most of the parameters in the kinetic steps of the model. Hence, we opted for a simpler model still based on the Monod equation, with a reduced number of unknown parameters. In our final model, a direct repressible effect of the inducer on cell growth is assumed, ignoring the intermediary kinetics of SgrS and ptsG. Our final model demonstrated a much better fit to the experimental data (Fig. 10). It could describe the key trends in the experimental data, such as the initial slowing down of growth in the first 1-3 hours, the growth reduction trend in the exponential phase, and the increase to the same final \(\mathit{OD}_{600}\) value over time.

Fig. 9: Fitting of the preliminiary model to the SgrS characterisation data, for varying inducer (aTc) concentrations used to induce SgrS.

Fig. 10: Fitting of the final model to the SgrS characterisation data, for varying inducer (aTc) concentrations used to induce SgrS.

Assumptions

1) The effect of the inducer on growth is modelled by a Hill equation.
2) The dilution and degradation of the inducer are accounted for in a rate equation.
3) The initial glucose concentration in the media is assumed to be 0.2% as per our experimental design, while the initial \(\mathit{OD}_{600}\) is assumed to be 0.14.

Equations of the Final Model

\[\mu = {\mu}_{\mathrm{max}} \times \frac{[\mathit{glucose}]}{[\mathit{glucose}] + k_G} \times \frac{{k_{aTc}^m}}{{[\mathit{aTc}]}^m + k_{aTc}^m} \] \[\frac{\mathrm{d}[\mathit{aTc}]}{\mathrm{dt}} = -(\mathit{deg_{aTc}} + \mu) \times [\mathit{aTc}] \] \[\frac{\mathrm{d}\mathit{OD}_{600}}{\mathrm{dt}} = \mathit{OD}_{600} \times \mu \] \[\frac{\mathrm{d}[\mathit{glucose}]}{\mathrm{dt}} = -\frac{\mathrm{d}\mathit{OD}_{600}}{\mathrm{dt}} \div \mathit{Yield} \]

Parameter	Description	Optimized Value (3 s.f.)
\(\mathit{deg_{aTc}}\)	Degradation rate of \([\mathit{aTc}]\)	\(2.15 \times 10^{-4} \text{ }\mathrm{hr}^{-1}\)
\(\mu_{\mathrm{max}}\)	Maximum growth rate of the E.coli strain in the absence of induction	\(1.05 \text{ }\mathrm{hr}^{-1}\)
\(k_G\)	Concentration of \(\mathit{glucose}\) required to achieve half of the maximum of cell growth rate	\(0.317 \%\)
\(k_{aTc}\)	Concentration of \(\mathit{aTc}\) required to achieve half of the maximum repression of cell growth rate	\(700 \text{ } \mathrm{nM}\)
\(\mathit{Yield}\)	Yield coefficient, biomass per mass of glucose utilized	\(4.48\)
\(m\)	Hill-coefficient for the repressive effect of Inducer(aTc) on cell growth rate	\(0.515\)
\(\mu\)	Rate of cell growth	-
\([\mathit{aTc}]\)	Concentration of the inducer (aTc) used to induce the transcription of SgrS mRNA	-
\(\mathit{OD}_{600}\)	OD₆₀₀	-
\([\mathit{glucose}]\)	Concentration of glucose in the growth media	-

SENSITIVITY ANALYSIS AND PREDICTIONS

The last and most exciting step of modelling was to use our models to gain useful insights into the systems, which helped shape the project’s design and set the direction for future work. After establishing the final models for the system, we carried out sensitivity analysis to identify parameters which could be tweaked during experiments to fulfill our experimental goals. As a part of the sensitivity analysis, the model parameters were varied by orders of magnitude and the change in the model output was studied.

The ability to run simulations to make predictions before carrying out the experiment is invaluable. It not only improves efficiency by reducing the need for trial and error but also provides insights that might otherwise go unnoticed. Thus, we have built our E.co Grow software with our models, to predict bacteria growth trends with user inputs like inducer concentration and other parameters.

Analysing Sensitivity

Results of our sensitivity analysis (Fig. 15) showed that manipulating parameters associated with the HicA toxin, as opposed to the HicB antitoxin, would result in a more significant difference in the growth trends observed. This implies that different extents of growth arrest and resumption could be more easily obtained by varying parameters pertaining to the production of toxin instead of the antitoxin. We ran simulations by varying a sensitive parameter like \(\mathit{deg_{HicA}}\) (Fig. 16) and found that a five-fold or a ten-told increase in this parameter can help in a precise control of the growth switch. These predictions proved to be valuable insights to the wet lab team as they could achieve these trends by doing something as simple as adding degradation tags to the HicA toxin.

Fig. 15: Sensitivity analysis of modifiable parameters of the model- \(\mathit{deg_{HicA}}\), \(\mathit{syn_{HicA}}\), \(\mathit{deg_{HicB}}\) and \(\mathit{syn_{HicB}}\). Sensitivity was analysed by studying the change in the growth profile from the optimized baseline model, as the different parameters were varied.

Fig. 16: Simulation of the effect of \(\mathit{deg_{HicA}}\) on cell growth. The black solid line represents the growth trend for uninduced sample; the red lines represent growth trends for induced HicA toxin with inducer (IPTG) concentration of 1mM. The solid, dashed and dotted lines are increasing values of \(\mathit{deg_{HicA}}\) , and the solid line is the optimized parameter value. A smaller \(\mathit{deg_{HicA}}\) results in greater growth arrest.

Optimising Experimental Design

For the Wet Lab to have a well-informed experimental design, we ran simulations and created a response surface that would allow us to identify the optimal combination of IPTG and Arabinose for inducing the toxin and antitoxin. The four considerations taken into account while choosing the best combination were:

a) Maximum arrest due to toxin

\[\Bigg(\frac{\mathit{OD_{ara}} + \mathit{OD_{\mathit{IPTG}+ \mathrm{2hr}}}}{\mathit{OD_{\mathit{IPTG} + \mathrm{2hr}}}}\Bigg)^{-1}\]

b) Maximum recovery due to antitoxin

\[\mathit{OD_{\mathrm{final}}} - \mathit{OD_{ara}}\]

c) Minimal dose of IPTG for toxin induction

\[\mathrm{Min}\Bigg(1, \frac{1}{\mathit{IPTG}}\Bigg)\]

d) Minimal dose of Arabinose for antitoxin induction

\[\mathrm{Min}\Bigg(1, \frac{0.2}{\mathit{ara}}\Bigg)\]

We wanted to ensure that the growth of the cells was truly arrested when there was an "OFF" mode induced by the Growth Switch. In order to account for this, we included an additional objective function (shown below) as a filter to penalize a growth of 10% and above during the arrested phase.

\[\frac{\mathit{OD_{ara}} - \mathit{OD_{\mathit{IPTG} + \mathrm{2hr}}}}{\mathit{OD_{\mathit{IPTG} + \mathrm{2hr}}}} < 0.1\]

These criteria were expressed as objective functions and the ideal combination was chosen to be one which maximises the value of these functions. The response surface generated (Fig. 17) allowed us to identify the area of concentrations we could use if we wanted to maintain the performance of the growth switch while minimizing the usage of inducers. Upon testing in the lab, it was found that the system still performed well for the optimum region of concentrations predicted by the model (Fig. 18).

Fig. 17: Response surface depicting the value of the objective functions for different concentrations of IPTG (mM) and Arabinose (%). Combination of IPTG and Arabinose concentrations in the red region as well as the surrounding blue region were tested by the wet lab team.

Fig. 18: Results of the characterisation experiment for the combinations belonging to the red and blue regions of the surface plot above. The trends observed in the experiment validated the predictions made by the model and the most optimal inducer combination was found to be 0.5mM IPTG (for toxin) and 0.4% Ara (for antitoxin).

Other Simulations

In addition to identifying sensitive parameters, we ran simulations by varying factors such as induction timepoints of toxin and antitoxin. Based on the results obtained, we understood that inducing the toxin or antitoxin at later time points do not allow for an observable growth arrest or growth resumption, respectively. As observed in the figure below (Figs. 19 and 20), timepoint of toxin induction is another parameter that can be controlled to achieve the desired extent of growth arrest.

Fig. 19: Simulations of bacterial growth with varying induction timepoints of toxin.

Fig. 20: Final OD₆₀₀ values for different induction timepoints of toxin. Inducing the HicA toxin at a later timepoint results in a less-pronounced growth arrest.

Fig. 21: Simulations of bacterial growth with varying induction timepoints of antitoxin. Inducing the HicB toxin at a later timepoint results in a less-pronounced growth resumption from the arrested mode.

Optimising Experimental Design

The most direct information the SgrS growth knob model could provide is the predicted cell growth and glucose concentration trends for a range of inducer concentrations (used to induce SgrS). From the simulations of growth with varying inducer concentrations, an inducer concentration for achieving a desired duration and rate of growth can be selected (Fig. 22 and Fig. 23). The amount of glucose remaining in the media at a timepoint can also be inferred from the simulated glucose trends in Fig. 24.

Fig. 22: Simulations of bacterial growth for varying inducer (aTc) concentrations used to induce SgrS. Increased inducer concentrations result in slower rate of growth, prolonging the duration of growth before reaching the final OD₆₀₀ value.

Fig. 23: Dose response of the maximum rate of cell growth with inducer (aTc) concentrations used to induce SgrS. Increased inducer concentrations result in the decrease of the maximum growth rate. The maximum growth rate drops steeply when inducer concentration is increased from 0nM to 400nM.

Fig. 24: Simulations of glucose remaining in media over time for varying inducer (aTc) concentrations used to induce SgrS. Increased inducer concentrations result in slower consumption of glucose, prolonging the duration of growth before all glucose is used up.

Analysing Sensitivity

In our sensitivity analysis, four parameters of our growth knob model, \(\mathit{deg_{aTc}}\), \(\mu_{\mathrm{max}}\), \(k_G\) and \(k_{aTc}\), that could be possibly be modified in our experimental design were selected. Growth was much more sensitive to \(\mu_{\mathrm{max}}\), \(k_G\) and \(k_{aTc}\) than to \(\mathit{deg_{aTc}}\) (Fig. 25). This meant that if we were to change to another inducer with a different degradation rate, the growth trends would not be significantly modified.

Fig. 25: Sensitivity analysis of the modifiable parameters of the model, \(\mathit{deg_{aTc}}\), \(\mu_{\mathrm{max}}\), \(k_G\) and \(k_{aTc}\). Sensitivity was analysed by studying the change in the growth profile from the optimized baseline model, as the different parameters were varied.

Other Simulations

To gain deeper insights into the effect of the three sensitive parameters on growth, we simulated growth trends of varying values of each parameter.

First we varied \(k_{\mathit{aTc}}\), the half-velocity constant for the inducer on cell growth. In our simulations, a smaller \(k_{\mathit{aTc}}\) value resulted in greater growth repression when SgrS is induced (Fig. 26). A smaller \(k_{\mathit{aTc}}\) corresponds a greater repressive effect of the inducer on growth, and can be achieved by using an inducer with a higher synthesis rate of SgrS.

Next we varied \(k_G\), the half-velocity constant for glucose on cell growth. In our simulations, a larger \(k_G\) value resulted in greater growth repression when SgrS is induced (Fig. 27). A larger \(k_G\) also means that the cell is less sensitive to glucose for growth, and can be achieved by using a strain that is less sensitive to glucose for growth.

Lastly, when we varied \(\mu_{\mathrm{max}}\), the maximum growth rate of the bacteria, we observed that a slower growth rate also resulted in a more significant reduction in growth when SgrS is induced (Fig. 28).

Therefore, our model predicts that using a stronger inducer to synthesize SgrS, choosing a strain of bacteria that is less sensitive to glucose for growth and has a slower growth rate, will result in a greater repressive effect of SgrS on cell growth.

Fig. 26: Simulation of the effect of \(k_{\mathit{aTc}}\) on cell growth. The black solid line represents the growth trend for uninduced SgrS; the red lines represent growth trends for induced SgrS with inducer (aTc) concentration of 1000nM. The solid, dashed and dotted lines are increasing values of \(k_{\mathit{aTc}}\), and the dashed line is the optimized parameter value. A smaller \(k_{\mathit{aTc}}\) value resulted in greater growth repression.

Fig. 27: Simulation of the effect of \(k_G\) on cell growth. The black solid lines represents the growth trend for uninduced SgrS; the red lines represent growth trends for induced SgrS at an inducer (aTc) concentration of 1000nM. The solid, dashed and dotted lines are increasing values of \(k_G\), and the dashed line is the optimized parameter value. Growth reduction due to induced SgrS (red) compared to uninduced SgrS (black) is more pronounced with a greater \(k_G\) value.

Fig. 28: Simulation of the effect of \(\mu_{\mathrm{max}}\) on cell growth. The black solid lines represents the growth trend for uninduced SgrS; the red lines represent growth trends for induced SgrS at an inducer (aTc) concentration of 1000nM. The solid, dashed and dotted lines are increasing values of \(\mu_{\mathrm{max}}\), and the dashed line is the optimized parameter value. Growth reduction due to induced SgrS (red) compared to uninduced SgrS (black) is more pronounced with a a slower growth rate (smaller \(\mu_{\mathrm{max}}\)).

LEARNINGS AND FUTURE WORK

Learnings

The insights that our models have provided to the experimental design have been summarized below.

Growth Switch Model	Growth Knob Model
The arresting of growth can be observed around 1hr after the induction of the HicA toxin. The resumption of growth from an arrested mode can be observed around 2hrs after the induction of the HicB antitoxin. The dose of toxin-antitoxin which give the best performance with minimum inducer usage was identified to be 0.5mM IPTG and 0.4% Ara.	The inducer (aTc) concentration required to slow growth is in the range of 10² and 10³nM. The slowing of growth due to the induction of SgrS can be observed in the first 3h of growth from 0.1 OD₆₀₀. Cell growth from 0.1 OD will reach final OD₆₀₀ at 16h, for an inducer (aTc) concentration of 1000nM.

Improvements

In order to make our predictions more well-informed, we decided to analyse the identifiability of the parameters in both our models. Our analysis involved examining the Posterior distributions (Kuzmanovska, 2012) of the different parameters. A larger distribution of the Posterior implies that the parameter is not converging and hence non-identifiable. We also created two-dimensional plots to identify potential correlation between parameters, which may contribute to their non-identifiability. The main aim of carrying out a deeper analysis of our models was to provide users with more representative predictions that would account for the uncertainty in the different parameters.

One-dimensional trace plots

According to the one-dimensional trace plots shown below, most of the parameters exhibit a Posterior that is distributed over a large range of values which pointed towards non-identifiability. These parameters are not converging to a single unique value but instead assuming a new fitted value for each run. However, some model parameters such as \(\mathit{deg_{HicA}}\) and \(\mathit{syn_{HicA}}\) display a relatively narrow range of Posterior which makes them identifiable to some extent.

Fig. 11: One-dimensional trace plots for parameters of the HicA-HicB growth switch model.

Two-dimensional trace plots

We also created two-dimensional trace plots for pairs of parameters that we suspected could be interdependent. Based on the trends in the trace plots, parameters such as \(\mathit{syn_{HicB}} and \mathit{k_{HicB}}\) appear to be highly correlated. The strongest correlation was observed between \(\mathit{OD_{\mathrm{min}}}\) and \(\mu_{\mathrm{max}}\) which is reasonable as one would expect the growth rate and inflection point on the growth curve to be related. As evident from the two-dimensional trace plots shown below, there are multiple sets of parameter values that are associated with a high Posterior. This implies that there is no one unique solution for the parameter values, essentially making them non-identifiable.

Fig. 12: Two-dimensional trace plots for parameters of the HicA-HicB growth control model.

Overall, the HicA-HicB growth switch model developed was of paramount importance to the success of our project as it provided useful insights pertaining to optimisation and predictability of the HicA-HicB growth switch . Upon deeper analysis, we found that the predictions woudl be more reliable if the identifiability aspect of the model could be improved. A way to improve the identifiability of our model is to have the wet lab team carry out experiments that allowed measurement of different variables, which would help introduce more constraints on the parameters. Another solution will be to develop an even simpler model that could still be representative of the biological system at but contain fewer unknown parameters.

One-dimensional trace plots

The one-dimensional trace plots in Fig. 13 showed a small Posterior distribution for all parameters of our growth knob system.

Fig. 13: One-dimensional trace plots for parameters of the SgrS growth control model.

Two-dimensional trace plots

In the two-dimensional trace plots below (Fig. 14), a narrow range for the sets of parameter values was observed. Although there was correlation observed between \(m\) and \(k_G\), their parameter values still converged within a small range.

Fig. 14: Two-dimensional trace plots for parameters of the SgrS growth control model.

Based on the convergence of parameter estimates and decent interdependence between parameters, we concluded that the SgrS growth knob model is fully identifiable. The results of this analysis strengthened the reliability of the predictions made using the SgrS growth knob model.

Future Work

Given that our technology aims to enhance the long-term productivity of bacteria, our next step will be to develop protein production models. Such models will allow us to predict the protein production trends of our growth switch and growth knob systems under specified conditions. We envision that users will one day be able to achieve desired bacteria protein production trends based on the predictions and recommendations of our robust models.

REFERENCES

Winter, A. J., Williams, C., Isupov, M. N., Crocker, H., Gromova, M., Marsh, P., … Crump, M. P. (2018). The molecular basis of protein toxin HicA– dependent binding of the protein antitoxin HicB to DNA. Journal of Biological Chemistry. https://doi.org/10.1074/jbc.RA118.005173

Rotem, E., Loinger, A., Ronin, I., Levin-Reisman, I., Gabay, C., Shoresh, N., … Balaban, N. Q. (2010). Regulation of phenotypic variability by a threshold-based mechanism underlies bacterial persistence. Proceedings of the National Academy of Sciences of the United States of America. https://doi.org/10.1073/pnas.1004333107

Fei, J., Singh, D., Zhang, Q., Park, S., Balasubramanian, D., Golding, I., … Ha, T. (2015). Determination of in vivo target search kinetics of regulatory noncoding RNA. Science. https://doi.org/10.1126/science.1258849

Kuzmanovska, I. (2012). Markov chain Monte Carlo methods in biological mechanistic models. ETH Zurich. https://doi.org/10.3929/ethz-a-007343293