Integrative Circuit-Host Modeling of a Genetic Switch in Varying Environments

Synthetic biology is advancing into a new phase where real-world applications are emphasized. There is hence an urgent need for mathematical modeling that can quantitatively describe the behaviors of genetic devices in natural, fluctuating environments. We utilize an integrative circuit-host modeling framework to examine the dynamics of a genetic switch and its host cell in varying environments. For both steady-state and transient cases, we find increasing nutrient reduces the bistability region of the phase space and eventually drives the switch from bistability to monostability. In response, cellular growth and proteome partitioning experience the same transition. Antibiotic perturbations cause the similar circuit and host responses as nutrient variations. However, one difference is the trend of growth rate, which augments with nutrient but declines with antibiotic levels. The framework provides a mechanistic scheme to account for both the dynamic and static characteristics of the circuit-host system upon environmental perturbations, underscoring the intimacy of gene circuits and their hosts and elucidating the complexity of circuit behaviors arising from environmental variations.

involves a mechanistic description of host physiology including cellular growth, translation, and proteome partitioning which are intrinsically subject to environments. Therefore, the framework has the potential to capture environmental regulation of gene circuit behaviors (Fig. 1).
In this study, we explore the potential of the integrative modeling for capturing the behaviors of gene circuits in complex environments. Specifically, using the genetic toggle switch 7 as our model circuit and E. coli as our cellular host, we examine both the steady-state behaviors and transient dynamics of the circuit and the host upon nutrient and antibiotic variations. Our results elucidate complex environmental dependence of circuit behaviors and showcase the power of the integrative framework for understanding circuits under complex settings.

