### Dynamics at the single cell level

A synthetic lactose-luciferase (lac-lux) promoter controls the transcription of the lac repressor LacI (*X*), the activator LuxR (*Y*) and the enzyme LuxI (*Z*). Here, we consider the bacterium *Escherichia coli* as cellular chassis, and all proteins are assumed to be ssrA-tagged for enzymatic degradation [26], which for a zeroth-order kinetics enhances robustness of oscillations [27]. In addition, the chemical IPTG (*I*) acts as inducer and binds to LacI inhibiting its repressive effect; the chemical AHL (*A*) is required by LuxR to activate transcription. The active form of LacI is a tetramer (*X*_{4}), whereas for LuxR it consists of a dimmer of the complex LuxR-AHL ((*Y:A*)_{2}). In addition, repression by LacI tetramer induces a DNA loop in the promoter region (1a) [28], which may introduce a delay (*T*) into the system. The processes of transcription, translation, folding and multimerization could also induce a delay, but this is neglected in this work. Indeed, transcription in eukaryotes dictates a delay because of splicing [16]. However, this is not the case in prokaryotes. Herein, we assume that the reaction between two tetramers for making the DNA loop is not reversible, given that the looping structure remains even for low levels of LacI [22].

Furthermore, here we assume a fast mRNA dynamics (as the quasi-steady state is roughly reached in 5 min) [29]. The total concentrations of the three proteins (*X*, *Y* and *Z*) remain proportional, then we reduce the model to a single variable. The dynamics of the regulatory model is given by

where *C* is the plasmid copy number, *α* the nominal transcription/translation rate, *δ* and Λ the kinetic constants of the ClpXP protease, *μ* the cell growth rate, and *σ* the ratio activator/repressor (*X* = *σY* = *σZ*, being (1 + 2*σ*)*X* the total protein amount). LacI-mediated DNA loop enhances the autorepression. This is incorporated into the model by introducing an additional repressive term proportional to a looping constant (*L*_{0}). In case of no loop, *L*_{0} = 0. Moreover, this repression can be modulated by (LuxR:AHL)_{2}, as it occurs with CRP_{2}:cAMP in the natural lac promoter [23]. The functional form of the regulatory factor *f* has been previously studied [23, 28] and is given by

where *θ*_{
X
}is the binding coefficient of LacI_{4}-DNA, *θ*_{
I
}of LacI-IPTG_{2}, *θ*_{
Y
}of (LuxR:AHL)_{2}-DNA, and 1/*φ* accounts for the basal transcription rate. The parameter Ω accounts for the potential interaction of LacI_{4} and (LuxR:AHL)_{2} in the DNA loop. In case of no interaction, Ω = 1. In addition, the synthesis rate of AHL is assumed proportional to LuxI, and we consider that cells express the enzyme AiiA which degradates AHL [30]

where *κ* is the synthesis constant of AHL by LuxI, *δ*_{
A
}the thermodynamic degradation constant of AHL, and *λ*_{
A
}the degradation rate by AiiA.

The multimerization kinetics (at a given total cellular amount of AHL) is given by

where *k*_{l/-l}are the forward/reverse kinetic coefficients (*l* = 1, 2, 3, 4). By exploiting the different time scales in the dynamics and neglecting the amount of DNA-bound protein, we can define the dimensionless variables as *x* = *X*_{1}/*θ* with , *i* = *I*/*θ*_{
I
}, *y* = *Y*_{1}/(*K*_{4}*θ*_{
Y
})^{1/2}, and *a* = *A/K*_{3}. Notice that *y* and *a* depend on *x*. The reactions for multimerization can be assumed much faster than the ones for transcription and translation. Then the system (4) is reduced to the steady state. Being that, we obtain , *Y:A* = *Y*_{1}*A/K*_{3}, and , with *K*_{
l
}= *k*_{-l}/*k*_{
l
}.

Furthermore, the total amounts of LacI and LuxR are *X* = *X*_{1} + 2*X*_{2} + 4*X*_{4} and *Y* = *Y*_{1} + *Y:A* + 2(*Y:A*)_{2} respectively. For simplicity in the notation, we denote *X* as a function of the dimensionless variable *x*, , and . In addition, we assume that the total amount of AHL is approximately equal to the free one, then it turns out that . Being that, *Y:A* = 0 when *A* = 0 and at very high levels of *A*. For the following, we denote *x*_{
T
}= *x*(*t* -*T*). Time is also re-scaled by cell division, *τ* = *μ*_{0}*t/ln*(2) = *t*/*τ*_{0}.

