Engineering self-organized criticality in living cells

Complex dynamical fluctuations, from molecular noise within cells, collective intelligence, brain dynamics or computer traffic have been shown to display noisy behaviour consistent with a critical state between order and disorder. Living close to the critical point can have a number of adaptive advantages and it has been conjectured that evolution could select (and even tend to) these critical states. One way of approaching such state is by means of so called self-organized criticality (SOC) where the system poises itself close to the critical point. Is this the case of living cells? It is difficult to test this idea given the enormous dimensionality associated with gene and metabolic webs. In this paper we present an alternative approach: to engineer synthetic gene networks displaying SOC behaviour. This is achieved by exploiting the presence of a saturation (congestion) phenomenon of the ClpXP protein degradation machinery in E. coli cells. Using a feedback design that detects and then reduces ClpXP congestion, a critical motif is built from a two-gene network system, where SOC can be successfully implemented. Both deterministic and stochastic models are used, consistently supporting the presence of criticality in intracellular traffic. The potential implications for both cellular dynamics and designed intracellular noise are discussed.


I. INTRODUCTION
In order to adapt to environmental challenges, biological systems exhibit a diverse array of response mechanisms grounded in sensors and actuators as well as information-processing units. Adaptive responses require dynamical features that combine low energetic costs along with fast changes to efficiently respond to environmental changes. Flocks of birds and fish schools widely fluctuate in time but rapidly reorganize when a perturbation (such as the presence of a predator) occurs. Within cells, noise was early identified as playing multiple roles affecting cell fate, population heterogeneity, signal amplification or response to stress [1][2][3]. Noise is both an inevitable outcome of stochastic molecular interactions and an essential ingredient in decision making [4].
It has been shown that many complex systems seem to be poised close to so called critical points separating ordered from disordered states [5][6][7][8]. In a nutshell, both living and non-living systems organize at the boundary separating regular (predictable) from random (disordered) behaviours. At this point, complex dynamics with scale-invariant properties emerge [9,10]. If s defines the total activity in one given event, such as number of firing neurons [11][12][13][14], gene expression [15][16][17][18], number of active ants in a colony [19,20], critical epidemic bursts [21] or the size of traffic jams Critical points can be reached by fine-tuning a given "control" (or bifurcation) parameter η ( Fig. 1a-b). This parameter (such as density of particles, temperature or reaction rate) directly influences the system's state, as described by the order parameter S (system's activity, for example). The way to criticality based on tuning key parameters is well illustrated by enzymatic queueing processes [26]. These authors used the framework of queueing theory (QTH) to study the dynamics of different proteins (the 'customers' in QTH) that are processed by a downstream set of enzymes that play the role of 'servers'. Specifically, they considered the native E. coli protease, ClpXP. The key concept here is that the protease complexes are a limited resource (hereafter we will thus consider the amount of ClpXP as a constant) that can only 'process' (degrade) a limited number of incoming proteins. In Fig. 1a we provide a basic diagram considering a protein σ being expressed at some given rate η. If the rate of protein production is low (queues are short), degradation is efficient since the proteases can process all incoming σ units (free phase). If production is too high, a long queue of molecules 'waiting' to be processed will be present (congested phase). The two regimes are separated by a critical point where an optimal balance is reached, along with wide fluctuations in concentrations [26].
The critical point is a rather unique one. Can these systems poise themselves into critical states with-FIG. 1: Paths to intracellular criticality. Tunable critical dynamics can be found in simple genetic circuits (a) where a given gene is constitutively expressed into a protein σ that decays and is also actively degraded by cell proteolytic machinery (ClpXP). By tuning expression rate η (d), a critical rate ηc is found to separate a phase of efficient degradation from another (light gray) involving congestion. In (b) the thick line indicates that few proteins are found for η < ηc (the proteolytic machinery efficiently degrades it) while it accumulates on the right side, due to congestion (ClpXP fails to degrade all the incoming proteins). An alternative, non-tuned path is self-organised criticality (SOC) which can emerge is provided by the sandpile (c). As grains of sand are slowly added at a rate η, the angle of the pile θ grows and only small avalanches will be observed. However, as the critical (maximum) θc is reached, avalanches of all sizes take place, reducing θ. The feedback between the order parameter S and the control parameter θ is summarised in (d). To facilitate the conditions enabling SOC, a two-gene circuit with negative feedback (e) allows mapping the sandpile feedback diagram (f). Here, both proteins compete for ClpXP (higher levels of σ1 also implies high values of σ2) and a repression feedback is mediated by σ2σ2 (the Lac repressor dimer) with σ1 and σ2 acting as order and control parameters, respectively. out fine tuning? An alternative mechanism to reach critical states is provided by self-organized criticality (SOC) [27][28][29]. In this case, control and order parameters interact in such a way that the system spontaneously self-organizes into a critical state [30,31]. The canonical example of SOC is the critical sandpile (Fig. 1c). By slowly adding grains of sand to the pile (at a rate η), its slope θ increases. At the beginning, only a few grains will fall down but the number s of grains in an avalanche rapidly grows as the angle of repose θ c is approached. Once we have θ = θ c , the interaction between θ and sand avalanches (the order parameter) will keep the system at criticality [29]. This is summarized in Fig. 1d where the nature of the feedback between control and order parameters is sketched. Is a SOC mechanism a possible way of generating critical gene expression dynamics? In this paper we propose a novel approach to create SOC dynamics by engineering the interaction between order and control parameters in a simple two-gene network design.

