DNA looping increases the range of bistability in a stochastic model of the lac genetic switch

The two-state model contains two transcriptional states for the DNA, Off and On where the operon is transcriptionally inactive and active, respectively. DNA looping in not possible in this two-state model. When the operon is in the On state, transcription proceeds as a first order reaction

$$\begin{equation} \fl \mathrm{On} \mathop{\rightarrow}^{k_{\rm ts}} \mathrm{On} + m. \end{equation} \tag{ 1 }$$

Protein (Y) is translated at a rate that is proportional to the mRNA (m) copy number

$$\begin{equation} \fl m \mathop{\rightarrow}^{k_{\rm tl}} m + Y \end{equation} \tag{ 2 }$$

and protein and mRNA degrade at a rate proportional to their respective abundances

$$\begin{equation} \fl m \mathop{\rightarrow}^{k_{\rm degm}} \varnothing , \end{equation} \tag{ 3 }$$

$$\begin{equation} \fl Y \mathop{\rightarrow}^{k_{\rm degp}} \varnothing . \end{equation} \tag{ 4 }$$

Due to interactions between the lac repressor and inducer molecules, switching between the active and inactive transcriptional states occurs at a rate dependent on the inducer concentration inside the cell,

The switching rate functions k_fn([I]) and k_nf([I]) are determined by the microscopic interactions between the repressor and inducer. It was previously shown [33] that a particular microscopic model for these interactions gave rise to switching rates that were first-order for a given constant inducer concentration and with functional dependencies on inducer concentration that could be well-described by Hill-like functions. We fit the simulation data from [33] (see supporting information section S1 available at stacks.iop.org/PB/10/026002/mmedia) to the functions

$$\begin{equation} \fl k_{\rm fn}([{\rm I}]) = k_{\rm fn}^0 + \big(k_{\rm fn}^1 - k_{\rm fn}^0\big) \frac{[{\rm I}]^{H_{\rm fn}}}{ {I_{\rm fn}}^{H_{\rm fn}} + [{\rm I}]^{H_{\rm fn}}}\, \end{equation} \tag{ 6 }$$

and

$$\begin{equation} \fl k_{\rm nf}([{\rm I}]) = k_{\rm nf}^0 + \big(k_{\rm nf}^1 - k_{\rm nf}^0\big) \frac{[{\rm I}]^{H_{\rm nf}}}{ {I_{\rm nf}}^{H_{\rm nf}} + [{\rm I}]^{H_{\rm nf}}}. \end{equation} \tag{ 7 }$$

and obtain Hill-like expressions for the rate functions in our model in terms of internal inducer concentration. Here the parameters k⁰_x and k¹_x determine the minimum and maximum transition rates, I_x is the internal inducer concentration that the transition rate is $\frac{1}{2}(k_{\rm \it x}^0 + k_{\rm \it x}^1)$, and H_x is the Hill coefficient, which is proportional to the slope at [I] = I_x.

However, including the inducer as a separate species would dramatically increase the complexity of the model. Instead we express the transition rates in terms of the number of LacY proteins in the cell membrane,

Consider that inducer (I) can enter the cell either through passive diffusion from the environment (with constant concentration [I]_ex),

or through active transport by the LacY protein (modeled using irreversible Michaelis–Menten kinetics),

If we assume that inducer responds very quickly to a change in the LacY concentration, we can solve for the steady-state concentration of internal inducer as a function of the LacY and external inducer concentrations:

$$\begin{equation} \fl [{\rm I}]_{}([{\rm Y}]_{},[{\rm I}]_{\rm ex}) = [{\rm I}]_{\rm ex}\left( 1 + \frac{k_{\rm it}}{k_{\rm id}}\cdot \frac{[{\rm Y}]_{}}{[{\rm I}]_{\rm ex} + K_M} \right) \end{equation} \tag{ 11 }$$

with

$$\begin{equation} \fl K_M = \frac{k_{\rm yioff} + k_{\rm it}}{k_{\rm yion}}. \end{equation} \tag{ 12 }$$

Substituting equation (11) into equations (6) and (7) we obtain expressions for k_fn(Y) and k_nf(Y) with a little algebra.