Accordingly, the regulatory and degradation factors read, respectively, as

### Spectral analysis

We define and . Then, the differential equation (1) reads

for a constant value of *i*. We denote , and , where these variables are evaluated in the steady state (*x*_{
ss
}). The steady state satisfies . Being Δ*x* = *x* -*x*_{
ss
}, using first-order perturbations we can write

Thereby, the equation of eigenvalues (*λ*) is *λ* = Φ -Γ + Ψ*e*^{-λT}[31]. Let us define the following variables and *U*_{2} =(Γ - Φ)*T* to simplify the results of the spectral analysis. The oscillatory boundary condition is given by *λ* = *jω*, being *j* the imaginary unit. Then we obtain , the analytical equation for the Hopf bifurcation. On the contrary, the bistability boundary condition implies *λ* = 0, then we obtain *U*_{1} = -1. To further analyze the bistability of the system, we simplify the model at high levels of IPTG and externally introduced AHL. Without lost of generality, a dynamics governed by captures the principal features of the system, where *α'* ~*Cα*, *θ'* ~*σ*^{2}*K*_{4}*θ*_{
Y
}/*θ*^{2}, *δ'* ~*δ*, and Λ*'* ~Λ.

### Numerical integration and stochasticity

To illustrate this point, let us consider a molecular system governed by the following general differential equation , where *X* is the protein amount. The transcription term is highly nonlinear and accounts for the system delay (*T*). In our case, the function *f* depends on both *X*(*t*) and *X*(*t* -*T*). For the enzymatic degradation, we have to notice that, as Λ is a low value, this becomes zeroth-order for high values of *X*, while first-order when *X* is close to 0.

The delay-based system can be numerically integrated following

To solve this equation we use the MATLAB routine dde23.

To account for molecular noise we use the Langevin approach [25]. We introduce an intrinsic white noise into the model

where *η* (*t*) is a random fluctuation with ⟨*η* (*t*)⟩ = 0 and ⟨*η* (*t*)*η* (*t'*)⟩ = *δ* (*t*- *t'*) (Dirac delta). This equation was solved using the specified routine and considering a constant noise in the integration interval (Δ *t*)

where *ξ* is a Gaussian-distributed random number with mean 0 and standard deviation .

### Dynamics at the population level

LuxI catalyzes the production of AHL which can be pumped to the medium facilitating a cell-to-cell communication (i.e., QS) [32]. We also account for the dynamics of the intracellular and extracellular (labeled with subscript *e*) concentrations of IPTG and AHL, and for the cell population (*N*, with ) together with the spatial role in solid medium. The diffusion, the intracellular and the extracellular dynamics are given by

where *η* the equilibrium constant of the membrane transport of IPTG and AHL [33], *v* = *V/V*_{
e
}the ratio of volumes, *D* the diffusion (assumed linear and homogeneous) coefficient, and *N*_{
m
}the maximum cell capacity of the medium. Being the transport through the cell membrane fast, it turns out that *I* ≃ *I*_{
e
}and *A*_{
e
}≃ *QA* with

Since AHL is quickly degraded and *Q* ≃ 1 in a large population, we can take the quasi-steady state *A* = *κZ/λ*_{
A
}when AHL is not externally introduced. In addition, we neglect the movement of cells when replicating because even for *τ* = 100 this displacement would be ~0.1 mm.

Then, integrating the molecular and population models, the reaction-diffusion dynamics of the dimensionless cellular system is governed by

where the space is normalized by *D*. Solving the partial differential equation for IPTG diffusion in time and radial space (*r*) [34], the signaling pattern reads

where *r*_{0} is the radius of the initial drop of IPTG (here assumed 1 mm) and *I*_{0} is the modified Bessel function. Since for a given concentration of IPTG (*i*_{
c
}) the system is not oscillatory, we obtain, using the equation (14), the spatio-temporal limit

where, for small values of *r*_{0} (*r*_{0} << *Dt*), we have assumed that .