II. RESULTS
A. Two-gene SOC motif model: deterministic and stochastic dynamics The importance of the queueing dynamics in the enzymatic processing is illustrated by the E. coli stress response to starvation, which is triggered by an excess of mistranslated proteins. Stress can cause a significant increase in the concentration of aberrant proteins, which must be degraded. When such an overload occurs, the concentration of the sigma factor (the master stress regulator) builds up, eventually triggering the stress response [34]. Recent theoretical work also suggests that queueing could be adaptive in parallel enzymatic networks when the input flux of substrates is balanced by the maximum processing capacity of the network [35]. Here we go a step further and show how a simple SOC circuit can be actually engineered in vivo.
Our goal in this work is to define the basic design principle to build a "genetic sandpile" system that can capture the feedback structure shown in Fig. 1e. First of all, consider the simple, two-gene network circuit shown in Fig. 1e. Two proteins σ 1 and σ 2 resulting from their expression will be used as the building blocks for the order and control parameters, respectively, thus implementing a SOC feedback loop (Fig. 1f). The aim is to exploit the topology of the gene-gene interaction in such a way that the system can detect the degree of congestion of the ClpXP system by using σ 2 as a sensor of the σ 1 levels. Our control protein σ 2 can form dimers, i. e. σ 2 + σ 2 → σ 2 σ 2 and this dimeric forms act as inhibitors (see Methods, eqs. (1-3)). If congestion occurs, the abundance of σ 2 increases and its negative feedback effects also do so. A standard Hill function will be used to model the dimers as transcriptional repressors. The construction of our circuit involves two steps: (a) engineering the "critical" motif and (b) adjusting the protein production levels. This might seem contradictory with the "self-organized" description of SOC, but an example of why this is required is given by rice piles [36] which exhibit SOC but for some given grain aspect ratios.
Along with the topology of the SOC motif, a separation of scales is known to be a characteristic of SOC dynamics [31]. While the control parameter has a slow dynamics (the angle of the sandpile) the system's response (the avalanche time scale) is fast. In our model, two additional parameters are used to favor the presence of criticality. These are the promoter efficiency for σ 2 , labelled η 2 ; and an extra inhibition acting on the repressor σ 2 σ 2 indicated as µ in Fig. 1e. We need to remember that degradation (and other dissipative events) affects σ 2 and thus a minimal concentration of this sensor is needed in order to effectively detect congested states. On the other hand, in order to experimentally validate our model, we need to tune the strength of the feedback (required to trigger a rapid decay of the intracellular concentration of σ 1 ). By tuning these two parameters, we include both the non-SOC design based on queueing as a special case [26] (see Section I in the SM) and a mechanism to achieve the SOC state.
In Fig. 2 we summarise the behaviour of the twogene system ( Fig. 1e; see Methods and Section II in the SM for mathematical details) using both deterministic (Section II.A, SM) and stochastic (Section II.B, SM) dynamics. A unique stable equilibrium (σ eq = (σ 1,eq , σ 2,eq ), indicated with a solid black circle in the (σ 1 , σ 2 ) phase portraits of Fig. 2) is found, with a characteristic structure of the orbits in the phase space, as shown in Fig. 2a for the unregulated domain (here η 2 = 10 −3 ). As the expression rate η 2 of the control parameter increases close to η 2 ∼ 10 −2 , it is easy to see the presence of a slow-fast dynamics in the distinct structure of the vector fields consistently with the SOC requirement of time scale separation. Here the critical motif allows for large fluctuations in σ 1 to occur (Fig. 2b) as shown by the compression of the trajectories in the phase portrait close to the fixed point, to be compared with the more homogeneous flow displayed in Fig. 2a for η = 10 −3 . The analysis of this system shows that, once close to the equilibrium point, small changes in the control σ 2 trigger marked population spikes in σ 1 (see Fig. S4, SM). Larger values for η 2 (Fig. 2c) do not exhibit such a time scales separation. The analytic and numerical investigation of the eigenvalues of the fixed point to study its stability properties reveal the presence of a maximum in the so called ratio index (Fig. S7, SM), indicating a remarkable change in the vector field of the phase portrait when η 2 ∼ 10 −2 (given our set of fixed parameters indicated in Fig. 2).
To see how this nonlinear flows behave under the presence of intrinsic noise, a stochastic numerical implementation of the two-gene circuit has been carried out using the Gillespie method [32]. In Fig. 2d the coefficient of variation CV = σ 2 1 − σ 1 2 / σ 1 of the generated time series is displayed against η 2 for three values of µ. This coefficient provides a statistical estimate of the variance of the fluctuations and a well-defined maximum is observed when η 2 ∼ 10 −2 . In Fig. 2e we have overlapped several stochastic realisations with the vector field close to the maximum CV (for η 2 = 10 −2 ). The density plot reveals that the stochastic system visits very frequently the fixed point σ eq (orange-red colours in Fig. 2e), but also wanders far away in the lower part of the phase portrait, where the vector field is faster and pushes the stochastic paths far away from the equilibrium point, then returning it back to the deterministic equilibrium. In Fig. 2f the resulting distribution is displayed. Specifically, if P (σ 1 ) indicates the probability distribution of σ 1 expression levels (activity), the cumulative distribution is defined as P > (σ 1 ) = σ1 0 P (σ) dP (σ) and helps smoothing the random noise exhibited by P (σ 1 ). If the original distribution follows a scaling P (σ) ∼ σ −γ 1 , the cumulative one gives P > (σ 1 ) ∼ σ −γ+1 1 . The stochastic model gives a value of γ ∼ 3 ( fig. 2f) estimated from the average of five different runs. The time series associated to this parameter combination is shown in the inset, revealing a characteristic bursting dynamics typical of the SOC state. These results can be compared with the smooth and Gaussian behaviour obtained with the unregulated dynamics setting µ = 0 ( Fig. 2g-h).