Since the lac operon has three operators to which LacI can bind, the transcriptional states are more complicated than On and Off alone. The simplest modification to the two-state model to handle this complexity is to add a third state—Loop—that describes the operon when it is in a looped conformation with the repressor. It was hypothesized that in a looped conformation the LacI–O₁ binding can fluctuate allowing RNAP to bind and transcribe while O₁ is unoccupied with LacI nearby [3]. We introduce the Loop state to model this effect. This third state should not be interpreted as a specific conformation of the operon/repressor complex. Instead it should be understood as a coarse-grained state that represents all DNA/repressor looped conformations (O₁–O₂ and O₁–O₃) as well as bound conformations in which transcription is not repressed (O₂, O₃, or both bound). This state allows for the slow leakage of transcripts. We model this leakage by allowing transcription from the Loop state at a rate times the normal On transcription rate by adding the reactions

and

$$\begin{equation} \fl \mathrm{Loop} \mathop{\rightarrow}^{\epsilon k_{\rm ts}} \mathrm{Loop} + m \end{equation} \tag{ 14 }$$

to our model. The leakage factor is conditional probability to have repressor bound to O₂, O₃, or both O₂ and O₃ while in the Loop state. The schematic representation of the two models are shown in figure 1. This three-state model is similar to another three-state promoter model, which has been investigated recently [34]. This work focused on the stochastic mutual repressor model and described the state of the operators using three transcriptional states. Their methods would be difficult to implement for our models due to the noise in the translation rate arising from fluctuating mRNA abundance.

We estimated from [3]. They measured the distribution of YFP counts in a mutant with the lac circuit modified to disable positive feedback while allowing DNA loops to form. This experiment was conducted by replacing LacY with a Tsr-YFP fusion protein, which also localizes in the cell membrane, but cannot import allolactose into the cell. These distributions were well-described by a gamma distribution parameterized by the burst rate a and the burst size b. This approximate result was derived from a gene–mRNA–protein model with one transcriptional state [26]. The burst rate a ≡ k_ts/k_degp can be interpreted as the number of transcripts produced per cell cycle and the burst size b ≡ k_tl/k_degm refers to the average number of protein transcribed from a single mRNA. The burst rate and the burst size were measured for TMG concentrations between 0 and 200 μM and found to be relatively constant over this range: 0.34 < a < 1.25 and 1.29 < b < 2.59. We used this measurement to estimate that 5.7 × 10⁻⁴ < < 2.1 × 10⁻³ and assumed = 8.3 × 10⁻⁴ for our model. The two transcription rates k_ts and k_ts are responsible for the different burst sizes observed in [3]. 'Small' bursts of transcription occur due to leakage whereas 'large' bursts occur due to full dissociation of the complex (Loop→Off→On).

Since the On↔Off switching rate functions were well-described by a Hill function, we assume that the Loop to Off rate function could also take that form. Thus

$$\begin{equation} \fl k_{\rm lf}([{\rm I}]) = k_{\rm lf}^0 + \big(k_{\rm lf}^1 - k_{\rm lf}^0\big) \frac{[{\rm I}]^{H_{\rm lf}}}{ {I_{\rm lf}}^{H_{\rm lf}} + [{\rm I}]^{H_{\rm lf}}}. \end{equation} \tag{ 15 }$$

We will assume that the Off to Loop rate is constant in inducer concentration, since transitions into the looped state must occur via thermal fluctuations that take the singularly bound inducer/DNA complex into a conformation in which the unoccupied binding site on the lac repressor can reach a free operator. It is not possible to determine these switching rates directly from the data available from the literature, so we search the parameter space to determine for which values of these parameters the system exhibits the desired response. Of the parameters presented in table 1, only the five parameters describing the constant Off→Loop rate and the Loop→Off Hill function (15) are predicted. The remaining parameters are taken from the literature [3, 14, 33].

The deterministic rate equations for mRNA and LacY abundance are

$$\begin{equation} \fl \frac{{\rm d}x}{{\rm d}t} = k_{\rm ts}\mathcal {F}(y) - k_{\rm degm}x \end{equation} \tag{ 16 }$$

$$\begin{equation} \fl \frac{{\rm d}y}{{\rm d}t} = k_{\rm tl}x - k_{\rm degp}y \end{equation} \tag{ 17 }$$

where x and y are (continuous) concentration variables for mRNA and LacY respectively. Here we have defined the total transcription probability,

$$\begin{equation} \fl \mathcal {F}(y) = \frac{P_{\rm on}(y) + \epsilon P_{\rm loop}(y)}{ P_{\rm on}(y) + P_{\rm off}(y) + P_{\rm loop}(y) }, \end{equation} \tag{ 18 }$$

