Spaceinduced bifurcation in repressionbased transcriptional circuits
 Amanda Lo Van^{1, 2},
 Hedi A Soula^{1, 3} and
 Hugues Berry^{1, 2}Email author
https://doi.org/10.1186/s129180140125z
© Lo Van et al.; licensee BioMed Central Ltd. 2014
Received: 20 July 2014
Accepted: 25 September 2014
Published: 8 November 2014
Abstract
Background
Albeit the molecular mechanisms of gene expression are well documented, our understanding of their dynamics is much less advanced. Recent experimental evidence has revealed that gene expression might be accurately organized in space, with several molecular actors localized to specific positions in the cell. However, the influence of this spatial localization on the dynamics of gene expression is unclear. This issue is also central in synthetic biology, where one usually considers the spatial localization in the cell of the genes of the inserted synthetic construct as irrelevant for its temporal dynamics.
Results
Here, we assessed the influence of the spatial distribution of the genes on the dynamics of 3gene transcriptional ring networks regulated by repression, i.e. repressilator circuits, using individualbased modelling to simulate their dynamics in two and three space dimensions. Our simulations suggest that variations of spatial parameters  namely the degree of demixing of the positions of the gene or the spatial range of the mRNA and proteins (i.e. the typical distance they travel before degradation)  have dramatic effects by switching the dynamical regime from spontaneous oscillations to a stationary state where each species fluctuates around a constant value. By analogy with the bifurcations arising from the variation of kinetic parameters, we referred to those transitions as spaceinduced bifurcations.
Conclusions
Taken together, our results strongly support the idea that the spatial organization of the molecular actors of transcriptional networks is crucial for the dynamics of gene expression and suggest that the spatial localization of the synthetic genes in the cell could be used as an additional toggle to control the dynamics of the inserted construct in synthetic biology experiments.
Keywords
Background
In a seminal paper published more than 50 years ago [1], F. Jacob and J. Monod proposed a generic mechanism for protein synthesis in cells, whereby a DNA gene produces a messenger RNA molecule (mRNA) which is then used to produce the corresponding protein. They also described how this gene expression process is controlled by cytosolic macromolecules called repressors, that stop the expression of a given gene by binding to it. Since then, it was discovered that repressors are actually proteins (or sometimes RNAs) produced by other genes and that positive effectors (activators) also exist [2]. Fifty years of molecular studies have unravelled the complexity of the underlying molecular machinery [3],[4], but the existence and biological significance of such repressionbased transcription systems have been thoroughly confirmed.
The design, insertion in cells and theoretical analysis of small synthetic transcriptional regulation networks have proven to be very useful to our understanding of the temporal dynamics of gene networks. From a purely theoretical standpoint [5][7], repressionbased transcriptional regulation loops (i.e. ring networks), are generically expected to exhibit bistability with hysteresis (even numbers of genes) or a regime of spontaneous periodic oscillations (odd numbers of genes). Accordingly, insertion of 2gene repressionbased synthetic loops in bacteria indeed may give rise to bistable dynamics [8], whereas 3gene repressionbased synthetic loops inserted into living bacteria are capable of spontaneous oscillations with a very long period (more than 2.5 hours, i.e. twice the cell doubling time) [9]. In both articles [8],[9], the construction of the synthetic networks and their insertion and study in living bacteria were accompanied by a short theoretical study explaining why the observed dynamics were to be expected. In both cases, the complex dynamics of interest arose in the models via a bifurcation from a unique stable steady state: varying one kinetic parameter, the model predicts the occurrence of a saddlenode bifurcation supporting bistability (2 genes) or a Hopf bifurcation explaining the appearance of a limit cycle and its spontaneous oscillations (3 genes). Interestingly, both models consisted of ordinary differential equations (ODE) that assume mass action kinetics [10] for all the reactions. This corresponds to a strong assumption about the internal medium of the cell that is supposed to be dilute, perfectlystirred and spatially homogeneous. Beyond its importance for switchlike or oscillator circuits, this view in fact has deep impact in the whole field of synthetic biology. When inserting synthetic gene network constructs in chassis cells, synthetic biologists usually do not consider that the spatial localisation of the corresponding plasmids or of the insertion points on the chromosome is an important parameter for the temporal dynamics of the synthetic construct. Whereas this viewpoint seems reasonable if the hypothesis of a perfectlymixed internal space is valid, it should be questioned if spatial homogeneity is violated.
Actually, the traditional view of the interior of the cell as a perfectlystirred and spatially homogeneous medium, where the concentration of each reactant would be identical wherever one measures it inside the cell, has increasingly been challenged by recent results, especially in bacterial cells. First, recent advances in the measurement of single particle trajectories inside living cells have unravelled that the bacterial cytoplasm is a complex, extremely crowded and dense medium that strongly affects molecule mobility in a spatially nonhomogeneous way [11][15]. Therefore, due to the intrinsic physical nature of the cytoplasm, mobility inside the cytoplasm may not guarantee perfect mixing and homogeneous spatial distribution of its constituents. Recent experimental results have shown that the position of chromosomes in the nucleus of eukaryotic cells or the location of the genes on the bacterial chromosome are not random or perfectlystirred but selforganized to sit on specific locations forming spatial maps that are stable over long time scales [16],[17]. Even in bacteria, transcription is believed to be organized spatially, with the major molecular actors assuming stable and defined intracellular locations [18][20].
Despite this accumulating evidence of the localization of the elements of transcriptional regulation networks inside the cell, the influence of spatial properties on the temporal dynamics of gene expression remains to be fully described. Indeed, in the simplest instance of classical biochemical reactions, the impact of spatial localization has only recently received attention in the case of enzyme complexes [21][23] and in the case of signalling clusters in membrane domains [24][26]. All these results show that spatial correlation strongly modifies the apparent chemical affinity involved in the pathway both for the transient and equilibrium behaviors. Additionally, some evidence suggest that this depends strongly on the actual diffusion values.
In this paper, we use computer simulations to explore if and how the localization of the genes can influence the dynamics of small repressionbased transcriptional regulation ring networks. We focus on repressionbased transcriptional regulation loops composed of three genes, i.e. repressilators. Using stochastic spatiallyexplicit individualbased computer simulations, we find that the localization in space of the genes is of crucial importance to the dynamics of the system since it controls even the global dynamics regime, i.e. whether the system fluctuates around stationary values or exhibits spontaneous oscillations. We show that when parameters related to the spatial organization of the genes change, the repressilator undergoes a sharp transition between the oscillatory and stationary regime. Effective control spatial parameters include the degree of demixing of the gene locations or the spatial range (i.e. the typical distance travelled before degradation) of the mRNA and proteins. Since this transition is very similar to the bifurcations along the kinetics parameters that are usually evoked to explain the appearance of the oscillatory regime, we refer to it as a spaceinduced bifurcation. Our results therefore unravel the importance of spatial properties in the dynamics of transcriptional regulation networks. Moreover, they suggest that the spatial localisation of the synthetic genes in the cell could be used as an additional toggle to control the dynamics of the inserted construct in synthetic biology experiments.
Methods
Our objective is to study the dynamics in time and intracellular space of a generic transcriptional circuit in a bacterial cell like E. coli. Extensively detailed realistic modelling of bacteria (wholecell models), with experimentallyderived values for most rate constants, realistic cell space volume, protein size and diffusion coefficients, and interaction with the metabolism, starts to be accessible to today computer power [27][30]. However, those models bear the limitation that many of the parameters still lack experimentallyderived values. Moreover, their computational cost still forbids to use them for thorough exploration of the parameter space with reasonable statistical significance, which is precisely the objective of the present work. We therefore opted for simpler models and restricted our focus to the study of the major processes implicated in transcriptional regulation (and described in the seminal Jacob and Monod paper [1]) rather than taking into account extensive details of the cell intracellular machinery.
Gene expression model
where ; α and β are the transcription and translation rates, respectively; τ _{ M } and τ _{ P } are the lifetimes of mRNA and proteins; and denote the singly and doublybound genes, respectively and o kon, o koff, o koff2 the corresponding reaction rate constants. Note that to keep the model as simple as possible, the values of those constants were chosen identical for all genes i.
Mass action kinetics approximation
where we used mass conservation for the different bound fractions of the genes, i.e. .
Parameter values and numerics
Numerical integration (4th order RungeKutta method) and bifurcation analysis were performed with Xppaut (available at http://www.math.pitt.edu/~bard/xpp/xpp.html). To keep the model as generic as possible, we expressed time in multiples of the integration time step (△ t=1), thus yielding reaction rate constants expressed as inverse of time steps. The bifurcation diagram was explored along the protein lifetime T _{ P } and the transcription rate α (Figure 1BC). Those parameters were thus varied over several orders of magnitude. The values of the translation rate β and the mRNA lifetime T _{ M } were set so as to ensure realistic copy numbers of proteins and mRNAs. mRNA copy number in single E. coli cells was quantified between <0.1 and 5 [31] or even up to 50 for strongly expressed promoters like the P_{ lac } promoter under strong induction [32]. Protein copy numbers vary over a wide range, from 10 to 10^{5}[31],[33], resulting in a protein/mRNA copy number ratio that is roughly between 10^{2} to 10^{4}[31]. In our model, the range of mRNA/Protein copy number M _{ i }/P _{ i } is given by the product of parameters t _{ P } β. Since T _{ P } was varied between 10^{3} and 10^{6} in the bifurcation diagram, we set β=0.1. In E. coli, the typical lifetime of mRNAs (minutes or less) is much smaller than that of proteins (larger than one cell cycle), especially fluorescent reporter proteins (like GFP) [31],[34]. Considering the range of values over which the protein lifetime T _{ P } was varied, we fixed T _{ M }=50 time steps.
The affinity of DNA regulation proteins (transcription factors) for their specific DNA site is also variable, but typically reported values are between 0.51.0 nM [35] and several tens (to few hundreds) of nM [36],[37]. Given a total E. coli cell volume of 3 fl [38], we set k _{ on }=10^{5} per molecule per time step and k _{ off }=10^{3} and k _{ o f f2}=10^{5} per time step. This corresponds to affinities of 55 and 0.55 nM for the free (unbound) and singlybound genes, respectively. Finally, we considered that each of the three genes can be present in the bacterial cell in multiple copies (i.e. ${G}_{i}^{T}>1$). Indeed, even with low copy number plasmids (as used in [9]), each gene is introduced in 34 to 10 copies [39]. Therefore, unless otherwise specified the number of copies per gene type was set to (∀i∈{0,1,2}).
Individualbased simulation of the spatial dynamics
A spatiallyexplicit stochastic individualbased (MonteCarlo) simulator of the set of biochemical reactions given by (1) was implemented as a latticebased algorithm programmed in C. Boundary conditions were periodic. Each of the 3 gene types were present in G _{ T } copies.
Spatial configuration of the genes
The initial locations of the 3G _{ T } genes were set according to one of the following three scenarios (Figure 2):
In the clustered configuration, we first positioned G _{ T } internal boxes of linear size r on the lattice. One copy of each gene i {0,1,2} per box was then positioned at a randomly chosen location inside the box.
In the segregated configuration, we first positioned 3 internal boxes of linear size r on the lattice. Each of the 3 boxes received all of the G _{ T } copies of gene 0,1 or 2, positioned at randomly chosen locations inside the box.
When the box size r converges toward the lattice size R, the clustered and segregated configurations are closing to the uniform one. The ratio r/R can thus be used to quantify the amount of demixing (segregation or clustering): vanishing demixing for r/R 1 and strong demixing for r/R 0. The thus defined gene locations were kept constant during the simulation, i.e. the genes were immobile.
Simulation algorithm
 1.
Each free (unbound) gene copy transcribes a new mRNA molecule M _{ i } with probability α.
 2.
Every mRNA molecule M _{ i } can be degraded (with probability 1/T _{ M }) or can translate a new protein P _{ i } (with probability β≤11/τ _{ M }). If not degraded, each mRNA then undergoes a random walk step (see below).
 3.
Each free (unbound) protein is degraded with probability 1/τ _{ P } or undergoes a random walk step (with probability 1 1/τ _{ P }).
 4.
Whenever a free protein P _{ i } shares its lattice site with a free target gene or singlybound target gene , binding occurs with probability o kon.
 5.
Bound geneprotein complexes that were formed before the current time step can unbind depending on their occupancy status, i.e. with probability o koff or o koff2, if singly or doublybound, respectively.
At each random walk step, the walker changes its current location with probability p _{move}, moving to one of its 4 (2D) or 6 (3D) randomlychosen nearest neighbors (with uniform probability). This corresponds to diffusion coefficient D=p _{move}/4 (2D) or D=p _{move}/6 (l a t t i c e s p a c i n g)^{2}/time step (3D). Note that all diffusive molecules (mRNA and proteins) have identical diffusion coefficient.
Parameter values and numerics
The internode distance of the space lattice was set to △ x=1 (arbitrary units, a.u.) and the lattice size was set to R×R=400×400 a.u.^{2} (2D) or R×R×R=60×60×60 a.u.^{3} (3D). Time was expressed in numbers of MonteCarlo (MC) time steps. Regarding the values of the other parameters, our goal here was to compare with the dynamics predicted by mass action kinetics. Therefore, unless otherwise specified, the standard set of parameters in the individualbased simulations was taken identical to the mass action kinetics model above, i.e. G _{ T }=5 copies, τ _{ M }=50 MC time steps, τ _{ P }=10^{3} MC time steps; probability rates in (MC time steps)^{1}: α=β=0.10, k off=10^{3} and k off2=10^{5}. Note that the value of k on for individualbased simulations is expected to differ from the value used in mass action law kinetics. Indeed in the former, k on is a reaction probability rate upon reactant encounter in space, while in the latter k on is a classical reaction rate constant that, in addition to the reaction probability upon reactant encounter also accounts for reactant encounter probability. Taking the size of the reaction volume into account (400×400 or 60×60×60) we estimated that k on=1 per MC time step per encounter in individualbased simulations was comparable with the value used in the mass action kinetics model. Finally, the standard value of the movement probability was set to p _{move}=1.0, yielding diffusion coefficients D=0.250 or 0.167 a.u.^{2}/MC time step in 2D and 3D, respectively.
Quantifying gene expression dynamics
where $\mathbb{\mathbb{E}}\left[\right]$ denotes expectation over the time points, is the mean of the time series over the window, ${\sigma}_{{P}_{i}}^{2}$ its variance, j∈[0,ll] and the lag time τ∈[0,l1]. We then averaged the A C F(τ) values across each segment of the time series. The first zero crossing of the thus averaged ACF (FZCA) is the smallest value τ _{0} such that A C F(τ _{0})=0. One FZCA was computed for each protein time series resulting from the individualbased simulation. For each simulation run, we retained the minimal value among the three estimated values (one per protein type). Note that the first 2×10^{6} MC time steps of each time series were rejected to discard transient behaviors
Results
Here, we study a generic transcriptional circuit consisting of three genes G _{ i }, i∈{0,1,2} forming a transcriptional ring network of repression: G _{0} represses the expression of G _{1}, G _{1} represses G _{2} and G _{2} closes the loop by repressing G _{0} (Figure 1A). Such repressionbased transcriptional circuits are sometimes referred to as 'repressilator" circuits in the literature [9], because circuits made of odd numbers of negative interactions are, generically, potential spontaneous oscillators [6],[7]. In agreement with the model proposed by Jacob and Monod [1], we assume that each gene G _{ i } is active in the absence of bound protein, thus producing mRNA transcripts M _{ i } (at rate α) then the corresponding protein P _{ i } (rate β). Each gene G _{ i } can bind up to two copies of its repressing protein, with cooperative binding. Once bound, genes stop transcribing mRNAs. mRNAs and proteins have finite lifetimes τ _{ M } and τ _{ P } respectively. See Methods for further details, in particular the set of elementary reaction schemes (1) that describes the system.
In the following, we study this generic repressionbased transcriptional circuit with two modelling approaches: the traditional mass action kinetics, that assume perfectlystirred conditions (infinitely fast mixing) thus neglecting the effects of spatial fluctuations and stochastic individualbased (MonteCarlo) simulations that explicitly take into account spatial fluctuations.
Mass action kinetics predict spontaneous oscillations
The traditional approach in biochemistry to model the kinetics of such generic repressionbased transcriptional circuit is based on the theory of elementary chemical reaction kinetics, usually referred to as "Mass action kinetics" [10]. Mass action kinetics are meanfield approximations assuming that the reaction medium is dilute, perfectlystirred and spatially homogeneous. Under those assumptions, one considers that the local fluctuations of reactant concentration (induced by e.g. the reaction itself) can be neglected and replaces local reactant concentrations by their average values over a large spatial domain (usually the whole reaction volume). For the generic repressionbased transcriptional circuit studied here, mass action kinetics yield the system of 3×4 coupled Ordinary Differential Equations (ODEs) shown as (2).
Figure 1BC shows a bifurcation analysis of these equations. In the twodimensional parameter space defined by the transcription rate α and the protein lifetime τ _{ P }, the system presents two regions delineated by two Hopf bifurcation branches. Outside the region enclosed by the Hopf bifurcation branches, (2) has a single, stable steadystate (the white region in Figure 1B). All the reactants are therefore predicted to converge at long times to stationary values. Inside the grey region of Figure 1B, this steadystate loses its stability and this stability loss is accompanied by the birth of a stable limit cycle (Figure 1C). The points in the parameter space where the steadystate changes stability and the limit cycle appears constitute the Hopf Bifurcation branches. The dynamics in the region delimited by the Hopf bifurcation branches thus consists in spontaneous oscillations where all the mRNA and protein species oscillate in time (Figure 1D). Because of the cyclic cooperative repression, protein oscillations are twobytwo antisynchronized: protein P _{ i } reaches its peak values when P _{ i1} is minimal.
Individualbased simulations predict strong dependence on the spatial locations of the genes
Massaction kinetics, albeit widely used, is actually just one methodology among several others to model the dynamics of the reaction set (1). To evaluate the possible effects of the spatial extension of these reactions, one has to employ alternative modelling methodology. Spatially explicit individualbased simulations (MonteCarlo simulations, see Methods) explicitly describe the diffusive motion of each individual mRNA and protein molecules and emulate the reaction steps of (1) as stochastic processes whenever the concerned reactants encounter in space along their respective random walk (diffusion). Just like for the mass action kinetics, each gene type i∈{0,1,2} is present in G _{ T } copies. In agreement with previous models of the influence of diffusion on gene expression [40],[41], we assumed that the amplitude and speed of the gene displacements in space can be neglected compared to the mRNA and protein, so in our simulations, the genes are immobile and only the mRNAs and proteins move via Brownian motion (with identical diffusion coefficients).

The uniform configuration corresponds to a perfect mixing of the genes: all 3G _{ T } gene copies are positioned at independent randomlychosen locations.

The clustered configuration corresponds to a first case of demixing: G _{ T } gene triplets composed of one G _{0}, one G _{1} and one G _{2} gene are restricted inside nonoverlapping subregions of space

The segregated configuration illustrates a different case of demixing: the 3 gene types i∈{0,1,2} are segregated into 3 subregions of space, each subregion containing all the G _{ T } copies of a given type.
The degree of demixing of the genes can be continuously adjusted by setting the ratio r/R between the size of the local segregation subregions r and that of the total reaction space R. Demixing is strong for r/R→0 (either for the clustered or segregated scenario) but disappears when r/R→1.
The results of the 2D individualbased stochastic simulations for the segregated gene configuration with large demixing however show a strikingly different behavior (Figure 3C). After a transient period, protein levels do not show evidence of oscillations any longer but reach a stationary level around which they fluctuate. In some cases, the stationary regime for one of the proteins can even correspond to a stably and completely repressed state (G _{2} in the example shown in Figure 3C). As a result, the dynamics consists of two proteins (P _{0} and P _{1} in the case shown in Figure 3C) fluctuating around a stationary level, while a third one (P _{2} in Figure 3C) has vanished. We insist here on the fact that the only difference between the oscillatory regime in Figure 3AB and the stationary one in Figure 3C are the spatial locations of the 5×3 immobile genes. All the other parameters (rate constants, species densities/concentrations, lifetimes...) are identical. This is a major result of our paper: changing a purely spatial characteristics (here the positions of the genes) is enough to change the system dynamics in a qualitative manner (from oscillatory to stationary). Since they do not include spatial characteristics, mass action kinetics cannot account for this kind of global change of dynamics. Even worse, the dynamics illustrated in Figure 3C cannot be predicted at all by mass action kinetics. Indeed, according to mass action kinetics ((2)), the only reachable regime consists in a stationary state in which the stationary quantities of protein P _{ i } are identical ∀i∈{0,1,2} (since the kinetic parameters are identical for the 3 genes/mRNA/protein types). In other words, the only accessible stationary state according to mass action kinetics should have P _{0}(∞)=P _{1}(∞)=P _{2}(∞), whereas Figure 3C evidences the existence of a stationary regime where P _{0}(∞)=P _{1}(∞)≫0 and P _{2}(∞)≈0.
To quantify the difference between those two regimes, we found that the first zero crossing of the ACF is a very good quantifier. Figure 3D shows the ACF for the segregated case of Figure 3C. This ACF is typical of stochastic time series fluctuating around a stationary mean: it decreases very rapidly with the autocorrelation delay and is roughly devoid of subsequent oscillations. It is very easy to distinguish it from the ACF of the uniform and clustered configuration that are typical of oscillating regimes, with large period. In particular the first zero crossing (FZCA) is found much smaller in the stationary segregated configuration (0.041×10^{6} MC time steps) than in the oscillatory uniform or clustered ones (>0.2×10^{6} MC time steps). In the following, the FZCA will thus be used as a quantifier to distinguish between stationary regimes (FZCA on the order of 10 thousands MC time steps) and oscillating ones (FZCA of several hundred thousands MC time steps).
A bifurcation based on the spatial localization of the genes
The segregated and clustered configurations illustrated in Figure 3 above correspond to high demixing (i.e. r/R<0.010). We then investigated how the observed effects depend on the degree of demixing.
Here again the segregated scenario yields a very different picture (Figure 4). For large enough mixing (i.e. $r/R\gtrsim 0.5$), the dynamics remain oscillatory, with FZCA values that are indistinguishable from the uniform case. For r/R<0.5, though, the average FZCA progressively decreases with increasing demixing, switching from oscillatory values (>0.1 million MC time steps) to values typical of stochastic fluctuations around a stationary state (<0.05 million MC time steps). The corresponding curve actually describes a bifurcation: as demixing crossovers the critical value (r/R)_{ c }≈0.5, the dynamics undergoes a global (qualitative) change from oscillatory to a stable stationary state. However, the bifurcation parameter here is not a kinetic parameter or density parameter, as usually the case, but a parameter related to the spatial locations of the genes. We refer to this behavior as a spaceinduced bifurcation.
Finally, Figure 5B shows the FZCA values obtained when the copy number of each gene (G _{ T }) is varied. In the spatial configurations with large demixing (clustered or segregated configurations), decreasing the number of genes does not qualitatively change the dynamics (see Additional file 1: Figure S1 b ^{"} and c ^{"} in the Supporting Material): even with a single copy of each gene type, the segregated configuration keeps stationary dynamics whereas the clustered one maintains its oscillatory regime (albeit with modified waveform). The behavior with the uniform configuration is more complex. By definition, with a unique copy of each gene (G _{ T }=1), the uniform and segregated spatial configurations are identical. One the other hand, with standard parameters (i.e. G _{ T }=5 copies of each gene), the above results show that the dynamics of the uniform configuration is essentially identical to the clustered configuration. Therefore, in the uniform configuration, one expects a transition from the oscillatory to the stationary regimes when the number of gene copy decreases from 5 to 1. Figure 5B shows that this transition occurs between G _{ T }=2 and G _{ T }=3. For G _{ T }≤2 copies the uniform spatial configuration display the stationary dynamics typical of the segregated configuration (Additional file 2 : Figure S2) whereas for G _{ T }≥3, the dynamics is oscillatory, similar to the clustered configuration.
Taken together, those results show that with spatially explicit dynamics, the system still displays bifurcations when the kinetic parameters are varied. Spatial parameters thus bring an additional dimension to the bifurcation diagrams, in addition to the kinetic ones.
Spaceinduced bifurcation in three dimensions
In the above results, the movement of the reactants in the stochastic individualbased simulations occurred along a twodimensional spatial domain. Latticebased simulations of Brownian diffusion in two dimensions are less demanding in terms of computation cost than in three dimensions, thus permitting exploration of the parameter space with reasonable accuracy and sampling. However, compared to two dimensions, Brownian diffusion in three dimensions has fundamentally different properties regarding how space is explored (compactness of the random walk) [42],[43]. In this section, we show that the occurrence of spaceinduced bifurcation is robust to those changes in compacity and is preserved in three dimensional systems.
When the diffusion coefficient decreases, the dynamics in those 3D simulations exhibits two phases: for intermediate diffusion coefficients, i.e. $D\gtrsim 1{0}^{2}$ a.u. ^{2}/MC time step (Figure 6B) the regime remains oscillatory for all the gene configurations (Figure 6B), but the period of the oscillations increases to recover the values observed for the oscillatory regimes in 2D above (i.e. close to 200,000 time steps, compare Figure 6C3 with Figure 3A). When the diffusion coefficient decreases further (D≲2×10^{2} a.u. ^{2}/MC time step), the behavior depends on the spatial configuration of the genes. With the clustered configuration, the oscillatory regime persists (with low frequency), at least in the limit of the values of D that we tested. With the uniform or segregated configurations however, the dynamics exhibit a qualitative change as D decreases below 2×10^{2}, switching from an oscillatory regime (Figure 6C2) to a stationary one where all species fluctuate around a constant steady state (Figure 6C1). Therefore, this is a further illustration that spaceinduced bifurcations are also observed in 3D when the spatial range decreases. In agreement with our 2D results above, the segregated configuration appears to be more sensitive than the uniform one: the oscillatorystationary bifurcation seems to necessitate less reduction of the diffusion coefficient in the segregated configuration than in the uniform one. This observation is a further example that the dynamic regime in our system is strongly controlled by the gene positions in space.
Discussion
Our goal in this study was to investigate whether space can influence the temporal dynamics of repressilator circuits i.e. 3gene repressionbased transcriptional regulatory loops, that exhibit prototypical spontaneous oscillations in certain ranges of kinetic parameters (transcription/translation rate, species lifetimes). We used spatiallyexplicit stochastic individualbased modelling to simulate repressilator circuits in 2D and 3D, with various degrees of demixing of the gene positions in space. Our main finding is that variations of some spatial parameters (degree of demixing of the genes, spatial range of the mRNA and proteins) can have dramatic effects on the system dynamics, switching it from the spontaneous oscillatory regime to a stationary regime where each species fluctuates around constant values. This effect is similar to the bifurcations along the kinetic parameters that are usually evoked to explain the appearance of the oscillatory regime in those systems. We thus referred to it as a spaceinduced bifurcation. Our study therefore indicates that spatial parameters should be considered as additional bifurcation dimensions to the kinetic parameters to predict the dynamics of those systems. This conclusion strongly supports the idea that the spatial organization of the molecular actors of transcriptional networks is crucial for the dynamics of gene expression.
The possibility that the positions in space of the elements of intracellular biochemical systems may control their dynamics has already been suggested in previous works (see e.g. [24][26],[44][46]). In the specific case of the dynamics of gene expression, the study of the influence of space on the dynamics of the expression of an isolated gene has recently started to be explored with computational or theoretical approaches. The main conclusion was that diffusion of mRNAs and proteins and the spatial correlations created by the coupling with reaction, can strongly increase the fluctuations of gene expression, i.e. noise [40],[41],[47]. Our results is a significant advance in this problem since we consider gene networks (albeit small ones), not a single isolated gene and we show that the alterations of the dynamics due to spatial parameters in those systems may be qualitative: beyond changing the mean value or the fluctuations, they may even alter the global regime of the dynamics (stationary vs oscillatory).
From our simulation results, two major experimental predictions emerge concerning the dynamics of 3gene repressionbased transcriptional ring networks. If each of the 3 genes is present in a single copy (see Figure 5B, G _{ T }=1), we expect the dynamics to be stationary in most cases, except if each of the repressor gene is located very close to its repressed target gene (i.e. the clustered scenario). In other words, we predict that spontaneous oscillations may be difficult to obtain in networks with a single copy of each gene. Moreover, our study suggests that the spatial range of the proteins and mRNAs, i.e. their Kuramoto length (the typical distance they travel before degradation), is a major determinant of the system dynamics. The precise location of the genes is likely to control the dynamics when the spatial range (i.e. the product of the protein or mRNA lifetime and its diffusion coefficient) is not too large. For large spatial ranges (i.e. Figure 5A, 1/τ _{ M }=0.005 or Figure 6B, D=0.167), the expected regime agrees with the prediction of mass action kinetics, i.e. spontaneous oscillations with the kinetic parameters we used. Therefore, in vivo, spaceinduced bifurcations should be significant in systems where the lifetimes and/or the diffusion coefficient of the mRNAs and/or proteins are limited.
Estimating the spatial range of intracellular proteins in vivo is still a challenge for experimental biology. Based on the measured diffusion coefficient of Fus3 MAP kinase in the yeast [48] or that of GFP in E. coli[49], Cottrell et al. [41] obtained coarse estimates leading to the conclusion that cytoplasmic proteins should have spatial ranges that are much larger than the cell itself. This leads to the widespread opinion that the spatial distribution of cytoplasmic proteins in the cell is uniform (wellmixed). However, even in bacteria like E. coli, this idea is questionable. First, consistent recent experimental evidence suggest that many proteins adopt localized distribution inside the cell, in opposition to a wellmixed situation [50][53]. The molecular actors of gene expression in C. crescentus or E.coli may even present a very specific spatiallyorganized intracellular structure [20], with, in particular, chromosomallyexpressed mRNAs that exhibit very low diffusion coefficients [19]. Another recent and very symbolic example is LacI, the repressor of the lac operon. Direct measurements of the diffusion coefficient of LacI in living E.coli cells indicate rapid diffusion (D of the order or 0.1–1 μm ^{2}/s) and consequently suggest a very large spatial range [54]. However, recent measurement of the steadystate distribution of LacI inside living E.coli revealed that depending on the location of the lacI gene on the bacterial chromosome, the distribution of the LacI protein can either be homogeneous or highly inhomogeneous and mostly localized [55]. In fact, the diffusion of proteins and mRNAs in living cells, even in simple and small cells like E.coli is far from a simple Brownian motion. This may lead to violations of the hypotheses that underly coarse estimations of the spatial range. First, for macromolecules that interact with DNA, the diffusion process itself is composite (facilitated diffusion) because part of it occurs in 3D in the bulk and part of it as a restricted sliding movement along the DNA molecule [56][58]. Moreover the bacterial cytoplasm itself does not present the usual properties of a liquid. It is a complex, extremely crowded and dense medium that strongly affects protein and mRNA mobility in a spatially nonhomogeneous way [11][15]. How exactly those complex properties arise from the intracellular elements and how exactly they alter molecule mobility is still unclear [15], but theoretical arguments indicate that such macromolecular crowding may strongly affect the dynamics of gene expression [47],[59]. We conclude that, in spite of the rather large values measured for the diffusion coefficient of some proteins or mRNAs, the physical nature of the cytoplasm is likely to considerably restrict the spatial range of proteins and mRNAs. Therefore, one would generically expect that spaceinduced bifurcations may be significant in living cells.
One possibility that we have left unexplored in the present work is that of an inverse transition: starting from a region in the parameter space where mass action kinetics predict a stable steady state (i.e. the white region in Figure 1B), can oscillations be induced by manipulation of some spatial parameter? We first stress that the observed shunting down of oscillations resulting from spatial segregation is not exactly a return to the steady state of the mass action kinetics, since in the latter, all three proteins are predicted to be expressed in equal amounts whereas in the former, one of the proteins is often expressed at a much lower level than the others. In previous works [24][26], we have showed that in the absence of bifurcation in the mass action kinetics, the unique equilibrium remains stable, globally and asymptotically. However, when a bifurcation does exist in the mass action kinetics approximation, the question remains open, especially when the parameters are very close to the bifurcation of the mass action kinetics system.
Conclusion
Our model for gene expression is limited to the main processes implied in gene expression and is limited to a 3gene network. We think that the basic mechanisms underlying spaceinduced bifurcations are simple enough that they should still be effective in other systems, in particular, in transcriptional regulation networks comprising larger numbers of gene types or positive regulation (activators). Likewise, we expect our conclusions to be robust to increasing molecular details and more realistic modelling of each subprocess (RNA polymerase, initiation, elongation, termination, ribosomes...). Therefore, our results bring a strong support to the view that the spatial organization of the molecular actors of transcriptional networks is crucial for the dynamics of gene expression. In a synthetic biology perspective, they suggest that the spatial localisation of the synthetic genes should be factored out in the global strategy used to shape the dynamics of the synthetic network to be inserted in the cell.
Additional files
Declarations
Authors’ Affiliations
References
 Jacob F, Monod J: Genetic regulatory mechanisms in the synthesis of proteins. J Mol Biol. 1961, 3 (3): 3180356.View ArticlePubMedGoogle Scholar
 Gann A: Jacob and monod: from operons to evodevo. Curr Biol. 2010, 20 (17): 718723.View ArticleGoogle Scholar
 Kwak H: Control of transcriptional elongation. Annu Rev Genet. 2013, 47 (1): 483508. PMID24050178,PubMed CentralView ArticlePubMedGoogle Scholar
 Dogini DB, Pascoal VDB, Avansini SH, Vieira AS, Pereira TC, LopesCendes I: The new world of rnas. Genet Mol Biol. 2014, 37 (1 Suppl): 285293.PubMed CentralView ArticlePubMedGoogle Scholar
 Hirsch MW: Convergent activation dynamics in continuous time networks. Neural Netw. 1989, 2 (5): 331349.View ArticleGoogle Scholar
 Thomas R, Thieffry D, Kaufman M: Dynamical behaviour of biological regulatory networks. i. biological role of feedback loops and practical use of the concept of the loopcharacteristic state. Bull Math Biol. 1995, 57 (2): 247276.View ArticlePubMedGoogle Scholar
 Gouzé JL: Positive and negative circuits in dynamical systems. J Biol Syst. 1998, 06 (01): 1115.View ArticleGoogle Scholar
 Gardner TS, Cantor CR, Collins JJ: Construction of a genetic toggle switch in Escherichia coli. Nature. 2000, 403: 339342.View ArticlePubMedGoogle Scholar
 Elowitz MB, Leibler S: A synthetic oscillatory network of transcriptional regulators. Nature. 2000, 403: 335338.View ArticlePubMedGoogle Scholar
 Chellaboina V, Bhat S, Haddad MM, Bernstein DS: Modeling and analysis of massaction kinetics. IEEE Contr Syst. 2009, 29 (4): 6078.View ArticleGoogle Scholar
 Golding I, Cox EC: Physical nature of bacterial cytoplasm. Phys Rev Lett. 2006, 96 (9): 098102View ArticlePubMedGoogle Scholar
 Weber SC, Spakowitz AJ, Theriot JA: Bacterial chromosomal loci move subdiffusively through a viscoelastic cytoplasm. Phys Rev Lett. 2010, 104 (23): 238102View ArticlePubMedGoogle Scholar
 English BP, Hauryliuk V, Sanamrad A, Tankov S, Dekker NH, Elf J: Singlemolecule investigations of the stringent response machinery in living bacterial cells. Proc Natl Acad Sci U S A. 2011, 108 (31): 365373.View ArticleGoogle Scholar
 Coquel A, Jacob J, Primet M, Demarez A, Dimiccoli M, Julou T, Moisan L, Lindner A, Berry H: Localization of protein aggregation in Escherichia coli is governed by diffusion and nucleoid macromolecular crowding effect. PLoS Comput Biol. 2013, 9 (4): 1003038View ArticleGoogle Scholar
 Parry BR, Surovtsev IV, Cabeen MT, O’Hern CS, Dufresne ER, JacobsWagner C: The bacterial cytoplasm has glasslike properties and is fluidized by metabolic activity. Cell. 2014, 156 (12): 183194.PubMed CentralView ArticlePubMedGoogle Scholar
 Cremer T, Cremer M: Chromosome territories. Cold Spring Harb Perspect Biol. 2010, 2 (3): 003889View ArticleGoogle Scholar
 Le TBK, Imakaev MV, Mirny LA, Laub MT: Highresolution mapping of the spatial organization of a bacterial chromosome. Science. 2013, 342 (6159): 731734.PubMed CentralView ArticlePubMedGoogle Scholar
 Képès F: Periodic transcriptional organization of the e.coli genome. J Mol Biol. 2004, 340 (5): 957964.View ArticlePubMedGoogle Scholar
 Llopis PM, Jackson AF, Sliusarenko O, Surovtsev I, Heinritz J, Emonet T, JacobsWagner C: Spatial organization of the flow of genetic information in bacteria. Nature. 2010, 466 (7302): 7781.PubMed CentralView ArticleGoogle Scholar
 Weng X, Xiao J: Spatial organization of transcription in bacterial cells. Trends in Genetics. 2014, 30 (7): 287297.View ArticlePubMedGoogle Scholar
 Norris V, Turnock G, Sigee D: The escherichia coli enzoskeleton. Mol Microbiol. 1996, 19 (2): 197204.View ArticlePubMedGoogle Scholar
 Vélot C, Srere PA: Reversible transdominant inhibition of a metabolic pathway. in vivo evidence of interaction between two sequential tricarboxylic acid cycle enzymes in yeast. J Biol Chem. 2000, 275 (17): 1292633.View ArticlePubMedGoogle Scholar
 Amar P, Legent G, Thellier M, Ripoll C, Bernot G, Nystrom T, Saier MHJr, Norris V: A stochastic automaton shows how enzyme assemblies may contribute to metabolic efficiency. BMC Syst Biol. 2008, 2: 27PubMed CentralView ArticlePubMedGoogle Scholar
 Caré BR, Soula HA: Impact of receptor clustering on ligand binding. BMC Syst Biol. 2011, 5: 48PubMed CentralView ArticlePubMedGoogle Scholar
 Caré BR, Soula HA: Receptor clustering affects signal transduction at the membrane level in the reactionlimited regime. Phys Rev E. 2013, 87: 012720View ArticleGoogle Scholar
 Soula H, Caré B, Beslon G, Berry H: Anomalous versus sloweddown brownian diffusion in the ligandbinding equilibrium. Biophys J. 2013, 105 (9): 20642073.PubMed CentralView ArticlePubMedGoogle Scholar
 Roberts E, Magis A, Ortiz JO, Baumeister W, LutheySchulten Z: Noise contributions in an inducible genetic switch: a wholecell simulation study. PLoS Comput Biol. 2011, 7 (3): 1002010View ArticleGoogle Scholar
 Karr JR, Sanghvi JC, Macklin DN, Gutschow MV, Jacobs JM, Bolival B, AssadGarcia N, Glass JI, Covert MW: A wholecell computational model predicts phenotype from genotype. Cell. 2012, 150 (2): 389401.PubMed CentralView ArticlePubMedGoogle Scholar
 Thiele I, Fleming RMT, Que R, Bordbar A, Diep D, Palsson BO: Multiscale modeling of metabolism and macromolecular synthesis in e. coli and its application to the evolution of codon usage. PLoS One. 2012, 7 (9): 45635View ArticleGoogle Scholar
 Roberts E, Stone JE, LutheySchulten Z: Lattice microbes: highperformance stochastic simulation method for the reactiondiffusion master equation. J Comput Chem. 2013, 34 (3): 245255.PubMed CentralView ArticlePubMedGoogle Scholar
 Taniguchi Y, Choi PJ, Li GW, Chen H, Babu M, Hearn J, Emili A, Xie XS: Quantifying E. coli proteome and transcriptome with singlemolecule sensitivity in single cells. Science. 2010, 329 (5991): 533538.PubMed CentralView ArticlePubMedGoogle Scholar
 So LH, Ghosh A, Zong C, Sepulveda LA, Segev R, Golding I: General properties of transcriptional time series in Escherichia coli. Nat Genet. 2011, 43 (6): 554560.PubMed CentralView ArticlePubMedGoogle Scholar
 Ishihama Y, Schmidt T, Rappsilber J, Mann M, Hartl FU, Kerner MJ, Frishman D: Protein abundance profiling of the Escherichia coli cytosol. BMC Genomics. 2008, 9: 102PubMed CentralView ArticlePubMedGoogle Scholar
 Andersen JB, Sternberg C, Poulsen LK, Bjorn SP, Givskov M, Molin S: New unstable variants of green fluorescent protein for studies of transient gene expression in bacteria. Appl Environ Microbiol. 1998, 64 (6): 22402246.PubMed CentralPubMedGoogle Scholar
 Oehler S, Eismann ER, Kramer H, MullerHill B: The three operators of the lac operon cooperate in repression. EMBO J. 1990, 9 (4): 973979.PubMed CentralPubMedGoogle Scholar
 Guille MJ, Kneale GG: Methods for the analysis of DNAprotein interactions. Mol Biotechnol. 1997, 8 (1): 3552.View ArticlePubMedGoogle Scholar
 Majka J, Speck C: Analysis of proteinDNA interactions using surface plasmon resonance. Adv Biochem Eng Biotechnol. 2007, 104: 1336.PubMedGoogle Scholar
 Volkmer B, Heinemann M: Conditiondependent cell volume and concentration of Escherichia coli to facilitate data conversion for systems biology modeling. PLoS One. 2011, 6 (7): 23126View ArticleGoogle Scholar
 Lutz R, Bujard H: Independent and tight regulation of transcriptional units in Escherichia coli via the LacR/O, the TetR/O and AraC/I1I2 regulatory elements. Nucleic Acids Res. 1997, 25 (6): 12031210.PubMed CentralView ArticlePubMedGoogle Scholar
 van Zon JS, Morelli MJ, TanaseNicola S, ten Wolde PR: Diffusion of transcription factors can drastically enhance the noise in gene expression. Biophys J. 2006, 91 (12): 43504367.PubMed CentralView ArticlePubMedGoogle Scholar
 Cottrell D, Swain PS, Tupper PF: Stochastic branchingdiffusion models for gene expression. Proc Natl Acad Sci U S A. 2012, 109 (25): 96999704.PubMed CentralView ArticlePubMedGoogle Scholar
 de Gennes PG: Kinetics of diffusioncontroled processes in dense polymer systems. i. nonentangled regimes. J Chem Phys. 1982, 76: 33163321.View ArticleGoogle Scholar
 Benichou O, Chevalier C, Klafter J, Meyer B, Voituriez R: Geometrycontrolled kinetics. Nat Chem. 2010, 2 (6): 472477.View ArticlePubMedGoogle Scholar
 Lipkow K, Andrews SS, Bray D: Simulated diffusion of phosphorylated CheY through the cytoplasm of Escherichia coli. J Bacteriol. 2005, 187 (1): 4553.PubMed CentralView ArticlePubMedGoogle Scholar
 Takahashi K, TanaseNicola S, ten Wolde PR: Spatiotemporal correlations can drastically change the response of a mapk pathway. Proc Natl Acad Sci U S A. 2010, 107 (6): 24732478.PubMed CentralView ArticlePubMedGoogle Scholar
 Maler O, Halasz AM, Lebeltel O, Maler O: Exploring the dynamics of mass action systems. In Proceedings Second International Workshop on Hybrid Systems and Biology (HSB 2013). Edited by Dang T, Piazza C. Taormina, Italy: Electronic Proceedings in Theoretical Computer Science; 2nd September 2013:84–91. ., [http://eptcs.web.cse.unsw.edu.au/paper.cgi?HSB2013.6]
 Pulkkinen O, Metzler R: Distance matters: the impact of gene proximity in bacterial gene regulation. Phys Rev Lett. 2013, 110 (19): 198101View ArticlePubMedGoogle Scholar
 Maeder CI, Hink MA, Kinkhabwala A, Mayr R, Bastiaens PIH, Knop M: Spatial regulation of Fus3 MAP kinase activity through a reactiondiffusion mechanism in yeast pheromone signalling. Nat Cell Biol. 2007, 9 (11): 13191326.View ArticlePubMedGoogle Scholar
 Elowitz MB, Surette MG, Wolf PE, Stock JB, Leibler S: Protein mobility in the cytoplasm of Escherichia coli. J Bacteriol. 1999, 181 (1): 197203.PubMed CentralPubMedGoogle Scholar
 Shapiro L, McAdams HH, Losick R: Why and how bacteria localize proteins. Science. 2009, 326 (5957): 12251228.View ArticlePubMedGoogle Scholar
 Mauriello EMF, Astling DP, Sliusarenko O, Zusman DR: Localization of a bacterial cytoplasmic receptor is dynamic and changes with cellcell contacts. Proc Natl Acad Sci U S A. 2009, 106 (12): 48524857.PubMed CentralView ArticlePubMedGoogle Scholar
 Sliusarenko O, Heinritz J, Emonet T, JacobsWagner C: Highthroughput, subpixel precision analysis of bacterial morphogenesis and intracellular spatiotemporal dynamics. Mol Microbiol. 2011, 80 (3): 612627.PubMed CentralView ArticlePubMedGoogle Scholar
 Landgraf D, Okumus B, Chien P, Baker TA, Paulsson J: Segregation of molecules at cell division reveals native protein localization. Nat Methods. 2012, 9: 480485.PubMed CentralView ArticlePubMedGoogle Scholar
 Elf J, Li GW, Xie XS: Probing transcription factor dynamics at the singlemolecule level in a living cell. Science. 2007, 316 (5828): 11911194.PubMed CentralView ArticlePubMedGoogle Scholar
 Kuhlman TE, Cox EC: Gene location and dna density determine transcription factor distributions in Escherichia coli. Mol Syst Biol. 2012, 8: 610PubMed CentralView ArticlePubMedGoogle Scholar
 Loverdo C, Benichou O, Moreau M, Voituriez R: Enhanced reaction kinetics in biological cells. Nature Phys. 2008, 4: 134137.View ArticleGoogle Scholar
 Benichou O, Chevalier C, Meyer B, Voituriez R: Facilitated diffusion of proteins on chromatin. Phys Rev Lett. 2011, 106 (3): 038102View ArticlePubMedGoogle Scholar
 Kuhlman TE, Cox EC: DNAbindingprotein inhomogeneity in E. coli modeled as biphasic facilitated diffusion. Phys Rev E. 2013, 88 (2): 022701View ArticleGoogle Scholar
 Morelli MJ, Allen RJ: Effects of macromolecular crowding on genetic networks. Biophys J. 2011, 101 (12): 28822891.PubMed CentralView ArticlePubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.