B. Engineering a synthetic SOC circuit in E. coli
The theoretical model predicts that the SOC feedback loop defined above (Fig. 1e-f) will display bursting dynamics with fat-tailed activity distributions associated to the σ 1 protein (order parameter). If the expression level (η 2 ) of the control protein σ 2 is large enough, its concentration will act as a congestion sen- The feedback control parameter is (a) η2 = 10 −3 , (b) η2 = 10 −2 and (c) η2 = 0.056, respectively. The nullclines are plotted in red (dσ1/dt = 0) and blue (dσ2/dt = 0). The orbits (shown with black curves) flow towards a single attractor (σeq). The vector field is indicated by unitary arrows, the color of which corresponds to their module (blue for small and red for the heavier). A background color scale is also used showing the arrival time required for each initial condition (σ1(0), σ2(0)).
Here yellow and violet indicate short and long arrival times, respectively. In (b) a combination of slow decays in σ2 and fast responses for σ1 are at work. The stochastic dynamics of the model reveals a maximum in the CV when η2 ∼ 10 −2 , as shown in panel (d) where the colours stand for µ = 0.5 (blue), µ = 1.0 (orange) and µ = 1.5 (green). The relative location of the deterministic flows are indicated by dashed lines. Three different values of the coupling parameter µ are also used to show the robust nature of the maximum, where the SOC motif has been tuned to generate fat-tailed behaviour, as shown in panel (e). In this panel the hot map is plotted on top of the phase space, showing a maximum close to the attractor as well as the fat-tailed scaling behaviour (f) with P>(σ1) ∼ σ −2 1 which gives γ ∼ 3 for P (σ1). Here five different runs are shown along with their average (dark line). By contrast, the flows and hot maps for the non-SOC circuit close to the queuing theory (QTH) transition (g) have a Gaussian pattern (h) with exponential tails as shown by the straight lines in the linear-log insets. The complete distribution is depicted with the upper violin plots (i) for different parameter values (as indicated). The effect on the coupling (µ) can be also observed in the different plotted simulation points. sor and will repress σ 1 , following the SOC motif design. Otherwise, the system will lack the feedback loop. Here we show how can we engineer the genes circuit incorporating some parameters that allow including the non-SOC phase transition as described above including intracellular queueing processes as a special case. Since the time scale of expression changes is comparable to the replication rate of individual cells, no individual time series can be gathered, but instead the collective response of the system will be analyzed to detect the presence of a SOC state. The underlying assumption is thus that we have a colony-level sampling of the dynamical states and the critical dynamics is assessed by looking at the distribution of cell states and the resulting aggregated statistics (as-suming ergodicity).
The explicit experimental design of our SOC circuit implementation is outlined in Fig. 3a-b. The order parameter is encoded by the green fluorescent protein (GFP), and the control parameter acting as the congestion sensor, by the LacI repressor protein, the relative expression of which will be estimated by means of the expression of the red fluorescent protein (RFP). In both cases we use the unstable variants GFP-lva and RFP-lva, respectively. The construct expresses GFP under the pLacI promoter (η 1 ) while the LacI repressor and RFP reporter protein are under the pBAD promoter, with non leaky tight regulation and highlevel expression inducible by Arabinose (η 2 ) [52][53][54]. All three proteins of the circuit are tagged with lva sequence to be degradable by the ClpXP proteolytic complex.
ClpXP is responsible for degrading proteins carrying the SsrA or YbaQ degron sequences, reducing the half-life of a tagged protein from hours to minutes. In an exponentially growing E. coli culture (Optical Density (OD) OD 660 from 0 to 2), the endogenous levels of ClpX and ClpP are constant and involve around 100 ClpXP molecules, which can degrade at least 10 5 molecules of GFPssrA per cell per replication cycle. However, due to the limited number of ClpXP protease complexes, the degrading capacity of ssrA-tagged proteins can be easily saturated by over production of a synthetic tagged protein [49][50][51].
The non-induced circuit depicted in Fig. 3a is thus a particular instance of our more complex motif (Fig. 3b) that would correspond in our case to the presence of endogenous LacI, with the GFP-lva expression being repressed. The presence of Arabinose leads to a strong expression of the repressor LacIlva (our control parameter) and the reporter protein RFP-lva. Only when high levels of the -lva tagged proteins are reached, and the degradation machinery is saturated, there is enough LacI-lva to repress the The gene construct SOC design in its induced state. The Arabinose (Ara) concentration controls the production of degradable LacI repressor and RFP, both, under the pBAD promoter. This design allows to use the red fluorescence as a proxy of LacI concentration. (c) Overlapped bright field and fluorescence images of bacterial culture induced with Ara (100 mM) and IPTG (10µM). Here, yellow bacteria expresses both GFP and RFP. (d) Flow cytometry dot plots (green channel vs red channel) of E. coli cultures exposed to different concentrations of IPTG and Ara. In the non-induced circuit (a), without Ara, the transition from non-congested proteolytic machinery phase (i e free, SUBCRIT) to congested phase (SPCRIT) depends on the tunable GFP-lva production. As IPTG increases, hence the de-repression of GFP expression, ClpXP is not able to degrade the excess of GFP and cells are mostly with high emission in green. When the circuit is induced (b), Ara triggers the expression of the LacI repressor and the RFP reporter, that are also degraded by the ClpXP complex. The increase of tagged proteins to be degraded contributes to the congestion of ClpXP but also the LacI repression helps to de-congest by reducing the tagged GFP expression. Thus, as Ara concentration is increased, we notice a shift towards higher RFP levels along with a dispersal and lower levels of GFP values. This defines the domain (gray window) where the feedback loop required for SOC can effectively operate. A SOC state is obtained in the presence of high Ara concentration around the IPTG values close to the queueing transition. Most cells (Q3 + Q4 ≈ 80%) are emitting in the red channel, but exhibit broad range of green fluorescence levels, since this state is characterized by fluctuations associated to large bursts of GFP expression (the heterogeneous GFP expression is apparent in the yellow cells of (c) and in the histogram of Fig. 4). Flow cytometry analysis of more IPTG-Ara combinations are shown in Section III.D, SM.
production of GFP-lva. The repression loop is consistently removed when the production of GFP-lva (our order parameter) is reduced enough as for desaturate the ClpXP protease that can degrade the repressor again. The addition of Isopropyl β-d-1thiogalactopyranoside (IPTG), as indicated by a negative input to the repression feedback (see Fig. 3b), switches on our SOC circuit and allows to control the level of GFP-lva expression. High levels of IPTG will lead to an overproduction of GFP-lva and the subsequent ClpXP complex congestion, thus reproducing the limit case that would correspond to the standard phase transition of the queueing process [26] (Fig. 1ab).
The SOC motif is contained in a single, high-copy plasmid to ensure a maximal concentration of the vector while maintaining the parity of the two parts of the circuit (Fig. 3a-b). The construct was transformed in XL1-Blue E. coli strain. Further details of the cloning process and sequences can be found in Section S.III, SM. Also, this design allows an easy tuning to obtain parameters involving SOC, in terms of the strength of the promoter (i. e. pBAD) and by adjusting the efficiency of the repressor (in our case, IPTG for LacI, see Section III, SM).
To perform the experiments, a single colony was inoculated in a volume of 4 ml, and grown at 37 • C until the exponential phase was reached with an OD 660 around 0.6. This homogeneous fresh culture was then used to inoculate all the conditions used in the experimental design. Each combination was inoculated with 1 µl of the starter culture to a final volume of 4 ml. Cells were grown for about 10 hours at 37 • C, until reaching an approximate OD 660 of 0.8-1. The output of the each condition was then analysed using both Fluorescence-Activated Cell Sorting (FACS) and fluorescence microscopy (Fig. 3c-d).
The results from the FACS are displayed in Fig. 3c, where a 4 × 4 array of different combinations of IPTG and Arabinose concentrations define our parameter space by means of dot plots. The range of con-centrations shown here are 0 ≤ mM [IPTG] ≤ 1 , 0 ≤ mM [Ara] ≤ 100 (see the specific values in Fig. 3c and in Section III, SM). As described above, these small molecules allow to explore a parameter space where we can move from a decoupling between the two genes to a full-fledged repression feedback required for criticality to occur. The different cell population responses to the tuning of both IPTG and Ara reveal the relative impact of each on the SOC motif. The target for a SOC state implies two requirements: (i) the expression of large enough levels of the control parameter to effectively perform its feedback; and (ii) a GFP expression characterised by bursts but displaying a low average activity. In the non-induced state, without arabinose (left column), increasing levels of IPTG concentration promote a standard transition from the free to the congested phase (sub-and super-critical phases, indicated in the two bottom and at the two top panels as SUBCRIT and SPCRIT, respectively). As IPTG grows, we effectively weaken the strength of the repression loop until a critical point is reached allowing congestion to rise. This is clearly observed from the displacement of the density dot plots from low to high levels of GFP.
As Arabinose concentration increases, we move in the other dimension of our parameter space, where the control molecule gets more common (cells emit in the red channel) but cannot always effectively act as a repressor. This clearly is a time point picture shot of a bacterial population that exhibits the fluctuating GFP levels characteristics of the SOC state. The same SOC behaviour is shown in the fluorescence microscope image of Fig 3c, where some bacteria do not have the ClpXP saturated (and thus do not display fluorescence), many are near the critical state of ClpXP saturation with lower levels of effective LacIlva to repress the GFP, and exhibiting a wide range of GFP-lva concentrations (bacteria in yellow) and few bacteria have enough Laci-lva to degrade the GFP (bright only in red).
The experiments with E. coli confirm SOC fluctua- thus leading to a scaling exponent γ ∼ 3, consistent with the stochastic simulations. The insets display the comparison of the raw histograms of the congested state (green dots (b)) and free-phase (blue dots,(c)), respectively and the SOC state (gray dots) corresponding to the same IPTG conditions. Additional distributions with a more detailed IPTG-Ara combinations are shown in Section III.D, the SM.
tions of GFP-lva, while the reporter of the control element (i. e. RFP-lva) remains basically stable in terms of the concentration levels and their dispersal. The experimental system successfully reproduces another important feature of criticality, namely the presence of a power law in the dynamics of the order parameter σ 1 . From the existence of a transition in the queueing process between the two phases described above, we can conjecture that the SOC motif should easily organize our gene network into this critical boundary provided that the control of the dimer concentration allows the loop to work properly. This is precisely what we found, as shown in Fig. 4a, where we display the same set of Arabinose-IPTG combinations shown in Fig. 3d. Here the statistics of expression are shown in Fig. 4a using cumulative distributions. The histograms of the non-induced E. coli colony with low Arabinose but tuned using IPTG (left panels) reveals a single-peak shape (a flat part followed by rapid decay in the cumulative plot). In general, as we increase the levels of both arabinose and IPTG, the distribution of our order parameter becomes fat-tailed once RFP levels become high, and well-defined power laws can be observed, as the two highlighted in Fig. 4b-c (the insets are linear-log plots to highlight the different behaviour of the two components of the SOC motif). A detailed parameter exploration is provided in Figs. S10-S17, SM.