where the P_z(y) functions are the probability to find the system in the transcriptional state z at fixed protein abundance y and are computed simply from considering the steady state of a Markov process consisting of the three transcriptional states alone. This process is represented by

which at steady state is governed by the master equation

$$\begin{eqnarray} \fl \mathbf {A}\cdot \mathbf {P} = {\left(\begin{array}{@{}ccc@{}}-k_{\rm nf}(y) & k_{\rm fn}(y) & 0 \\ k_{\rm nf}(y) & -(k_{\rm fn}(y)+k_{\rm fl}) & k_{\rm lf}(y) \\ 0 & k_{\rm fl} & -k_{\rm lf}(y) \\ \end{array}\right)} {\left(\begin{array}{@{}c@{}}P_{\rm on} \\ P_{\rm off} \\ P_{\rm loop}\end{array}\right)} = 0, \nonumber\\ \end{eqnarray} \tag{ 20 }$$

whose normalized solution is

$$\begin{eqnarray} &&P_{\rm on} = \mathcal {C}^{-1} k_{\rm fn}(y)k_{\rm lf}(y) \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} &&P_{\rm off} = \mathcal {C}^{-1} k_{\rm nf}(y)k_{\rm lf}(y) \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} &&P_{\rm loop} = \mathcal {C}^{-1} k_{\rm fl}k_{\rm nf}(y) \end{eqnarray} \tag{ 23 }$$

with $\mathcal {C} = k_{\rm fn}(y) k_{\rm lf}(y) + k_{\rm nf}(y) \big ( k_{\rm fl} + k_{\rm lf}(y) \big )$, and y is computed from [I] by (11).

To investigate the fixed points of equations (16)–(17), we take the time derivatives to be zero and solve for y to get

$$\begin{equation} \fl \frac{{\rm d}y}{{\rm d}t} = k_{\rm degp}[\mathcal {N}\,\mathcal {F}(y) - y] = 0 \end{equation} \tag{ 24 }$$

where

$$\begin{equation} \fl \mathcal {N} = \frac{k_{\rm ts}k_{\rm tl}}{k_{\rm degm}k_{\rm degp}} \end{equation} \tag{ 25 }$$

is approximately the mean population of the induced state. For a set of parameters exhibiting bistability over some range of [I]_ex, the form of $\mathcal {F}(y)$ allows the dynamical system to exist in three distinct phases. For low values of [I]_ex, the system has a single stable fixed point at n_lo corresponding to the uninduced phenotype. For high values of [I]_ex the system has a single fixed point at n_hi corresponding to the induced phenotype. For intermediate values of [I]_ex, there is a range in which the dynamical system has three fixed points: two stable fixed points n_lo and n_hi and an unstable fixed point n₀. This external inducer regime corresponds to a heterogeneous population containing both phenotypes. The locations of these fixed points are not fixed, but can change with the inducer concentration. For the transcription state transition functions considered here, (24) is generally not solvable analytically due to the form of transcriptional state rate functions used.

The probability to be in a particular transcriptional state at a specific mRNA and protein abundance can be determined by solving the CME for the system [35]. These equations govern the time evolution of the probability, P^s_mn, to find m mRNA molecules and n proteins, when being in transcriptional state s. Writing out the transcriptional states explicitly we have

$$\begin{eqnarray} &&\fl \frac{{\rm d}P^{\rm on}_{mn}}{{\rm d}t} = k_{\rm fn}(n)P^{\rm off}_{mn} - k_{\rm nf}(n)P^{\rm on}_{mn} \nonumber\\ &&+\, \big [k_{\rm ts} \big(\mathsf {E}^{-1}_{M}-1\big) + k_{\rm degm}\big(\mathsf {E}^{+1}_{M} - 1\big)m \nonumber \\ &&+\, k_{\rm tl} \big(\mathsf {E}^{-1}_{N}-1\big)m - k_{\rm degp} \big(\mathsf {E}^{+1}_{N} -1\big)n\big ]P^{\rm on}_{mn} \end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray} &&\fl \frac{{\rm d}P^{\rm off}_{mn}}{{\rm d}t} =k_{\rm nf}(n)P^{\rm on}_{mn} + k_{\rm lf}(n) P^{\rm loop}_{mn} - (k_{\rm fn}(n)+k_{\rm fl})P^{\rm off}_{mn}\nonumber\\ &&+\, \big [ k_{\rm degm}\big(\mathsf {E}^{+1}_{M} - 1\big)m \nonumber \\ &&+\, k_{\rm tl} \big(\mathsf {E}^{-1}_{N}-1\big)m - k_{\rm degp}\big(\mathsf {E}^{+1}_{N} -1\big)n\big ]P^{\rm off}_{mn} \end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray} &&\fl \frac{{\rm d}P^{\rm loop}_{mn}}{{\rm d}t} = k_{\rm fl}P^{\rm off}_{mn} - k_{\rm lf}(n)P^{\rm loop}_{mn} \nonumber\\ &&+\, \big [\epsilon k_{\rm ts} \big(\mathsf {E}^{-1}_{M}-1\big) + k_{\rm degm}\big(\mathsf {E}^{+1}_{M} - 1\big)m \nonumber \\ &&+\, k_{\rm tl} \big(\mathsf {E}^{-1}_{N}-1\big)m - k_{\rm degp}\big(\mathsf {E}^{+1}_{N} -1\big)n\big ]P^{\rm loop}_{mn} \end{eqnarray} \tag{ 28 }$$

