ARTICLE https://doi.org/10.1038/s41467-021-24695-4 OPEN Engineering self-organized criticality in living cells Blai Vidiella 1,2, Antoni Guillamon 3,4,5, Josep Sardanyés5, Victor Maull1,2, Jordi Pla 1,2, Nuria Conde 1,2 ✉ & Ricard Solé 1,2,6 ✉ Complex dynamical fluctuations, from intracellular noise, brain dynamics or computer traffic display bursting dynamics consistent with a critical state between order and disorder. Living 1234567890():,; close to the critical point has adaptive advantages and it has been conjectured that evolution could select these critical states. Is this the case of living cells? A system can poise itself close to the critical point by means of the so-called self-organized criticality (SOC). In this paper we present an engineered gene network displaying SOC behaviour. This is achieved by exploiting the saturation of the proteolytic degradation machinery in E. coli cells by means of a negative feedback loop that reduces congestion. Our critical motif is built from a two-gene circuit, where SOC can be successfully implemented. The potential implications for both cellular dynamics and behaviour are discussed. 1 ICREA-Complex Systems Lab, Universitat Pompeu Fabra, Barcelona, Spain. 2 Institut de Biologia Evolutiva (CSIC-UPF), Barcelona, Spain. 3 Departament de Matemàtiques, EPSEB, Universitat Politècnica de Catalunya, Barcelona, Spain. 4 Institut de Matemàtiques de la UPC-BarcelonaTech (IMTech), Universitat Politècnica de Catalunya, Barcelona, Spain. 5 Centre de Recerca Matemàtica, Edifici C, Campus de Bellaterra, Barcelona, Spain. 6 Santa Fe Institute, Santa Fe NM, USA. ✉email: nuria.conde@upf.edu; ricard.sole@upf.edu NATURE COMMUNICATIONS | (2021)12:4415 | https://doi.org/10.1038/s41467-021-24695-4 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-24695-4 I n order to adapt to environmental challenges, biological sys- feedback between control and order parameters is sketched. tems exhibit a diverse array of response mechanisms grounded However, the concept and its implementation have been con- in sensors and actuators as well as in information-processing troversial and even sandpiles have been found to achieve criti- units. Adaptive responses require dynamical features that com- cality only under very slow driving and when some microscopic bine low energetic costs along with fast changes to efficiently properties are properly tuned8,34–36. In this paper, we introduce respond to environmental changes. Flocks of birds and fish these minimal conditions of SOC dynamics in living cells by schools widely fluctuate in time but rapidly reorganize when a engineering the interaction between order and control parameters perturbation (such as the presence of a predator) occurs. Within in a simple two-gene network design. Given the consensus that cells, the noise was early identified as playing multiple roles the presence of this feedback loop is a pre-condition for SOC, the affecting cell fate, population heterogeneity, signal amplification approach taken here requires locating the SOC motif in the right or response to stress1–3. Noise is both an inevitable outcome of parameter space (not only a given point) where the scale-free stochastic molecular interactions and an essential ingredient in behaviour will be the robust outcome. As shown below, the SOC decision making4. motif effectively allows driving the system into a bursting dyna- It has been shown that many complex systems seem to be mical state where gene expression levels (our queue size dis- poised close to so-called critical points separating ordered from tributions) follow power laws. disordered states5–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 Results dynamics with scale-invariant properties emerge9,10. If s defines Two-gene SOC motif model: deterministic and stochastic the total activity in one given event, such as number of firing dynamics. The importance of the queueing dynamics in the neurons11–14, gene expression15–18, number of active ants in a enzymatic processing is illustrated by the E. coli stress response to colony19,20, critical epidemic bursts21 or the size of traffic starvation, which is triggered by an excess of mistranslated pro- jams22–24, the resulting distribution P(s) is a fat-tailed one, fol- teins. Stress can cause a significant increase in the concentration lowing a power-law of the form P(s) ~ s−γ, with a scaling expo- of aberrant proteins, which must be degraded. When such an nent γ usually located within the interval 2 ≤ γ ≤ 38,25. overload occurs, the concentration of the sigma factor (the master Critical points can be reached by fine-tuning a given control stress regulator) builds up, eventually triggering the stress (or bifurcation) parameter. This parameter (such as density of response37. Recent theoretical work also suggests that queueing particles, temperature or reaction rate) directly influences the could be adaptive in parallel enzymatic networks when the input system’s state, as described by the order parameter S (system’s flux of substrates is balanced by the maximum processing capa- activity, for example). One way to criticality based on tuning key city of the network38. Here we go a step further and show how a parameters is well illustrated by enzymatic queueing processes in simple SOC circuit can be actually engineered in vivo. Ref. 26. These authors used the framework of queueing theory Our goal in this work is to define the basic design principle to (QTH) to study the dynamics of different proteins (the customers build a genetic sandpile system that captures the feedback in QTH) that are processed by a downstream set of enzymes that structure shown in Fig. 1e and generates bursting, fat-tailed play the role of servers. They consider the native E. coli protease avalanches of activity. First of all, consider the simple, two-gene complexes (ClpXP) which are a limited resource that can only network circuit shown in Fig. 1e. Two proteins σ1 and σ2 resulting process (degrade) a limited number of incoming proteins. In from their expression will be used as the building blocks for the Fig. 1a we provide a basic diagram considering a protein σ being order and control parameters, respectively, thus implementing a expressed at some given rate η (our control parameter). If the rate SOC feedback loop (Fig. 1f). The aim is to exploit the topology of of protein production is low (queues are short), degradation is the gene-gene interaction in such a way that the system can detect efficient since the proteases can process all incoming σ units (free the degree of congestion of the ClpXP system by using σ2 as a phase). If production is too high, a long queue of molecules sensor of the σ1 levels. Our control protein σ2 can form dimers, i. waiting to be processed will be present (congested phase). The e. σ2 + σ2 → σ2σ2 and these dimeric forms act as inhibitors (see two regimes are separated by a narrow parameter domain Methods, Eqs. (1–3)). If congestion occurs, the abundance of σ2 (Fig. 1b; see also Supplementary Fig. 2) where an optimal balance increases and its negative feedback effects also do so. A standard is reached, along with broad fluctuations in concentrations26. Hill function will be used to model the dimers as transcriptional However, they do not follow power laws27 but exponential-tailed repressors. The construction of our circuit involves two steps: (a) forms PðsÞ  expðs=sc Þ, see Supplementary Fig. 2b–f, with sc engineering the critical motif and (b) adjusting the protein rapidly increasing as we approach criticality28. Here scaling is production levels. This might seem contradictory with the self- found to occur instead in the distribution of latencies, i. e. the organized description of SOC, but an example of why this is time required from the production to the final processing of each required is given by rice piles34 which exhibit SOC for some given particle23. grain aspect ratios. The critical point is a rather unique one. How can these sys- Along with the topology of the SOC motif, separation of scales tems poise themselves into critical states without fine-tuning? An is known to be a characteristic of SOC33. While the control alternative mechanism to reach criticality is provided by self- parameter has slow dynamics (the angle of the sandpile) the organized criticality (SOC)29–31. In this case, control and order system’s response (the avalanche time scale) is fast. In our model, parameters interact in such a way that the system spontaneously two additional parameters are used to favour the presence of self-organizes into a critical state32,33. The canonical example of criticality. These are the promoter efficiency for σ2, labelled η2; SOC is the critical sandpile (Fig. 1c). By slowly adding grains of and an extra inhibition acting on the repressor σ2σ2 indicated as μ sand to the pile (at a rate η), its slope θ increases. In the begin- in Fig. 1e. We need to remember that degradation (and other ning, only a few grains will fall down but the number s of grains dissipative events) affects σ2 and thus a minimal concentration of in an avalanche rapidly grows as the angle of repose θc is this sensor is needed in order to effectively detect congested approached. Once we have θ = θc, the interaction between θ and states. On the other hand, in order to experimentally validate our sand avalanches (the order parameter) will keep the system at model, we need to tune the strength of the feedback (required to criticality31. This is summarized in Fig. 1d where the nature of the trigger a rapid decay of the intracellular concentration of σ1). By 2 NATURE COMMUNICATIONS | (2021)12:4415 | https://doi.org/10.1038/s41467-021-24695-4 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-24695-4 ARTICLE Fig. 1 Paths to intracellular criticality. Tunable critical dynamics can be found in simple genetic circuits (a) where a given gene (here coding for GFP-Iva) is constitutively expressed into a protein σ that decays and is also actively degraded by the cell proteolytic machinery (ClpXP). By tuning expression rate η (d), a critical rate ηc is found to separate a phase of efficient degradation from another 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 (light green; ClpXP fails to degrade all the incoming proteins). An alternative, non-externally tuned path is self-organised criticality (SOC), provided by the sandpile (c, adapted from @ricard_sole Oct. 2017). 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 repression feedback is mediated by σ2σ2 (the Lac repressor dimer) with σ1 and σ2 acting as order and control parameters, respectively. Protein models generated using Pymol Software. tuning these two parameters, we include both the non-SOC Larger values for η2 (Fig. 2c) do not exhibit such a time scales design based on queueing as a special case26 (see Section I in the separation. The analytic and numerical investigation of the SM) and a mechanism to achieve the SOC state. eigenvalues of the fixed point to study its stability properties In Fig. 2 we summarize the behaviour of the two-gene system reveal the presence of a maximum in the ratio of imaginary to real (Fig. 1e; see Methods and Section II in the SM for mathematical parts (See Supplementary Fig. 7), indicating a remarkable change details) using both deterministic (Section II.A, SM) and stochastic in the vector field of the phase portrait when η2 ~ 10−2 (given our (Section II.B, SM) dynamics. A unique stable equilibrium (σeq = set of fixed parameters indicated in Fig. 2). (σ1,eq, σ2,eq), indicated with a solid black circle in the (σ1, σ2) phase To see how these nonlinear flows behave under the presence of portraits of Fig. 2) is found, with a characteristic structure of the intrinsic noise, a stochastic numerical implementation of the two- orbits in the phase space, as shown in Fig. 2a for the unregulated gene circuit has been carried out using thep Gillespie method39. In ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi domain (here η2 = 10−3). As the expression rate η2 of the control Fig. 2d the coefficient of variation CV ¼ hσ 1 i  hσ 1 i2 =hσ 1 i of 2 parameter increases close to η2 ~ 10−2, it is easy to see the the generated time series is displayed against η2 for three values of presence of slow-fast dynamics in the distinct structure of the μ. This coefficient provides a statistical estimate of the variance of vector fields consistently with the SOC requirement of time scale the fluctuations and a well-defined maximum is observed when separation. Here the critical motif allows for large fluctuations in η2 ~ 10−2. In Fig. 2e we have overlapped several stochastic σ1 to occur (Fig. 2b) as shown by the compression of the realisations with the vector field close to the maximum CV (for trajectories in the phase portrait close to the fixed point, to be η2 = 10−2). The density plot reveals that the stochastic system compared with the more homogeneous flow displayed in Fig. 2a visits very frequently the fixed point σeq (orange-red colours in for η = 10−3. The analysis of this system shows that, once close to Fig. 2e), but also wanders far away in the lower part of the phase the equilibrium point, small changes in the control σ2 trigger portrait, where the vector field is faster and pushes the stochastic marked population spikes in σ1 (see Supplementary figure 4). paths far away from the equilibrium point, then returning it back NATURE COMMUNICATIONS | (2021)12:4415 | https://doi.org/10.1038/s41467-021-24695-4 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-24695-4 Fig. 2 Nonlinear dynamics of the two-gene critical motif. (a–c) Orbits for Eqs. (1–2) in the (σ1, σ2) space with η2 = 10−3 (a), η2 = 10−2 (b), and η2 = 0.056 (c), setting η1 = 10−2, δ1,2 = 5 × 10−2, δc C = 10−2, K = θ = 10−3. The nullclines are plotted in red (dσ1/dt = 0) and blue (dσ2/dt = 0). The colour of the arrows of the vector field corresponds to their module (blue: small; red: large). The background colour shows the arrival times to the attractor (shorter times in yellow; longer times in violet). 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 is indicated by dashed lines. Three different values of the coupling parameter μ are used to show the robust nature of the maximum of CV, where the SOC motif has been tuned to generate fat-tailed behaviour, as shown in panel (e), where the hot map is overlapped to the phase space, showing a larger density close to the deterministic attractor as well as the fat-tailed scaling behaviour (f) with P > ðσ 1 Þ  σ 2 1 which gives γ ~ 3 for P(σ1). Here five distributions are shown for 5 independent runs along with their average (dark line). (g) Transition between free and congested phases (green curve, in linear-log scale), computed using the same parameter values as in (e) and (f), and located at η1 ≈ 0.02. The solid black line with dots shows how the value of the σ1 inhibition function, here labelled as ^ηðtÞ, is tuned by the system itself driving it close to the transition value (see also Supplementary Figures 9–11). Here P(η1) is the probability of ^ ηðtÞ to take different values of η1. By contrast, the flows and hot maps for the non-SOC circuit close to the queuing theory transition (h) have a Gaussian pattern (i) with exponential tails as shown by the straight lines in the linear-log insets (here η2 = 10−2 and μ = 0). The distributions are depicted with the violin plots in panel (j) for the indicated parameter values. to the deterministic equilibrium. In Fig. 2f the resulting system bounces back and forth between fluid and congested states distribution is displayed. Specifically, if P(σ1) indicates the as it tends to get close to criticality (see also Section II.A.3, SM). probability distribution of σ1 expression levels R(activity), the By exploiting the SOC motif, which allows reducing the rate of σ1 σ cumulative distribution is defined as P > ðσ 1 Þ ¼ 0 1 PðσÞ dPðσÞ production, the Poissonian dynamics of the original, non- and helps smoothing the random noise exhibited by P(σ1). If the regulated queueing dynamics is transformed into a bursting γ signal with fat-tailed statistics. original distribution follows a scaling Pðσ 1 Þ  σ 1 , the cumulative γþ1 one gives P > ðσ 1 Þ  σ 1 . The stochastic model gives a value of γ ~ 3 (Fig. 2f) estimated from the average of five different runs. Engineering a synthetic SOC circuit in E. coli. The theoretical The time series associated with this parameter combination is and computational modelling predicts that the SOC feedback shown in the inset, revealing a characteristic bursting dynamics loop defined above (Fig. 1e–f) will display bursting dynamics with typical of the SOC state, to be compared with the smooth, single- fat-tailed activity distributions associated with the σ1 protein scale behaviour of the unregulated dynamics setting μ = 0 (order parameter). If the expression level (η2) of the control (Fig. 2h–i) (see also Supplementary Figure 2). protein σ2 is large enough, its concentration will act as a con- The distinct nature of the SOC motif action is captured by gestion sensor and will repress σ1, following the SOC motif looking at the dynamics of the effective driving term (Eq. (3) in design. Otherwise, the system will lack the feedback loop. Here we Methods) that we label as ^η1 (Eq. (S.17), SM), as summarized in show how can we engineer the genes circuit incorporating some Fig. 2g (see also SM). We can appreciate how this dynamic parameters that allow including the non-SOC phase transition as driving behaves by computing the probability (labeled as Pð^ η1 Þ) of described above including intracellular queueing processes as a finding the system at some given driving value ^ η1 . This special case. Since the time scale of expression changes can not be distribution is peaked close to the free-congested transition (see captured at the single-cell level because we may need long time also Supplementary Figures 10–11 where parameters η1,2 are series and cells replicate in short times, no individual time series tuned to check the robustness of the identified SOC behaviour to can be gathered. Instead, the collective response of the system will these parameters). However, there are no true phases now: the be analyzed to detect the presence of a SOC state. The underlying 4 NATURE COMMUNICATIONS | (2021)12:4415 | https://doi.org/10.1038/s41467-021-24695-4 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-24695-4 ARTICLE assumption is thus that we have a colony-level sampling of the encoded by the green fluorescent protein (GFP), and the control dynamical states and the critical dynamics is assessed by looking parameter acting as the congestion sensor, by the LacI repressor at the distribution of cell states and the resulting aggregated protein, the relative expression of which will be estimated by statistics (assuming ergodicity). means of the expression of the red fluorescent protein (RFP). In The explicit experimental design of our SOC circuit imple- both cases, we use the unstable variants GFP-lva and RFP-lva, mentation is outlined in Fig. 3a–b. The order parameter is respectively. The construct expresses GFP under the pLacI NATURE COMMUNICATIONS | (2021)12:4415 | https://doi.org/10.1038/s41467-021-24695-4 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-24695-4 Fig. 3 Engineered gene circuit implementing the SOC motif in E. coli. Gene constructs considering non-induced (a) and induced (b) states. (c) Overlapped bright field and fluorescence images of bacteria-induced with Ara (100 mM) and IPTG (10 μM). Yellow bacteria express both GFP and RFP. d Flow cytometry dot plots (green 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 noncongested proteolytic machinery phase (i.e. free phase, subcritical) to congested phase (supercritical) depends on the tunable GFP-lva production. As IPTG increases, so do the de-repression of GFP expression, ClpXP is not able to degrade the excess of GFP and most cells exhibit high green emission. 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. Hence, as Ara concentration is increased, a shift towards higher RFP levels along with dispersal and lower levels of GFP values is observed. This defines the parameter space domain (grey window) where the feedback loop required for SOC is at work. For the larger concentrations of IPTG and Ara, the SOC loop remains effective, at least within the limits of the experimentally explored parameter space. 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 a broad range of green fluorescence levels, since this state is characterized by fluctuations associated with large bursts of GFP expression (the heterogeneous GFP expression is apparent in the yellow cells of (c) and in the histogram of Fig. 4). Microscope images of this experiment and FACS of more IPTG-Ara combinations are shown in Section III.D of SM and in the source data file. promoter (η1) while the LacI repressor and RFP reporter protein output of each condition was then analysed using both are under the pBAD promoter, with nonleaky tight regulation and Fluorescence-Activated Cell Sorting (FACS) and fluorescence high-level expression inducible by Arabinose (η2)40–42. All three microscopy (Fig. 3c–d). proteins of the circuit are tagged with lva sequence to be The results from the FACS are displayed in Fig. 3d, where a degradable by the ClpXP proteolytic complex. 4 × 4 array of different combinations of IPTG and Arabinose ClpXP is responsible for degrading proteins carrying the SsrA concentrations define our parameter space by means of dot plots. or YbaQ degron sequences, reducing the half-life of a tagged The range of concentrations shown here is 0 ≤ mM [IPTG] ≤ 1, protein from hours to minutes. In an exponentially growing E. 0 ≤ mM [Ara] ≤ 100. As described above, these small molecules coli culture (Optical Density (OD) OD660 from 0–2), the allow to explore a parameter space where we can move from a endogenous levels of ClpX and ClpP are constant and involve decoupling between the two genes to full-fledged repression around 100 ClpXP molecules, which can degrade at least 105 feedback required for criticality to occur. The different cell molecules of GFPssrA per cell per replication cycle. However, due population responses to the tuning of both IPTG and Ara reveal to the limited number of ClpXP protease complexes, the the relative impact of each on the SOC motif. The target for a degrading capacity of ssrA-tagged proteins can be easily saturated SOC state implies two requirements: (i) the expression of large by the overproduction of a synthetic tagged protein43–45. enough levels of the control parameter to effectively perform its The non-induced circuit depicted in Fig. 3a is thus a particular feedback; and (ii) a GFP expression characterised by bursts but instance of our more complex motif (Fig. 3b) that would displaying a low average activity. In the non-induced state, correspond in our case to the presence of endogenous LacI, with without arabinose (left column), increasing levels of IPTG the GFP-lva expression being repressed. The presence of concentration promote a standard transition from the free to Arabinose leads to a strong expression of the repressor LacI-lva the congested phase (bottom and at the two top panels, (our control parameter) and the reporter protein RFP-lva. Only respectively). As IPTG grows, we effectively weaken the strength when high levels of the -lva tagged proteins are reached, and the of the repression loop until a critical point is reached allowing degradation machinery is saturated, there is enough LacI-lva to congestion to rise. This is observed from the displacement of the repress the production of GFP-lva. The repression loop is density dot plots from low to high levels of GFP. Notice that here consistently removed when the production of GFP-lva (our order distributions appear peaked, as expected from the single-scale parameter) is reduced enough as for desaturate the ClpXP theoretical prediction for queueing processes. protease that can degrade the repressor again. The addition of In order to generate bursting patterns, Arabinose concentra- Isopropyl β-d-1-thiogalactopyranoside (IPTG), as indicated by a tion needs to be increased. To achieve this, we move in the other negative input to the repression feedback (see Fig. 3b), switches dimension (horizontal direction in Fig. 3d) of our parameter on our SOC circuit and allows to control the level of GFP-lva space, where the control molecule gets more common (cells emit expression. High levels of IPTG will lead to an overproduction of in the red channel) but cannot always effectively act as a GFP-lva and the subsequent ClpXP complex congestion, thus repressor. This clearly is a time point picture shot of a bacterial reproducing the limit case that would correspond to the standard population that exhibits the fluctuating GFP levels characteristics phase transition of the queueing process26 (Fig. 1a–b). of the SOC state. We emphasize that the distributions in the SOC The SOC motif is contained in a single, high-copy plasmid to zone (grey rectangle in Fig. 3d) are well delineated in the ensure a maximal concentration of the vector while maintaining parameter space explored experimentally. Outside this region, the the parity of the two parts of the circuit (Fig. 3a–b). The construct SOC loop for larger IPTG and Ara concentrations seems to be was transformed in the XL1-Blue E. coli strain. Further details of still effective (bursting dynamics is expected to occur here too). the cloning process and sequences can be found in Section III, An example of SOC behavior is shown in the fluorescence SM. Also, this design allows easy tuning to obtain parameters microscope image of Fig. 3c, where some bacteria do not have the involving SOC, in terms of the strength of the promoter (i. e. ClpXP saturated (and thus do not display fluorescence), many are pBAD) and by adjusting the efficiency of the repressor (in our near the critical state of ClpXP saturation with lower levels of case, IPTG for LacI, see Section III, SM). To perform the effective LacI-lva to repress the GFP, and exhibiting a wide range experiments, a single colony was inoculated in a volume of 4 ml, of GFP-lva concentrations (bacteria in yellow) and few bacteria and grown at 37 ∘C until the exponential phase was reached with have enough LacI-lva to degrade the GFP (bright only in red). an OD660 around 0.6. This homogeneous fresh culture was then The experiments with E. coli confirm SOC fluctuations of GFP- used to inoculate all the conditions used in the experimental lva, while the reporter of the control element (i. e. RFP-lva) design. Each combination was inoculated with 1 μl of the starter remains basically stable in terms of the concentration levels and culture to a final volume of 4 ml. Cells were grown for about 10 h their dispersal. The experimental system successfully reproduces at 37 ∘C, until reaching an approximate OD660 of 0.8–1. The another important feature of criticality, namely the presence of a 6 NATURE COMMUNICATIONS | (2021)12:4415 | https://doi.org/10.1038/s41467-021-24695-4 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-24695-4 ARTICLE Fig. 4 Statistics and power-law distributions for SOC from the experiments. a Cumulative (non-normalized) distributions P>(σ1,2) of GFP and RFP fluorescence levels, here plotted using green and red lines, respectively, for the same set of conditions shown in Fig. 3d. The candidate combinations leading to the SOC state (grey panels) are characterized by a broad range of GFP expression revealed by the tail associated with large bursts. In (b) and (c) two cumulative histograms are shown for (10 μM IPTG, 50 mM arabinose (Ara)) and (7.5 μM IPTG, 50 mM Ara), respectively. Both distributions are close to a scaling law P > ðσ 1 Þ  σ 2 1 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 (grey dots) corresponding to the same IPTG conditions. Notice that, unlike the larger images, the x-axes of the insets are linear (so that exponential laws appear as straight lines). Cumulative distributions are shown for fluorescence levels above 102.5 (Q2, Q3, Q4 from Fig. 3). All cumulative distributions have obtained from a constant population size of 105 cells. Additional distributions with more detailed IPTG-Ara combinations are shown in Section III.D in the SM and source data are provided as a Source Data file. power-law in the dynamics of the order parameter σ1. From the behaviour of the two components of the SOC motif). A detailed existence of a transition in the queueing process between the two parameter exploration is provided in Supplementary Fig. 14–19. phases described above, we can conjecture that the SOC motif should easily organize our gene network into this critical Discussion boundary provided that the control of the dimer concentration Self-organized criticality (SOC) has a seemingly paradoxical allows the loop to work properly. This is precisely what we found, nature: it involves steady states that are always on the edge of as shown in Fig. 4a, where we display the same set of Arabinose- instability. Are there intracellular processes poised close to critical IPTG combinations shown in Fig. 3d. Here the statistics of points? Traffic dynamics in other contexts suggests that optimal expression are shown in Fig. 4a using cumulative distributions. flows occur close to criticality along with very broad The histograms of the non-induced E. coli colony with low fluctuations22,23. Within cells, theoretical work suggests that Arabinose but tuned using IPTG (left panels) reveals a single- enzymatic networks might be poised to criticality when the peak shape (a flat part followed by a rapid decay in the cumulative substrate input rate is balanced by the processing capacity of the plot). In general, as we increase the levels of both arabinose and enzymatic network38 and that SOC states might pervade optimal IPTG, the distribution of our order parameter becomes fat-tailed growth46,47. Such critical balance would be a source of adaptation. once RFP levels become high, and well-defined power laws can be In this paper, we have followed a constructive approach by observed for a wide range of parameters, as we would expect from building a type of network motif implementing the logic of SOC our SOC motif. By contrast with the standard critical points, processes on a two-gene network by following the basic design where marked differences are apparent in the distributions above principle of linking order and control parameters33. As the and below criticality, all the statistical patterns within our domain activity level (σ1) grows due to overloaded proteolytic machinery, of parameters where the SOC motif is effective display scaling the competition for the ClpXP pool also increases the levels of the laws. The two SOC states highlighted in Fig. 3d are shown in control component σ2 that can dimerize to perform negative Fig. 4b–c (the insets are linear-log plots to highlight the different feedback on the emission of σ1, thus effectively reducing activity. NATURE COMMUNICATIONS | (2021)12:4415 | https://doi.org/10.1038/s41467-021-24695-4 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-24695-4 Using the SOC motif architecture, it was possible to create a particles was set using forward and side scattering channels (see Supplementary Fig. 13). separation of time scales driving to highly fluctuating, critical 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 dynamics. This work shows that this class of unstable attractors expression above 102.5. All data were analysed and plotted using FlowJo (v7) software can be engineered in living cells. Being at the critical state has and customized Phyton code. The regression line and slopes of the histograms were important consequences linked with optimality and might be calculated and ploted using Numpy 1.20.3, Matplotlib 3.4.2, Scipy 1.6.3, FlowCyto- relevant for information-processing tasks. Several authors early metryTools 0.5.1, Pandas 1.2.4, DiffEqJump v6.14.1 and Catalyst v6.12.1 For micro- suggested that biological computation could occur close to phase scopic images, the cells were harvested at the same time than the cytometry analysis and pictures were collected with an inverted microscope Leica DMI6000 (Leica LASX transitions19,48 and given the potential effects of a critical motif v3.3 software), using a 40x oil objective. Bright field, red and green fluorescent images on other cellular systems performing given tasks, our results were taken and then merged using ImageJ 1.8.0_172. could give support to this conjecture at the cellular level. In this context, Hasty and co-workers have shown that proper engi- Mathematical modelling. The mathematical model used here is a two- neering of the proteolytic machinery can be used to achieve dimensional system of nonlinear ordinary differential equations describing the relevant functionalities, including tunable post-translational coupling between the order (σ1) and the control (σ2) parameters required to obtain coupling49 or in vivo drug delivery based on pulses of bacterial criticality: lysis against colorectal tumours50,51. Our SOC motif could further dσ 1 ¼ f ðσ 2 Þ  δ1 σ 1  σ 1 Γðσ 1 ; σ 2 Þ; ð1Þ enhance some of these applications (wide fluctuations and rapid dt responses to external signals). An obvious extension of the critical dσ 2 motif could be a multicellular circuit able to trigger population- ¼ η2  δ 2 σ 2  σ 2 Γðσ 1 ; σ 2 Þ; ð2Þ dt level avalanches by exploiting the quorum sensing machinery. where the following Hill function response58 is used: Similarly, the fat-tailed behaviour could be wired to a diverse η1 range of functionalities, such as search paths with fat-tailed sta- f ðσ 2 Þ ¼ ; ð3Þ θ þ μ2 σ 22 tistics where bursting dynamics have adaptive value52,53. for the repression mediated by σ2σ2 dimers. The parameter μ ∈ [0, 1] weights the Critical states are known to be part of the cognitive equipment effect of IPTG on the strength of the negative control. When σ2 is small (the ClpXP of multicellular organisms, from the simple, non-neural placozoans system is working far from congestion) we have a production rate f (σ2 → 0) ≈ η1/θ. to neural systems and animal collectives54,55. The SOC motif might The inhibition function has a threshold value θ representing the concentration σ 2 be an efficient way of generating phenotypic diversity in a 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 microbial population and can be relevant to expand the space of synthetic biology computational designs56 into collective δc C Γðσ 1 ; σ 2 Þ ¼ ; ð4Þ intelligence57. A missing point here is the lack of a time dimension K þ σ1 þ σ2 that could help confirming our results and further develop a the- introduces the competition of both proteins for the proteolytic machinery. Here, as oretical framework. This can be achieved by constructing a similar 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 SOC motif within a eukaryotic cell, where the time scale of the ClpXP units, i. e. Γ(0, 0) = δcC/K. The expression of σ1 gives the behaviour of the resulting time series would be smaller than the cell division cycle. GFP-lva, whereas σ2 stands for LacI-lva. Thus f(σ2) is the expression for the Finally, given the analogies between our system and critical traffic response of pLac (the promoter controlling the expression of GFP) to the LacI in parallel computer networks, an extension of our approach could protein. In this function, η1 is defined as the production rate, θ the promoter sensitivity, and finally, μ weights how effective is the repression of LacI (effec- involve a 3D spatially explicit system and the development of tiveness being altered by IPTG: the more IPTG the lower the μ value). The pro- statistical physics models of critical intracellular traffic. duction 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 Methods respectively. Finally, both proteins are degraded by the ClpXP system, which can be Plasmids construction. Plasmid construction and DNA manipulations were saturated if there are enough proteins to be degraded. For this reason, there is a performed following standard cloning techniques. The LacI-lva (BBa_ C0012), sigmoid function, with degradation rate δc, C standing for ClpXP concentration RFP-lva (BBa_ K1399001) and GFP-lva (BBa_ K082003) genes were amplified and sensitivity K. Notice that both proteins compete for the degradation from the parts registry collection (2016). The forward primers were synthesized to machinery, thus inhibiting each other (being added in the denominator). Stochastic contain the proper promoter and/or RBS sequences: the pBAD promoter and the simulations of the previous deterministic model were also implemented using the RBS30 for the LacI gene, RBS34 for RFP gene and pLacIQ promoter with RBS34 for Gillespie method39 (see Section II.B, SM). GFP gene. The PCR products pBad-RBS30-Lacil-lva and RBS3-RFP-lva were joined together by assembly PCR, and cloned to pBluescript plasmid in the restriction Reporting summary. Further information on research design is available in the Nature sites EcoRI and XbaI. The PCR product pLacIQ-RBS34-GFP-lva was cloned to a Research Reporting Summary linked to this article. 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 Data availability origin). The sequence of primers is shown in the supplementary information (see All relevant data are available from the authors. Nonetheless, the source data is available Section III.A and Figure 12 of Supplementary Materials). in the open source repository OFS https://osf.io/h5cew/?view_only=33c12a4780954fcd bd9f1c0986adfdc6 Source data are provided with this paper. Strains and growth conditions. Plasmid cloning and evaluation of the circuit behaviour was performed in the E. coli XL1-Blue strain. All characterisation Code availability experiments were done in lysogeny broth (LB) Lennox media (10 g/L Tryptone 5 g/ The Code needed to do the analysis use the functions precoded in the indicated packages. L Yeast Extract, 5 g/L NaCl) with a final ampicillin concentration of 125 μg/mL. However, software and custom code used are available upon request. Single colonies were inoculated in 4 ml and grown at 37 ∘C with shaking (200 r.p. m.) during 4 h, to reach an approximate OD660 of 0.6. One microliter of the culture was re-inoculated in 4ml of fresh media, supplemented with ampicillin, and the Received: 30 November 2020; Accepted: 22 June 2021; corresponding Arabinose and IPTG concentrations. The cultures were grown overnight (10–14 h) at 37 ∘C with shaking. Once they were at OD660 of 0.8–1, were used for fluorescence measures. References Imaging of single-cell gene expression. The output of the SOC circuit was ana- 1. Rao, C. V., Wolf, D. M. & Arkin, A. P. Control, exploitation and tolerance of lyzed after 10 h of incubation at 37 ∘C with different combinations of inputs. Samples intracellular noise. Nature 420, 231–237 (2002). were diluted in PBS and analyzed using flow cytometry (BD LSRFortessa; lasers:405- 2. Raj, A. & Van Oudenaarden, A. Nature, nurture, or chance: stochastic gene 488-561-633, detectors: FSC\SSC + 14PTM, Facs Diva software 9). A total of 104 cells expression and its consequences. Cell 135, 216–226 (2008). were collected from each sample. Specific emission fluorescence channels for GFP 3. Eldar, A. & Elowitz, M. B. Functional roles for noise in genetic circuits. Nature (FITC-H) and RFP (PE-H) were measured. A proper gate to subtract the debris 467, 167–173 (2010). 8 NATURE COMMUNICATIONS | (2021)12:4415 | https://doi.org/10.1038/s41467-021-24695-4 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-24695-4 ARTICLE 4. Balázsi, G., van Oudenaarden, A. & Collins, J. J. Cellular decision making and mixed populations. Proc. Natl Acad. Sci. USA 94, 8168–8172 biological noise: from microbes to mammals. Cell 144, 910–925 (2011). (1997). 5. Solé, R., Manrubia, S. C., Benton, M., Kauffman, S. & Bak, P. Criticality and 42. Megerle, J. A., Fritz, G., Gerland, U., Jung, K. & Rädler, J. O. Timing and scaling in evolutionary biology. Trends Ecol. Evol. 14, 156 (1999). dynamics of single cell gene expression in the arabinose utilization system. 6. Mora, T. & Bialek, W. Are biological systems posied at criticality? J. Stat. Biophys. J. 95, 2103–2115 (2008). Mech. 144, 268–302 (2011). 43. Gottesman, S., Roche, E., Zhou, Y. & Sauer, R. T. The ClpXP and ClpAP 7. Solé, R. Phase Transitions. (Princeton U. Press, 2011). proteases degrade proteins with carboxy-terminal peptide tails added by the 8. Muñoz, M. A. Colloquium: criticality and dynamical scaling in living systems. SsrA-tagging system. Genes Dev. 12, 1338–1347 (1998). Rev. Mod. Phys. 90, 031001 (2018). 44. Flynn, J. M., Neher, S. B., Kim, Y. I., Sauer, R. T. & Baker, T. A. Proteomic 9. Stanley, H. E. Phase Transitions and Critical Phenomena. (Clarendon Press, 1971). discovery of cellular substrates of the ClpXP protease reveals five classes of 10. Goldenfeld, N., Lectures on phase transitions and the renormalization group. ClpX-recognition signals. Mol. Cell 11, 671–683 (2003). Addison-Wesley. Reading, MA (2018). 45. Farrell, C. M., Grossman, A. D. & Sauer, R. T. Cytoplasmic degradation of 11. Corral, A., Pérez, C. J., Diaz-Guilera, A. & Arenas, A. Self-organized criticality ssrA-tagged proteins. Mol. Microbiol. 57, 1750–1761 (2005). and synchronization in a lattice model of integrate-and-fire oscillators. Phys. 46. Furusawa, C. & Kaneko, K. Adaptation to optimal cell growth through self- Rev. Lett. 74, 118 (1995). organized criticality. Phys. Rev. Lett. 108, 208103 (2012). 12. Chialvo, D. R. Emergent complex neural dynamics. Nat. Phys. 6, 744–750 (2010). 47. Valverde, S. & Solé, R. Self-organized critical traffic in parallel computer 13. Plenz, D., Niebur, E. & Schuster, H. G. Criticality in neural systems. (Wiley- networks. Phys. A 312, 636–648 (2002). VCH, Weinheim, Germany, 2014). 48. Langton, C. G. Computation at the edge of chaos: phase transitions and 14. Hesse, J. & Gross, T. Self-organized criticality as a fundamental property of emergent computation. Phys. D. 42, 12–37 (1990). neural systems. Front. Neurosci. 8, 166 (2014). 49. Prindle, A., Selimkhanov, J., Li, H., Razinkov, I., Tsimring, L. S. & Hasty, J. 15. Serra, R., Villani, M., Graudenzi, A. & Kauffman, S. Why a simple model of Rapid and tunable post-translational coupling of genetic circuits. Nature 508, genetic regulatory networks describes the distribution of avalanches in gene 387–391 (2014). expression data. J. Theor. Biol. 246, 449 (2007). 50. Din, M. O., Danino, T. & Prindle, A. et al. Synchronized cycles of bacterial 16. Balleza, E. Critical dynamics in genetic regulatory networks: examples from lysis for in vivo delivery. Nature 536, 81–85 (2016). four kingdoms. PLoS One 3, e2456 (2008). 51. Faisal, A. A., Selen, L. P. & Wolpert, D. M. Noise in the nervous system. Nat. 17. Valverde, S., Ohse, S., Turalska, M., West, B. J. & Garcia-Ojalvo, J. Structural Rev. Neurosci. 9, 292–303 (2008). determinants of criticality in biological networks. Front. Phys. 6, 127 (2015). 52. Bartumeus, F. Lévy processes in animal movement: an evolutionary 18. Daniels, B. C., Kim, H. & Moore, D. et al. Criticality distinguishes the ensemble of hypothesis. Fractals 15, 151–162 (2007). biological regulatory networks. Phys. Rev. Lett. 121, 138102 (2018). 53. Hintsche, M. A polar bundle of flagella can drive bacterial swimming by 19. Solé, R. & Miramontes, O. Information at the edge of chaos in fluid neural pushing, pulling, or coiling around the cell body. Sci. Rep. 7, 1–10 (2017). networks. Phys. D. 80, 171–180 (1995). 54. Bialek, W. et al. Social interactions dominate speed control in poising natural 20. Miramontes, O. Order-disorder transitions in the behaviour of ant societies. flocks near criticality. Proc. Natl Acad. Sci. USA 111, 7212–7217 (2014). Complexity 1, 56–60 (1995). 55. Couzin, I. D. Synchronization: the key to effective communication in animal 21. Stollenwerk, N. & Jansen, V., Population biology and criticality: from critical collectives. Trends Cogn. Sci. 22, 844–846 (2018). birth-death processes to self-organized criticality in mutation pathogen 56. Grozinger, L. Pathways to cellular supremacy in biocomputing. Nat. Comun. systems. World Scientific (2011). 10, 1–11 (2019). 22. Nagel, K. & Schreckenberg, M. A cellular automaton model for freeway traffic. 57. Solé, R. et al. Synthetic collective intelligence. BioSystems 148, 47–61 (2016). J. Phys. I 2, 2221–2229 (1992). 58. Alon, U., An Introduction to Systems Biology: Design Principles of Biological 23. Solé, R. & Valverde, S. Information transfer and phase transitions in a model Circuits. (CRC press, 2019). of internet traffic. Phys. A 289, 595–605 (2001). 24. Guimerà, R., Diaz-Guilera, A., Vega-Redondo, F., Cabrales, A. & Arenas, A. Acknowledgements Optimal network topologies for local search with congestion. Phys. Rev. Lett. The authors thank Jordi Garcia-Ojalvo as well as the members of the Complex Systems 8, 248701 (2002). Lab for fruitful discussions. Special thanks to Arianna Bruguera for her help in the 25. Christensen, K. & Moloney, N.Complexity and Criticality. (Imperial College experiments. RS thanks S. Kauffman, S. Manrubia, Jordi Piñero and the late Per Bak for Press, 2005). many discussions on criticality. This work was supported by the Botin Foundation by 26. Cookson, N. A., Mather, W. H., Danino, T. & Mondragon, O. Queueing up for Banco Santander through its Santander Universities Global Division, the Spanish Min- enzymatic processing: correlated signaling through coupled degradation. Mol. istry of Economy and Competitiveness, grant PID2019-111680GB-I00, a MICIN grant Syst. Biol. 7, 561 (2011). PID2019-111680GB-I00, an AGAUR FI 2018 grant, and the Santa Fe Institute (where the 27. Gross, D., Fundamentals of Queueing Theory. (John Wiley and Sons, 2008). key idea was conceived). J.S. and A.G. have been partially funded by the CERCA Pro- 28. Mukherjee, G. & Manna, S. S. Phase transition in a directed traffic flow gramme of the Generalitat de Catalunya. J.S. acknowledges support from Agencia Estatal network. Phys. Rev. E 71, 066108 (2005). de Investigación grant RTI2018-098322-B-I00 and Ramón y Cajal contract RYC-2017- 29. Bak, P., Tang, C. & Wiesenfeld, K. Self-organized criticality. Phys. Rev. A 38, 22243. A.G. has been funded by the AGAUR grant 2017-SGR-1049 and by the 364–374 (1988). MINECO-FEDER-UE grants PGC-2018-098676-B-100 and RTI2018-093860-B-C21. 30. Bak, P. How Nature Works: The Science of Self-Organized Criticality. (Copernicus, New York, 1996). Author contributions 31. Jensen, H.J., Self-Organized Criticality: Emergent Complex Behaviour in R.S. conceived the initial idea, supervised the project and wrote the main paper in Physical and Biological Systems. (Cambridge U. Press, 1998). collaboration with all the authors. N.C., R.S. and B.V. designed the synthetic cellular 32. Sornette, D., Johansen, A. & Dornic, I. Mapping self-organized criticality onto circuit. R.S., and B.V. built the two-gene mathematical model. A.G., B.V., R.S., and J.S. criticality. J. Phys. I 5, 325–335 (1995). performed the mathematical calculations. B.V., A.G., J.S., V.M. and J.P. carried out the 33. Bonachela, J.A. & Muñoz, M.A., Self-organization without conservation: true simulations and numerical analysis. B.V., V.M. and N.C. performed the experiments. or just apparent scale-invariance?. J. Stat. Mech. 2009, P09009 (2009). 34. Frette, V. et al. Avalanche dynamics in a pile of rice. Nature 379, 49–52 (1996). Competing interests 35. Solé, R. V., Alonso, D. & McKane, A. Self-organized instability in complex The authors declare no competing interests. ecosystems. Philos. Tras, R. Soc. Lond. B 357, 667–681 (2002). 36. Watkins, N. W., Pruessner, G., Chapman, S. C., Crosby, N. B. & Jensen, H. J. 25 years of self-organized criticality: concepts and controversies. Space Sci. Additional information Rev. 198, 3–44 (2016). Supplementary information The online version contains supplementary material 37. Gottesman, S. Proteolysis in bacterial regulatory circuits. Ann. Rev. Cell Dev. available at https://doi.org/10.1038/s41467-021-24695-4. Biol. 19, 565–587 (2003). 38. Steiner, P. J., Williams, R. J., Hasty, J. & Tsimring, L. S. Criticality and Correspondence and requests for materials should be addressed to N.C. or R.Sé. adaptivity in enzymatic networks. Biophys. J. 111, 1078–1087 (2016). 39. Gillespie, D. T. A general method for numerically simulating the stochastic time Peer review information Nature Communications thanks the anonymous reviewer(s) for evolution of coupled chemical reactions. J. Comp. Phys. 22, 403–434 (1976). their contribution to the peer review of this work. Peer reviewer reports are available. 40. Guzman, L.-M., Carson, D. & Beckwith, M. J. Tight regulation, modulation, and high-level expression by vectors containing the arabinose pBAD Reprints and permission information is available at http://www.nature.com/reprints promoter. J. Bacteriol. 177, 4121–4130 (1995). 41. Siegele, D. A. & Hu, J. C. Gene expression from plasmids containing the Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in araBAD promoter at subsaturating inducer concentrations represents published maps and institutional affiliations. NATURE COMMUNICATIONS | (2021)12:4415 | https://doi.org/10.1038/s41467-021-24695-4 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-24695-4 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/. © The Author(s) 2021 10 NATURE COMMUNICATIONS | (2021)12:4415 | https://doi.org/10.1038/s41467-021-24695-4 | www.nature.com/naturecommunications