III. DISCUSSION
Self-organized criticality (SOC) has a seemingly paradoxical nature: it involves steady states that are always on the edge of instability. Are there intracellular processes poised close to critical points? Traffic dynamics in other contexts suggests that optimal flows occur close to criticality along with very broad fluctuations exhibiting optimal flow at critical points [22,23]. Within cells, theoretical work suggests that enzymatic networks might be poised to criticality when the substrate input rate is balanced by the processing capacity of the enzymatic network [35] and that SOC states might pervade optimal growth [37]. Such critical bal-ance would be a source of adaptation. In this paper we have followed a constructive approach by building a new type of network motif implementing the logic of SOC processes on a two-gene network by following the basic design principle of linking order and control parameters [31]. As the activity level (σ 1 ) grows due to an overloaded proteolytic machinery, the competition for the ClpXP pool also increases the levels of the control component σ 2 that can dimerize to perform a negative control on "emission" of σ 1 thus effectively reducing activity.
Using the SOC motif architecture, it was possible to create a separation of time scales driving to highly fluctuating, critical dynamics. This work shows for the first time that this class of "unstable attractor" can be engineered in living cells. Interestingly, this is similar to the behaviour of computer traffic models where packet production is regulated by the amount of actual congestion [38], which leads to fat-tailed distributions of packets. Being at the critical state has important consequences linked with optimality and might be relevant for information-processing tasks. Several authors early suggested that biological computation could occur close to phase transitions [19,39] and given the potential effects of a critical motif on other cellular systems performing given tasks, our results could give support to this conjecture at the cellular level. In this context, Hasty and co-workers have shown that a proper engineering of the proteolytic machinery can be used to achieve relevant functionalities, including tunable post-translational coupling [40] or in vivo drug delivery based on pulses of bacterial lysis against colorectal tumors [41]. Our SOC motif could further enhance some of these applications (wide fluctuations and rapid responses to external signals). An obvious extension of the critical motif could be a multicellular circuit able to trigger population-level avalanches by exploiting the quorum sensing machinery. Similarly, the fat-tailed behaviour could be wired to a diverse range of functionalities, such as search paths with fat-tailed statistics where bursting dynamics have adaptive value [43,44].
Critical states are known to be part of the cognitive equipment of multicellular organisms, from the simple, non-neural placozoans to neural systems and animal collectives [45,46]. The SOC motif might be a very efficient way of generating the highest phenotypic diversity in a microbial population and can be relevant to expand space of synthetic biology computational designs [47] into collective intelligence [48]. A missing point here is the lack of a time dimension that could help confirming our results and further develop a theoretical framework. This can be achieved by constructing a similar SOC motif within an eukaryotic cell, where the time scale of the resulting time series would be smaller than the cell division cycle. Finally, given the analogies between our system and critical traffic in parallel computer networks, an extension of our approach could involve a 3D spatially explicit system and the development of statistical physics models of critical intracellular traffic.