where $\mathsf {E}^{\pm i}_{M,N}$ are step operators that increment/decrement the mRNA index (M) or protein index (N) by i. The rates k_tl, k_degm, and k_ts,k_degp are the mRNA birth and death rates and protein birth and death rates respectively. The functions k_nf(n), k_fn(n), k_fl, k_lf(n), are the transcriptional state switching rates for On→Off, Off→On, Off→Loop, and Loop→Off.

In bacteria the ratio γ is typically large. We can use this fact to eliminate the mRNA degree of freedom from the CMEs. It is well-known that it is not sufficient to replace the mRNA dynamics with the average mRNA abundance while studying the switching times [24, 40, 41]. Transcriptional noise will lead to increased noise in protein abundance, which will affect the switching rates between induced and uninduced phenotypes. Another option would be to apply the quasi-steady-state approximation to the CMEs [42], however due to the stochastic switching of the transcriptional states, the transcription rate is not constant and thus the mRNA distribution is not at steady-state. We need to apply a different tack.

To remove the mRNA degree of freedom while still accounting for transcriptional noise we will exploit the fact that the number of protein molecules translated from a single mRNA is distributed geometrically. In the limit γ → ∞ an effective CME can be written without mRNA dependence by assuming that translation occurs in bursts with sizes that are geometrically distributed [43]. Figure 2(b) describes this pictorially. Since the probability that a single mRNA molecule is translated n times before decaying is

$$\begin{eqnarray} &&\fl P_{\rm tl.\,dec.}(n) \,{=}\, \left( \frac{k_{\rm tl}}{k_{\rm tl}{+}k_{\rm degm}} \right)^{\!n} \left(1- \frac{k_{\rm tl}}{k_{\rm tl}+k_{\rm degm}} \right) \,{=}\, P_{\rm tl}^n(1-P_{\rm tl}), \nonumber\\ \end{eqnarray} \tag{ 34 }$$

the translation terms in (26)–(28) can be rewritten as

$$\begin{equation} \fl k_{\rm tl}(\mathsf {E}^{-1}_{N} - 1)m \rightarrow k_{\rm ts}\left[ \left(1-P_{\rm tl}\right) \sum _{r=0}^n \left( P_{\rm tl} \mathsf {E}^{-1}_{N} \right)^{\!r} -1 \right], \end{equation} \tag{ 35 }$$

and we can drop the transcription terms, mRNA degradation terms, and mRNA indices. For example, consider the one-state model with the master equation

$$\begin{eqnarray} &&\fl \frac{{\rm d}P^{\rm }_{mn}}{{\rm d}t} = k_{\rm ts} \big(\mathsf {E}^{-1}_{M}-1\big)P^{\rm }_{mn} + k_{\rm degm}\big(\mathsf {E}^{+1}_{M} - 1\big)m P^{\rm }_{mn} \nonumber \\ &&+\, k_{\rm tl} \big(\mathsf {E}^{-1}_{N}-1\big)m P^{\rm }_{mn} - k_{\rm degp} \big(\mathsf {E}^{+1}_{N} -1\big)nP^{\rm }_{mn}. \end{eqnarray} \tag{ 36 }$$

This transforms to