Steady-state switch and host behaviors at various nutrient levels. We began our exploration by
investigating the steady-state behaviors of the switch and its host under the alteration of nutrient availability. Specifically, we leveraged an integrative circuit-host modeling framework we recently developed 22 , which contains three fundamental parts including a coarse-grained host physiology description, a detailed kinetic module of synthetic circuits and layered circuit-host coupling. The host part utilizes the bacterium E. coli as an example and focuses on carbon fluxes from extracellular nutrient to cellular biomass using a coarse-grained approach [24][25][26] . The underlying processes include carbon uptake, ATP generation, amino acid synthesis, transcription and translation. The associated molecular species are carbon sources, ATP, amino acids, RNAs and proteins. Importantly, RNAs are binned based on their functions into three sectors, including mRNAs encoding proteins, tRNAs delivering amino acids to ribosomes, and ribosomal RNAs (rRNA) that are involved in ribosomes. Proteins and associated mRNAs are classified into three categories including ribosomal and affiliated proteins (R sector), metabolic proteins (E sector) and other proteins (Z sector). Additionally, the cell has an ability to adjust its proteome upon amino acid shortage, which is described by a kinetic description of the alarmone ppGpp 27,28 . The circuit part corresponds to biomolecules generated by synthetic gene circuits. In this study, the genetic toggle switch is the circuit of interest, which brings an additional sector, H, for both proteins and RNAs from gene circuits. The circuit-host coupling includes elementary and ppGpp-induced host-to-circuit interactions as well as load-and functionality-associated circuit-to-host interactions. In this study, we consider only the elementary host-to-circuit interaction and load-associated circuit-to-host interaction.
Using the integrative model, we first examined how nutrient alters the phase diagram of the circuit by simulating the circuit-host system and comparing its steady states under different nutrient levels ([n]) ( Supplementary  Information, Section 2). The values of the model parameters, except the induction strengths of Protein 1 and Protein 2 of the switch (k 1 and k 2 respectively), are all listed in Tables S2 and S3, with the former containing parameters adopted from the previous model 22 and the latter containing those introduced in this study. For each nutrient level, the phase diagram was numerically drawn in the k 1 − k 2 space by comparing the circuit's steady states from the simulations starting with two initial conditions. Initial condition 1 (I.C. 1) corresponds to a high level of Protein 1 ([P 1 ]), a high level of mRNA 1, a low level of P 2 ([P 2 ]), and a low level of mRNA 2 while initial condition 2 (I.C. 2) is the opposite (Supplementary Information, Section 2). The system is defined bistable if the steady states from the two initial conditions do not converge and monostable otherwise. Figure 2a shows the overlaid phase diagrams for the nutrient levels of 10 (light blue), 33 (blue) and 1000 (dark blue) μM. We found that increasing the nutrient level shifts the bistability diagram to the upper right, which is consistent with previous reports 29 . The reason is following: nutrient augment increases the host's growth rate, which enhances the dilution of the proteins; as a result, the circuit production needs to be induced at a higher level to overcome the growth-induced protein dilution for the switch to maintain stable steady states. For the same reason, for a given induction of k 1 and k 2 (e.g., the dot in Fig. 2a at k 1 = 25 hr −1 and k 2 = 28 hr −1 ), the circuit, which is originally bistable at a low nutrient, can become monostable at a high nutrient.
In addition to phase diagrams, we investigated the responses of specific key circuit and host variables upon nutrient variations. Figure 2b,c shows the level of circuit proteins (Proteins 1 and 2) as a function of nutrient for a given set of production induction levels (k 1 = 25 hr −1 and k 2 = 28 hr −1 ). Notably, the levels of both proteins have two stable states at a low level of nutrient, while beyond the critical level (45.4 μM), the circuit becomes monostable. Correspondingly, host variables also show the same bistability-to-monostability transition. Figure 2d,e illustrates the growth rate of the host upon nutrient alteration and the corresponding partition of cellular proteome,  Fig. 2d). This observation is owed to the larger metabolic load of Protein 2 than that of Protein 1. Associated with the growth rate, proteome allocation (Fig. 2e) exhibits a similar correspondence where, in the steady state of low [P 1 ] and high [P 2 ], the host has a greater heterologous H sector (blue lines) but reduced R (orange lines), E (red lines), and Z (green lines) sectors.
To intuitively illustrate how nutrient level alters circuit behaviors, we simulated the temporal evolution of key circuit and host variables for different initial conditions and varied nutrient levels. To implement our simulations, the initial [P 1 ] was varied from 0 to 2000 nM while the initial values of other variables were kept invariant. Figure 3 shows the transient behaviors of the H, R, E, and Z sectors (top to bottom rows) of cellular proteome at the nutrient levels of 10 (left column), 33 (middle column) and 1000 μM (right column). Accordingly, Fig. S1 shows the behaviors of the corresponding [P 1 ], [P 2 ] and growth rate at each of these nutrient levels. Notably, the circuit and host variables diverge to two distinct stable states for different initial conditions at the nutrient levels of 10 μM and 33 μM, but converge into single stable states at the nutrient of 1000 μM. Consistent with our findings in Fig. 2, these results demonstrate dramatic nutrient modulations of the circuit-host system. Importantly, the above results reveal not only the shift of the circuit's phase diagram but also the steady-state characteristics of the host and its modulation by nutrient, which elucidates the intricate dependence between the circuit and the host and showcases the power of integrative modeling.