A. Plasmids construction
Plasmid construction and DNA manipulations were performed following standard cloning techniques. The LacI-lva (BBa C0012), RFP-lva (BBa K1399001) and GFP-lva (BBa K082003) genes were amplified from the parts registry collection (2016). The forward primers were synthesized to contain the proper promoter and/or RBS sequences: the pBAD promoter and the RBS30 for the LacI gene, RBS34 for RFP gene and pLacIQ promoter with RBS34 for GFP gene. The PCR products pBad-RBS30-Lacil-va and RBS3-RFP-lva were joined together by assembly PCR, and cloned to pBluescript plasmid in the restriction sites EcoRI and XbaI. The PCR product pLacIQ-RBS34-GFP-lva was cloned to a Bluescript plasmid by SpeI and PstI. The resulting plasmids were joined together by ScaI and the blund ends of Eco53kI and EcoRV. The clonings were realized in the pBluescript II SK(+) plasmid backbone (ColE1 high copy number replication origin). See sequence of primers in the supplementary table I (see also Fig. S9, SM).

B. Strains and growth conditions
Plasmid cloning and evaluation of the circuit behaviour was performed in E. coli XL1-Blue strain. All characterisation experiments were done in lysogeny broth (LB) Lennox media (10 g/L Tryptone 5 g/L Yeast Extract, 5 g/L NaCl ) with a final ampicillin concentration of 125 µg/mL. Single colonies were inoculated in 4 ml and grown at 37 • C with shaking (200 r.p.m.) during 4 hours, to reach an approximate OD 660 of 0.6. One microliter of the culture was re-inoculated in 4ml of fresh media, supplemented with ampicillin, and the corresponding Arabinose and IPTG concentrations. The cultures were grown overnight (10-14 hours) at 37 • C with shaking. Once they were at OD 660 of 0.8-1, were used for fluorescence measures.