$$\begin{eqnarray} &&\fl \frac{{\rm d}P_n}{{\rm d}t} = k_{\rm ts}\left[ \left(1-P_{\rm tl}\right) \sum _{r=0}^n \left( P_{\rm tl} \mathsf {E}^{-1}_{} \right)^{\!r} -1 \right] P_n \nonumber \\ &&-\, k_{\rm degp} (\mathsf {E}^{+1}_{} -1)nP_n. \end{eqnarray} \tag{ 37 }$$

All mRNA dependence is removed and now the translation term is that of a multi-step process where for a protein state n, transitions into that state can come from any protein abundance state from 0 to n − 1. Transitions out of the state n for translation can go to any state greater than n. These transitions are accounted for by the 1 in the translation term since the geometric distribution is normalized and the probability for a single transcription is P_tl. The protein decay term is not affected by this manipulation. This substitution is an adiabatic approximation that assumes that all mRNA is created or degraded between each protein reaction. This approximation is only valid when there is a large difference in time scales between the mRNA and protein dynamics. For E. coli, γ ≈ 53 based on degradation rates assumed for our model (cf table 1). This ratio is sufficient large that the mRNA and protein dynamics are well-separated. For the three-state model, this approximation reduces the dimensionality of the underlying system of equations by a factor of m_max (≈50) and significantly decreases the computational time. Analytically, this reduces (26)–(28) to a single dimension. A similar geometric bursting approximation for mRNA was used in a previous calculation [44] however only the average and variance of the protein abundance at steady-state could be calculated. To compute the full probability distribution of protein abundances and the switching rates between the induced and uninduced phenotypes a different method of solution must be employed.

Zoom In Zoom Out Reset image size

We numerically integrate the CMEs using the FSP to compute an approximate solution [32, 45]. This method is significantly faster than directly sampling the CMEs using the stochastic simulation algorithm (SSA). Computational time can be over 100 times shorter using this method compared to using the SSA with a sufficient number of realizations to be comparable to the error achieved from the FSP. To achieve this level of performance, this method truncates the state space at copy numbers that are not expected to be well-populated. The error due to the truncation is controlled adding fictitious states denoted sinks that accumulate probability that leaves the projection. Any transitions to states outside the projection are rerouted to these states. For example see figure 2(a), where the sink states are labeled η. The total probability lost due to the projection is the sum of probability in the sinks, which estimates the total error due to the truncation of state space. To find the optimal size of the projection, the system of ODEs is integrated while monitoring the accumulated probability in the sinks. If the sink probability increases past a threshold, the calculation is restarted with a larger projection space. The method for increasing the projection space depends on the type of problem being considered, for examples see [32, 46, 47]. For our simple 1D CME, we choose for our state projection $n\in [0,\alpha \mathcal {N}]$ where α is 1.25 initially and is increased when the calculation must be restarted.

Writing the CMEs in matrix form would lead to a large non-sparse matrix due to the tail of the burst distribution, because each state (n, s) can transition to (m, s) where $n<m\le \alpha \mathcal {N}$. Considering protein burst reactions alone, this situation leads to a block diagonal matrix where each of the three blocks are lower triangular. This problem can be solved by truncating the geometric series used in the burst approximation leading to a generator with a band structure. Any transitions from translation where the change in protein number is larger than a threshold are rerouted to a sink. Error due to the truncation of the geometric expansion is monitored by watching the accumulation in this sink. We found that 120 terms in the expansion were sufficient to achieve a total error of less than 10⁻⁷. Since (26)–(28) are not explicitly dependent on time, we integrated the equations using an off-the-shelf matrix exponentiation package [48].

We are interested in the mean switching rates between the induced and uninduced metastable states. These rates can be computed from the mean first passage time (MFPT) for the system to evolve from a stable fixed point of the deterministic rate equations (24)—either n_lo or n_hi—to the unstable fixed point n₀. Although there are elegant techniques available for estimating the stability of metastable states [49–51], we can compute these rates quickly and accurately using the geometric burst approximation coupled with the FSP. We compute the MFPTs by adding an additional sink at n₀, and integrating the equations with the initial condition that the copy number probability at the initial stable fixed point is unity [52].

The probability accumulated in this sink, η(t), is the cumulative distribution function (CDF) of first passage times. Aside from an initial transient period, the time evolution of η(t) is well-described by η(t) = 1 − e^−t/τ. We could then integrate the equations out to a time t_f where η(t_f) ≈ 1 and extract the MFPT τ by curve fitting. To save computational resources we use a different approach which does not require the equations to be integrated to t_f.