transient responses of the switch-host system upon nutrient shifts. Our steady-state results
showed that both the circuit stability and host physiology can be influenced by nutrient conditions. This inspired us to examine the transient properties of the circuit-host system upon nutrient shifts because, in nature, nutrient level varies constantly. A unique feature of our integrative model is the ability to quantitatively capture the dynamic behaviors of both circuit and host 22 . Thus, we utilized the framework to determine how the circuit and the host respond transiently when the environment is switched between rich and poor nutrient conditions. We first performed an in silico nutrient upshift. Specifically, at t < 0, the system was settled in either of the two steady states (solid and dashed lines) at a low nutrient level (10 μM); then, at t = 0, the nutrient was switched to a high level (1000 μM) (Fig. 4a). The experiment showed that the upshift causes [P 1 ] of the circuit, which is originally either in the stable low (dashed lines) or high (solid lines) states, to adapt to a single new steady state after a brief transient (Fig. 4b). This transition is due to the fact that, at a nutrient level of 10 μM, the system is bistable and both the circuit and host variables possesses two stable states but at 1000 μM the system becomes monostable. Accordingly, the growth rate, initially in one of the two stable states, converges to a single rate ( Fig. 4c) with the increase of nutrient level. The cellular proteome also transits from two stable partitions to a single, new partition as shown in Fig. 4d. www.nature.com/scientificreports www.nature.com/scientificreports/ We then tested the transient responses of the circuit-host system upon a nutrient downshift from 1000 to 10 μM, which corresponds to the alteration of the system from monostability to bistability. In principle, in a bistability regime a system in different initial states may evolve to different steady states. We thus speculated that, depending on the downshift time, the system can exhibit distinct transient patterns as well as steady-state behaviors. To test this speculation, we conducted in silico nutrient downshifts at different time points after the initiation of temporal evolution (Fig. 4e-h). Indeed, when the nutrient downshift occurs early (e.g., the transition time T s < 0.1 hr), [P 1 ] evolves to two distinct stable states for different initial conditions (Fig. 4e). With the increase of the shifting time (5 hr and 10 hr), fewer trajectories lead to the high [P 1 ] state ( Fig. 4f-g). Eventually, when the shift occurs after a sufficiently long time (T s > 20 hr), the system converges into the same steady state regardless of its initial states (Fig. 4h). The results demonstrated that a fully dynamic, mechanistic model is capable of elucidating transient circuit and host dynamics for nutrient shifts. The results also showed that the shifting time can be critical to the final state of the switch-host system.

Steady-state switch and host behaviors under varied antibiotic stresses.
Like nutrient availability, environmental stresses often vary in natural settings and have prominent impacts on cell physiology. The integrative model incorporates key cellular processes such as translation, which allows us to mechanistically describe cellular stresses such as the effects of chloramphenicol, a ribosome-inactivating antibiotic 30 . We thus used the framework to explore how chloramphenicol concentration ([C m ]) affects the steady-state switch and host behaviors in a rich nutrient (1000 μM) environment.  Figure 5a shows the bistability phase space of the circuit for varied chloramphenicol levels (0 (dark blue), 4 (dark green) and 12 (green) μM). Here, the bistability region for [C m ] = 0 μM is identical to the case of the nutrient level 1000 μM in Fig. 2a. These results show that the bistability region reduces with increasing [C m ], which can be understood as follows: the percentage of deactivated ribosomes increases with enhancing [C m ], which effectively reduces the overall translational capacity of the cell. To remain bistable, the toggle switch requires higher induction levels (k 1 and k 2 ) to overcome the reduction in translational capacity. Consequently, the circuit transitions from bistability to monostability with increasing [C m ] for a given set of k 1 and k 2 such as the dot in www.nature.com/scientificreports www.nature.com/scientificreports/ Fig. 5a (k 1 = 29 hr −1 and k 2 = 31.6 hr −1 ). Additionally, the results show that the lower edge of the phase boundary is affected differently from the upper edge, which is due to the differential loads of the two circuit proteins.
Along with the bistability analysis, we also examined the responses of circuit and host variables to chloramphenicol variations. Figure 5b,c shows steady-state circuit protein behaviors for a given set of induction strengths, corresponding to the dot in Fig. 5a, when [C m ] is systematically altered. Clearly, the circuit proteins change from being bistable to monostable when [C m ] is beyond a critical level (4.2 μM). Interestingly, the maximal [P 1 ] and the minimal [P 2 ] from the initial condition 1 (solid lines in Fig. 5b,c) occur at a non-zero level of chloramphenicol ([C m ] = 1.2 μM), which arises from the balance between the reduction in translational capacity and reduced nutrient uptake. In concert with circuit protein profiles, the steady-state growth rate displays a transition from bistability to monostability (Fig. 5d) as is the case for nutrient shift (Fig. 2). However, the difference is that the overall growth rate increases with nutrient but declines monotonically with chloramphenicol due to the differential roles of the two environmental parameters. Similarly to the growth rate, the steady-state proteome allocations among the four coarse-grained sectors are modulated (Fig. 5e).
We further simulated the temporal dynamics of circuit and host variables at different initial settings where the initial [P 1 ] is altered from 0 to 2000 nM while others remain constant ( Fig. 6 and Fig. S2). Consistent with our findings in Fig. 5, both the circuit (H sector) and the host parts (R, E, and Z sectors) of the proteome show bistability at low [C m ] (the left and middle columns) but become monostable at high [C m ] (the right column).