C. Imaging of single cell gene expression
The output of the SOC circuit was analyzed after 10h of incubation at 37 • C with different combination of inputs. Samples were diluted in PBS and analyzed using flow cytometry (BD LSRFortessa). A total of 10 4 cells were collected from each sample. Specific emission fluorescence channels for GFP (FITC-H) and RFP (PE-H) were measured. A proper gate to subtract the debris particles was set using forward and side scattering channels. For the FACS graphics, the GFP and RFP fluorescence of cells inside the gate were plotted in adjacent axes. The cumulative distributions depict all bacteria with a FITC-H expression above 10 2.5 . All data was analysed and plotted using FlowJo (v7) software and customized Phyton code. The regression line and slopes of the histograms were calculated using Numpy and ploted with Matplotlib.
For microscopic images, the cells were harvested at the same time than the cytometry analysis and pictures were collected with a inverted microscope Leica DMI6000, using a 40x oil objective. Bright field, red and green fluorescent images were taken, and then merged using ImageJ.

D. Mathematical modelling
The mathematical model used here is a twodimensional system of nonlinear ordinary differential equations describing the coupling between the order (σ 1 ) and the control (σ 2 ) parameters required to obtain criticality: where the following Hill function response [33] is used: for the repression mediated by σ 2 σ 2 dimers. The parameter µ ∈ [0, 1] weights the effect of IPTG on the strength of the negative control. When σ 2 is small (the ClpXP system is working far from congestion) we have a production rate f (σ 2 → 0) ≈ η 1 /θ. The inhibition function has a threshold value θ representing the concentration σ * 2 at which the rate drops to half its maximum value i. e. f (σ * 2 ) = η 1 /2θ. For larger values, it rapidly decays to zero. The saturation function, namely Γ(σ 1 , σ 2 ) = δ c C K + σ 1 + σ 2 , introduces the competition of both proteins for the proteolytic machinery. Here, as well, the limit case when no congestion occurs (due to low concentrations of both σ 1 and σ 2 ) gives a constant removal rate proportional to the concentration of ClpXP units, i. e. Γ(0, 0) = δ c C/K. The expression of σ 1 gives the behaviour of the GFP-lva, whereas σ 2 stands for LacIlva. Thus f (σ 2 ) is the expression for the response of pLac (the promoter controlling the expression of GFP) to the LacI protein. In this function, η 1 is defined as the production rate, θ the promoter sensitivity, and finally, µ weights how effective is the repression of LacI (effectiveness being altered by IPTG: the more IPTG the lower the µ value). The production rate of LacI (σ 2 ) is controlled by the pBaD promoter, which will trigger a heavier production when there is Arabinose in the medium: the more Arabinose, the higher the value of η 2 . Both proteins are diluted and degraded at rates δ 1 and δ 2 respectively. Finally, both proteins are degraded by the ClpXP system, that can be saturated if there are enough proteins to be degraded. For this reason there is a sigmoid function, with degradation rate δ c , C standing for ClpXP concentration and sensitivity K. Notice that both proteins compete for the degradation machinery, thus inhibiting each other (being added in the denominator). Stochastic simulations of the previous deterministic model where also implemented using the Gillespie method [32] (see Section II.B, SM).