The mean of a distribution can be computed from its CDF by integrating the complementary CDF over its full domain. We will use this fact to estimate the MFPT by integrating over the known part of the full domain and adding to this value an estimate of the rest of the integral. At each time t_s during the integration of the CMEs, we integrate the complementary CDF 1 − η(t) numerically up to the current time step t_s and estimate the contribution from the integral from t_s to ∞ by assuming the complementary CDF decays exponentially for times greater than t_s. Thus at each time step we compute

$$\begin{equation} \fl \tau (t_s) = \int _0^{t_s}(1-\eta (t))\,\mathrm{d}t+\frac{1}{r(t_s)}\eta (t_s) \end{equation} \tag{ 38 }$$

where

$$\begin{equation} \fl r(t_s) = \left.-\frac{{\rm d}}{{\rm d}t}\ln \left(1-\eta (t)\right)\right|_{t=t_s}. \end{equation} \tag{ 39 }$$

We stop when the change in τ(t_s) per time step is less than a threshold:

$$\begin{equation} \fl \left| \frac{{\rm d}\tau (t)}{{\rm d}t} \right|_{t=t_s} < \eta _{\rm converge}. \end{equation} \tag{ 40 }$$

For a tolerance of total error less than 5× 10⁻⁴, the algorithm requires ≈300 steps to converge.

We define the macroscopically bistable range to be the range of inducer concentrations where the probability to be either lo or hi is greater than 10%. We find this range for a given set of parameters by first finding the range of external inducer concentrations that lead to three fixed points in the deterministic rate equations by counting the roots of (24). This range of concentrations is then searched, looking for where the probability to be in the induced state

$$\begin{equation} \fl P(\mathrm{hi}) = \frac{\tau _{\mathrm{hi}}}{ \tau _{\mathrm{hi}} + \tau _{\mathrm{lo}} } \end{equation} \tag{ 41 }$$

is 0.1 and 0.9. Here τ_x is the MFPT out of the phenotype x. We then search this bistable range for the concentration that either induction state is likely (P(hi) = 0.5), and ensure that the rate at the candidate concentration is greater the minimum acceptable value. The search is performed using the MFPTs computed from geometric burst/FSP method using standard root finding algorithms.

To test our results from the FSP calculations, we used the SSA [53] to compute the probability distribution functions (PDFs) and MFPTs. To look at the microscopic dynamics of switching events, the trajectories from the SSA simulations were segmented into induced, uninduced, and switching sections using thresholding heuristics that allow for arbitrarily long switching trajectories.

The heuristics first mark regions of the trajectory as lo if the protein number is lower than the unstable fixed point, and hi otherwise. An acceptable switching event is a segment of time t₀ long that does not leave lo (hi) followed by an intermediate segment less than δt_max long that fluctuates between the two states followed by a segment of time t₀ long that stays in hi (lo). The lead-in and lead-out times will not affect the results since the model is Markovian. We chose t₀ = 10 cell cycles. The maximum switching time δt_max was chosen so greater than 99.9% of the trajectories with long enough lead-in and lead-out times were accepted. The necessary δt_max depends on the particular parameter set but was generally of the order of 50 cell cycles.

A primary goal for the work is to determine if a minimal three-state model is able to recover the full range of bistability that was experimentally observed for the lac circuit with looping, compared to the behavior observed for the circuit where looping was prohibited. However, the model's looped state represents a coarse-grained state composed of many different microscopic states and no experimental data regarding microscopic transition rates are known. We therefore searched the space of all possible parameters related to the loop state to characterize how its addition changes the behavior of model. The Off→Loop transition has a single parameter k_fl and the Loop→Off transition function has four parameters: k⁰_lf, k¹_lf, I_lf, and H_lf. We randomly choose parameter values from this five-dimensional space while keeping all other model parameters fixed and tested the model's bistability properties.

To limit the size of the search, we bounded the parameter space such that only biologically reasonable parameter sets were sampled. The Off→Loop transition models the physical process of forming a DNA loop. If it is too slow to observe, the three-state model reverts the two-state model so we set the minimum rate of k_fl to be one event per 100 cell cycles. Likewise, loop formation requires the binding of a repressor–operator complex to a second operator and must be slower than repressor binding, so we set the maximum of k_fl to be one-half the rate of repressor binding, k⁰_nf/2.