Dynamic responses of the switch and the host to antibiotic shifts.
We also investigated the transient dynamics of the switch-host system upon antibiotic alterations. Similar to the case of nutrient shifts, the circuit-host system was first allowed to settle into one of the two steady states in the absence of chloramphenicol (i.e., [C m ] = 0 μM). Then, chloramphenicol level was shifted from 0 to 12 μM (Fig. 7a) at time 0 and, subsequently, transient circuit behaviors and the host physiology were examined (Fig. 7b-d). Because the system is bistable at [C m ] = 0 μM but monostable at [C m ] = 12 μM, the upshift of the antibiotic leads to the convergence of [P 1 ]. The system, initially in either of the two steady states, develops into a single steady state after a long transient time (Fig. 7b). Accordingly, the host variables, including both the growth rate (Fig. 7c) and proteome partition (Fig. 7d), experience the similar stability change. We also tested how chloramphenicol relaxation modulates the switch and host dynamics. Figure 7e-h shows the transient behaviors of the circuit's Protein 1 when chloramphenicol is reduced from 12 to 0 μM at different time points. Because the system is monostable at [C m ] = 12 μM but bistable at [C m ] = 0 μM, Protein 1 can evolve to different steady states, depending on differential initial conditions at time 0, when the downshift occurs relatively early (Fig. 7e-g). However, when the transition occurs late, Protein 1 converges to a single steady state regardless of its initial condition (Fig. 7h).

Discussion
After successful creation of a wide array of engineered gene circuits, synthetic biology is now advancing into a new era wherein circuits are deployed in the real world. To aid in this transition, mathematical models are needed to quantitatively understand genetic devices in natural environments which fluctuate constantly. Our study shows that an integrative approach has the potential to mathematically describe environmental dependence of circuit-host dynamics. Specifically, the framework provides a mechanistic scheme to account for both the steady-state characteristics and transient responses of the circuit-host system upon nutrient and antibiotic variations. Our results also illustrate the complexity of gene circuit behavior in changing environments through intimate circuit-host coordination.
There are few remarks for our study. First, the proteomic fraction dedicated to heterologous proteins can be high. One key feature of our framework is being coarse-grained, which means proteins with similar biological functions or similar dynamic patterns can be grouped into a single effective protein as a way to simplify representation. Following this spirit, the proteins P 1 and P 2 here represent not only the specific transcriptional factors in the toggle switch (i.e., CI and LacI respectively) 7 but also heterologous proteins whose productions are controlled by the switch. As the switch can be used to drive a number of genes with desired functions such as complex biosynthetic pathways 31 , the proteins encoded by the switch genes and the genes driven by the switch can be coarse-grained into effective P 1 and P 2 . As a result, heterologous proteins can take up a large sector of the total proteome. To support the argument, a recent experiment showed that proteins over-expressed from a single gene (lacZ) can constitute about 30% of the total proteome 25 . Second, our results seem different from conclusions from a different recent coarse-grained approach 32 . Extended from the two-variable ODE representation in the original switch model 7 , the coarse-grained model in that study includes circuit-host competition by revising the production terms of the two proteins. However, it does not involve a detailed mechanistic description of host variables or resource allocation dynamics like our framework does. Additionally, the form of the revised production terms suggests that the other model primarily considers the competition of machineries such as RNAP and ribosomes; in contrast, our model includes the competition of machineries, energy and building blocks. Because of fundamental differences of model assumptions, the corresponding results may not be directly translatable in a www.nature.com/scientificreports www.nature.com/scientificreports/ strict one-to-one way with ours. One more key difference is that, in our study, we explore circuit dynamics under varied nutrient and antibiotic levels instead of metabolic load, because the focus of our paper is to investigate how environmental factors shape circuit dynamics. Third, while our study focuses on the behaviors of the toggle switch in E. coli only, one may be curious whether the results apply to different circuits and cellular hosts. For different synthetic circuits in the same host, we speculate that system behaviors can be described by revising the detailed circuit kinetics and recalibrating the circuit-host coupling while keeping intact for the host physiology. By contrast, for the same circuit in different organisms, the circuit module can remain unchanged but the description of the host physiology and circuit-host coupling need to be updated. Generally, we speculate that the degree of model revision correlates positively with the phylogenetic relevance between E. coli and the host of interest. As supported by our previous work 22 , our original model was successfully adapted to describe key physiological variables of Salmonella typhimurium, a species of proteobacteria closely related to E. coli, without any alterations to the model parameters. For comparison, to capture the more distantly-related, Gram-positive Streptomyces coeliocolour, the value of one key parameter needed to be revised while keeping the model structure and other parameters invariant.
As an initial exploration, this work demonstrates the promise of integrative modeling for understanding circuit behaviors in simple varying environments. Moving forward, it is valuable to examine circuit dynamics in more complex settings, such as environments with periodic or random external fluctuations 33 , intracellular conditions where biomolecules fluctuate intrinsically 34 and habitats where multiple species of microbes coexist and their interactions vary 12,[35][36][37] . More importantly, experiments need to be performed to test model predictions and, through experiment-modeling iterations, to achieve a quantitative understanding of circuit dynamics in complex environments.

Methods
construction of a circuit-host model for the genetic toggle switch. To model circuit dynamics in varying environments, we employed an integrative circuit-host modeling framework we recently developed 22 . Different from traditional models focusing on circuit alone, the framework involves three fundamental parts: a coarse-grained description of host physiology, a detailed module for exogenous circuit kinetics and multiple layers of circuit-host coupling. The framework possesses a mechanistic description of host physiology, which is intrinsically subject to environmental variations, and an explicit consideration of circuit-host coupling which organically links circuit dynamics to host behaviors. Thus, the model is uniquely suited to model gene circuits in naturally fluctuating environments. For this study, the host part can be captured using a 12-variable system of ordinary differential equations (ODEs) focusing on amino acid flux. The circuit part can be modeled using four ODEs corresponding to the kinetics of the mRNAs and proteins of the two genes of the switch. Notably, one of the circuit proteins (Protein 2, P 2 ) is assumed to consume more resources than the other (Protein 1, P 1 ), as the case of experimental study 7 . Additionally, for simplicity, the two proteins are assumed functionally neutral. We primarily considered variations of nutrient and chloramphenicol levels in environments, because they are two key parameters regulating environments and essential to cellular growth and metabolism. A mathematical model of 16 variables was constructed to study the toggle switch and its host, which is detailed in Supplementary Information (Section 1). The ODEs were numerically simulated in MATLAB (MathWorks) using the ode15s solver with a relative and an absolute tolerance of 10 −6 . parameters and initial conditions. Detailed information relating to the parameters, simulation methods and initial conditions of the study is provided in Supplementary Information (Section 2 and 3).