For the Loop→Off transition Hill function k_lf(n), the parameter I_lf determines the approximate inducer concentration at which the transition rate will start to increase, i.e. at what inducer concentration the loop starts to become destabilized. We bounded this parameter such that the denominator in (15) is ${\mathcal {O}}(1)$ for inducer concentrations between 1 and 100 μM. This bound ensures that the transition rate function is able to switch between its low and high ranges over the experimentally relevant range of inducer concentrations. If I_lf was too large or too small, k_lf(n) would be effectively constant over the full range of LacY abundances and not provide feedback. The Hill coefficient H_lf determines the sharpness of the response and was restricted to be between 2 and 4 since inducer–repressor binding is cooperative with four binding sites. An in vitro study of DNA loop stability for the wild-type lac operon reported the mean lifetime of the loop to be on the order of the cell cycle [54]. In light of this observation, for the minimum transition rate k⁰_lf (in the absence of inducer) we generously bounded this value by setting the minimum rate for the search to be one dissociation per ten cell cycles and the maximum rate to be one-half of the repressor unbinding, k⁰_fn/2, since by definition loop dissociation must be slower than repressor unbinding. The maximum transition rate k¹_lf (in the presence of saturating inducer) was restricted to be faster than twice the zero inducer rate k⁰_lf but slower than 100 events per cell cycle.

We scanned 125 113 unique parameter sets by uniformly drawing a random value for each parameter from the above ranges. For each set we calculated the range of macroscopic bistability, defined at the low end by the inducer concentration where 10% of the steady state population is induced and at the high end by the inducer concentration where 90% of the population is induced (see section 3.3). Whereas the two-state model exhibits bistability near 10 μM (10.1–10.9 μM), the three-state model can produce bistability at inducer concentrations anywhere from 10.6–518 μM for sets of parameters constrained to the range above. Additionally, the range of macroscopic bistability for the two-state model is very small, 0.8 μM, while the addition of a third state can produce much wider ranges, up to 300 μM. As the overall switching rate to Loop decreases, the bistability range increases on average.

For each set we calculated the mean switching time at the inducer concentration where 50% of the steady-state population is uninduced and 50% is induced. At this concentration the time to switch to lo or to hi is identical and we use this time (τ_50%) as a measure of a parameter set's characteristic switching time. For the two-state model, τ_50% was 55 cell cycles, but with the three-state model τ_50% values could be obtained anywhere from 26–9.5× 10⁵ cell cycles. Decreasing I_fl has the effect on average of increasing the switching time τ_50% (see supporting information section S3 available at stacks.iop.org/PB/10/026002/mmedia).

Having performed a parameter search for the looped state, we then selected for further study only those parameter sets that reproduced the range of bistability seen in one experiment [3]. A parameter set was deemed acceptable if its range of bistability contained the range 50–60 μM and τ_50% was less than 100 cell cycles. Approximately 0.9% of the random samples satisfied these criteria. Figure 3 shows how the accepted parameters were distributed in the search space.

Zoom In Zoom Out Reset image size

The absolute value for several of the parameters appear to be unimportant to the switching properties. The Hill coefficient H_lf and the basal Loop→Off rate k⁰_lf are well-distributed within their allowed ranges. Likewise, the Off→Loop rate k_fl is uniformly distributed below 45 cell cycles. On the other hand, the inducer transition concentration I_lf and the maximum Loop→Off rate k¹_lf are more likely to take on values in certain ranges of the search space, high and low, respectively.

Some parameters showed a degree of mutual dependence. The basal Loop→Off rate k⁰_lf and the Off→Loop rate k_fl have a strong linear dependence (r² ≈ 0.6). Together with their uniform distribution, this dependence suggest that the total fraction of the time spent in the looped state in the absence of inducer, which determines the basal permease copy number, is an important property for switching, but the rate of cycling between the looped and unlooped states is relatively unimportant.

The inducer transition concentration I_lf and the maximum Loop→Off rate also appear to be dependent. All of the accepted sets fall above a convex region of the I_lf–k¹_lf phase space: k¹_lf ≲ I_lf⁴. If the maximum Loop→Off transition rate is high, then the transition concentration must also be high and if the transition concentration is low the Loop→Off transition rate must be low. A similar but less pronounced dependence is seen between the maximum Loop→Off rate and the Hill coefficient H_lf. This dependence is simply because there appears to be a maximum Loop→Off rate. When I_lf is large, k_lf(n) will not be able to reach its saturating level so values of k¹_lf larger than the maximum rate can be used.

The distribution of calculated switching properties for the accepted parameter sets is shown in figures 4(a)–(e). The onset of bistability, the external inducer concentration necessary for 10% induction probability can begin, as low as 30 μM and can last through 120 μM. The width of the bistable range can be from 20–80 μM. τ_50% values fall between 25–100 cell cycles. The properties of LacY in these cells is also of interest. The uninduced states were centered at zero LacY for all acceptable sets and the induced state means were found at 2350 ± 10.

Zoom In Zoom Out Reset image size

We tested the sensitivity of the model to changing the leakage factor . Over the likely range of from 5.7× 10⁻⁴ to 2.1× 10⁻³, the range of bistability only changed from 23.0 to 22.75 μM. The switching lifetime τ_50% only changed from 54.7 to 53.8 cell cycles.

We chose as a representative parameter set the mean of the acceptable parameter set distribution since that approximates the centroid of the hypervolume that encloses all acceptable parameters. In figure 3 this parameter set is indicated by a white point, the details of these parameters are enumerated in table 1, and the transcription state transition functions are plotted in figure 5. The PDF computed for this parameter set and the two-state model are shown in figure 6 compared to the PDF computed from the SSA. The combined FSP and geometric burst approximation predict the probability distributions with a high degree of agreement. The total probability lost from transitions out of the projection was less than 10⁻⁴, and the probability lost from truncating the geometric distribution was less than 10⁻⁷. The deviations of the FSP PDF from the SSA PDF are due to the adiabaticity assumption of the geometric burst approximation, not from probability loss from truncation. A numerical solution of the CMEs for the two-state model with explicit mRNA dependence showed exact agreement with the simulated distribution (see supporting information section S2 available at stacks.iop.org/PB/10/026002/mmedia) but with drastically increased computational cost.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The dependence of induction state lifetime on the external inducer concentration is generally exponential as shown in figure 7. Since changing the inducer concentration just affects the effective transcription rate, this functional dependence is in agreement with what is known about the dependence of switching rates in genetic networks [40]. This parameter set is bistable within the external inducer concentration range of 42–66 μM using 10% minority as the bistability metric. The switching time at equal induction probability is 63 cell cycles. For the two-state model, the bistability range was 10.1–11.0 μM with a switching time of 52 cell cycles. These values are in agreement with the results presented in [3] however the bistability of the non-looping mutant was only tested down to 20 μM external inducer so a transition from the induced phenotype to uninduced could not observed.

Zoom In Zoom Out Reset image size

The mean protein number in the uninduced states appear to be in agreement with experiment. This parameter set exhibits an uninduced state centered at 0 LacY whereas the two-state model has an uninduced state centered at 41. The location of the modes is in agreement with [3] since they showed for mutants unable to produce permease that the LacY distribution was centered at 0 for low inducer concentrations and the LacY distribution for mutants that cannot form DNA loops or produce permease was centered around 75. They did not report copy numbers for the induced state, however we found that adding DNA looping resolved the differences between the induced and uninduced states by moving the induced maximum from 1778 to 2351 and decreasing the width of the distribution by a factor of 2.7.

To investigate how the system transitions between the two phenotypes, we plotted protein and mRNA abundance histograms at regular time steps during a switching event in figure 8. We computed multiple trajectories using the SSA and searched them for segments in which the system switched between metastable states. All segments were then translated in time such that the first transcriptional state change event before switching is at t = 0.

Zoom In Zoom Out Reset image size

For the lo→hi transition, the three-state model shows quasi-deterministic switching behavior with an average switching time of 2.9 cell cycles. The duration and shapes of the protein trajectory during a switching event are all uniform. The two-state model does not exhibit this behavior. Instead the transitions are stochastic and irregular. The reason for this behavior appears to be that the three-state model stays in the Loop state in the uninduced phenotype and very rarely switches into Off. Prior to switching the mRNA population increases leading to more frequent bursts of protein, which can carry the system from the Off transcriptional state into On. Once On is reached switching back to Off is rare, leading to the deterministic switching behavior.

The hi→lo transition is also quasi-deterministic and is slower than the lo→hi transition (8.3 cell cycles). However the mechanism triggering the transition out of the induced state is different. Here random fluctuations of repressor binding are what drive the transition. The transcriptional state fluctuates into Loop by quickly passing through Off. If the system stays too long in the Loop state, the mRNA number will have decayed back to the distribution found in lo and the protein number will fall deterministically.