## Abstract

Inspired by the Ising model, we introduce a gene regulatory network that induces a phase transition that coordinates robustly the behaviour of cell ensembles. The building blocks of the design are the so-called toggle switch interfaced with two quorum sensing modules, Las and Lux. We show that as a function of the transport rate of signalling molecules across the cell membrane the population undergoes a spontaneous symmetry breaking from cells individually switching their phenotypes to a global collective phenotypic organization. By characterizing the critical behaviour, we reveal some properties, such as phenotypic memory and hypersensitivity, with relevance in the field of synthetic biology. We argue that our results can be extrapolated to other multicellular systems and be a generic framework for collective decision-making processes.

## 1. Introduction

The spontaneous collective order/coordination arising from the interactions between individual constituents, i.e. self-organization, is a driving principle in biology and Nature [1–3]. Recurrent examples at different spatio-temporal scales include the formation of the lipid bilayer, tissue and organ morphogenesis or flocking/herd behaviour. Focusing on self-organizing *critical* phenomena, phase transitions appear to be a driving force in some biological processes at the level of the organization of the cell. Thus, the membrane organization/patterning [4–6] or the segregation of P granules during asymmetric cell division [7–9] have been explained using liquid–liquid and liquid–gas analogies. However, bridging scales has been found to be arduous and, to the best of our knowledge, there are no examples of critical behaviour, à la physics, during the *collective* organization of living matter, i.e. in cell populations. Possible examples include the pancreatic islet excitability [10] or the behaviour of cancerous aggregates [11]. We stress that by ‘critical behaviour’ we mean the proved existence of formal properties such as critical points/exponents and universality classes [12]. Notably, in the context of cell populations these features of critical behaviour would imply the possibility of applying unifying levels of description, and thereby mechanisms. This would be especially relevant in cancer or developmental processes for understanding the behaviour of cell aggregates and fate determination. Moreover, additional characteristics, such as the resilience to noise or, in the case of continuous phase transition, memory effects (hysteresis), would provide interesting advantages in synthetic biology. Consequently, identifying and designing regulatory motifs leading to the collective organization of cell populations via *canonical* phase transitions is of major interest for both fundamental and applied research.

Here, as a proof of concept, we introduce a gene regulatory circuit that induces a continuous phase transition in bacterial cell colonies. We show by stochastic numerical simulations that this motif presents all the features à la physics of critical behaviour, including a spontaneous symmetry breaking that is robust to molecular noise and also shows phenotypic memory effects, i.e. hysteresis. Our design rationale is inspired by the Ising model, the prototypical example of continuous phase transitions in statistical mechanics to explain ferromagnetism [12]. Yet, the applicability of the Ising model to biological phenomena has been proved to be fruitful. Thus, the two exclusive states associated with the magnetic moment in the Ising model, spin up/down, can actually be extrapolated to the possible states of DNA regions occupied/non-occupied by transcription factors [13–15]. This has led to a number of applications in the context of gene regulation to explain the patterning of developing embryos [16], the role played by chromatin in gene regulation [17] or cell pluripotency [18,19]. Also, in the context of neuronal dynamics, mapping the spiking activity of a cell to a spin variable has led to studies for understanding the collective response of neurons [20]. Moreover, generalizations of the Ising model, such as the Potts model, are used to simulate the behaviour of cell aggregates/tissues [21]. In our case, we mimic the spin dynamics by making use of the toggle switch synthetic circuit. The latter implements two coupled negative feedback loops that drive the mutual repression of two genes: *lacI* and *λ-cI*. This circuit illustrates the cellular decision-making process and has been extensively studied [22–28]. The toggle switch exhibits a bistable behaviour over a wide range of parameter values with two distinct, not simultaneous, stable states: one with a ‘high’ expression level of *lacI* and ‘low’ expression level of *λ-cI* and vice versa. As for the switching dynamics between states, in the original experimental implementation of the circuit, these states show robustness to fluctuations and transitions between them are rare [22]. Thus, phenotypic flipping between states is achieved by using transient external signals (inducers) that change the activity of one of the two regulatory proteins [24]. Here, on the contrary, we exploit the noise-induced switching features of the toggle circuit when subjected to large levels of molecular noise, that is, when there is a low copy number of the two regulatory proteins. The latter can be experimentally achieved by different means such as modifying the affinity of the ribosome binding sites and/or modulating the so-called bursting size [29]. In fact, an intriguing question is how cells can robustly coordinate in a population in the presence of noise [30–32].

On the other hand, we couple the genetic dipoles, i.e. the toggle switches, at the multicellular level by means of the bacterial quorum sensing (QS) communication mechanism. The latter relies on diffusive signalling molecules, auto-inducers (AIs), which are transported, either actively or passively, in and out of the cell and drive the expression of target genes in a cell density-dependent manner [33]. When the cell density is ‘low’, AI concentration is kept low in the cells and in the extracellular medium and the positive feedback loop that ultimately produces more AI is not active; in contrast, when the cell density is ‘high’, the concentration of signalling molecules inside the cell increases, the positive feedback loop is active and there is both a reinforcement of the cell communication and the expression of downstream target genes. Many synthetic and natural systems based on this mechanism exhibit coordinated behaviour. Synthetic gene circuits have been synchronized by QS leading to the emergence of collective oscillations [34,35]. Also, by coupling cell density and cell motility using the QS mechanism cells exhibit spatio-temporal patterns [36]. In addition, previous work suggests that QS can coordinate a population of toggle switches [37], and experiments have demonstrated that the QS pathway can be used to select one state of the toggle switch at the population level [24]. Here, we extend these ideas and interface the toggle switch with two modified QS modules. Thus, we propose variants of the Lux and Las circuits where the AI synthase genes are placed under the control of the toggle switch repressors.

Importantly, the proposed circuit design rationale can be extrapolated to eukaryotic cells and other cell–cell communication mechanisms, e.g. ligand–receptor dynamics (see Discussion). Consequently, rather than a particular circuit, we introduce here a novel framework for understanding phase transitions in multicellular environments. Our results show that cells driven by this Ising-like regulatory circuit undergo a robust, noise-resilient, *canonical* continuous phase transition with well-defined critical properties. In addition, the cell colony exhibits phenotypic memory effects and hypersensitivity when external inducers unbalance protein production. Altogether, we introduce here a mechanism for collective cell organization that might also lead to applications in the field of synthetic biology.

The paper is organized as follows. In §2, we introduce the building blocks of the circuit design, their effective mathematical description and the observables for characterizing the critical behaviour. In §3, by using the analogies of our system with the Ising model, we show that the cell population exhibits a *canonical* continuous phase transition, phenotypic memory effects and hypersensitivity to external inducers. Finally, in §4, we elaborate on the limitations, extensions and consequences of our results and methodology, and summarize the main conclusions.

## 2. Modelling

### 2.1. The toggle switch as a genetic dipole

In dimensionless time units (protein lifetime units), the toggle switch dynamics for proteins *U _{i}* (LacI) and

*V*(

_{i}*λ*-cI) at cell

*i*(

*i*= 1, 2, …,

*N*) can be described by the reactions [22,27,38] 2.1where is the gene regulatory function of species

*Z*( being the species concentration);

*α*,

*β*and

*K*stand for the basal rate of expression, the maximal unrepressed regulated production rate and the concentrations of repressor protein leading to half-maximum repression, respectively. By setting the parameter values to

*α*= 0.5 nM,

*β*= 10 nM and

*K*= 0.5 nM, the two possible stable phenotypic states are and the corresponding symmetric situation, . Taking into account that, in

*Escherichia coli*cells , these concentrations correspond to approximately 0.5 and approximately 5 molecules per cell on average. These large intrinsic noise conditions lead to noise-driven transitions between

*s*

_{1}and

*s*

_{2}on a reasonable time scale [25] (figure 1). As for the value of the Hill exponent used in our model, typical values used previously vary between 2 [38] and 4 [27]. Importantly, as long as the cooperativity is larger than 1, a bistable behaviour develops [22] and our results would remain qualitatively the same. The original experimental implementation of the toggle switch shows cooperativity exponents between 2 and 3 [22]. In our model, we have chosen the upper bound within those limits, 3.

That is, the system behaves as a genetic ‘dipole’ switching between phenotypes *s*_{1} and *s*_{2} as time progresses. Bias between states can be achieved by means of external inducers that modify the production rate of one of the proteins [24]. This effect is mathematically described by modulating the basal production rate of *U* that can be either increased or decreased depending on an external signal, *h*: *α _{u}* =

*α*+

*h*, where

*h*≥ −

*α*. Unless otherwise stated, we assume a symmetric switch such that

*h*= 0.

### 2.2. Interfacing the toggle switch with quorum sensing circuits

The translation of an Ising-like interaction mechanism into the logic of genetic regulation would mean that cells are able to ‘copy’ their phenotypic state into other cells. This is actually what the QS mechanism promotes by means of signalling molecules. Thus, we interface the toggle switch with modified Lux and Las circuits that place the AI synthase genes, *luxI* and *lasI*, under the control of the repressor molecules, *U* and *V* (figure 2). Note that our design assumes a null Lux/Las pathways' crosstalk. We acknowledge the fact that recent results have shown that the latter exists [39]. However, experimental implementations of synthetic circuits containing both Lux and Las signalling mechanisms, and relying on the orthogonality of their pathways, show functionality [35]. Consequently, Lux/Las crosstalk is, at most, ‘mild’ and we neglect cross-reactivity in our design (see Discussion).

Note that in this circuit *U* promotes the expression of AI* _{U}* by repressing

*V*and AI

*. By the same token,*

_{V}*V*promotes the expression of AI

*by repressing*

_{V}*U*and AI

*. At the same time, AI*

_{U}*stimulates*

_{U}*U*expression and, consequently,

*V*repression (the opposite is true for AI

*). For example, cells in the high*

_{V}*U*/low

*V*state would produce and release the AI

*signal which, upon reaching other cells, would repress*

_{U}*V*and induce the expression of

*U*, therefore promoting the same state in other cells. This enables a simplified and tractable mathematical description since the system behaves, effectively, as if

*U*and

*V*act both as auto-inducers and as responsible for the switching dynamics. By using this effective approach, we describe mathematically the regulatory motif by the reactions (in units of the degradation rates of

*U*and

*V*) 2.2 2.3where

*D*is the effective transport rate across the cell membrane,

*r*=

*V*

_{cell}/

*V*

_{ext}is the ratio between the cell and the external volume, and

*U*and

_{e}*V*stand for the signalling molecules outside the cytoplasmic region. Note that, just by considering viscosity effects, the diffusion of signalling molecules through the cell membrane is, at least, one order of magnitude slower than the diffusion in the extracellular environment (much slower if transport is due to an active mechanism). Consequently, within the time scale of the transport through the cell membrane the spatial effects can be disregarded in our modelling as we can assume a homogeneous concentration in the extracellular environment. We point out that this approximation is valid in small cell colonies. However, in systems where signalling molecules have to cover long distances, then spatial effects are indeed relevant as revealed by the wave propagation phenomenon in synthetic oscillators synchronized by QS signalling [34].

_{e}In order to simplify, we have considered the degradation rate of the external molecules to be the same as that of the internal molecules. Moreover, in our mathematical description we assume that both QS modules and signalling molecules behave similarly. These simplifications would surely lead to quantitative disagreements when compared with actual experimental implementations of the circuit but, importantly, not to qualitative ones with regard to the underlying mechanism/phenomenology (see Discussion).

### 2.3. Control parameter

The usual way to reveal the existence of a phase transition is by plotting critical observables as a function of a control parameter. In the Ising model, the latter is, usually, the temperature: below a critical temperature there is a spontaneous symmetry breaking leading to an ordered phase where dipoles align parallel. Alternatively, the coupling intensity of spins can also be used as the control parameter: for a given fixed temperature, above a critical value of the coupling the phase transition develops [12]. The corresponding situations in the gene circuit would be either to decrease the intrinsic noise (the ‘temperature’) or to increase the transport rate across the cell membrane (the ‘coupling’). Here, we choose this second option for the sake of characterizing the critical value of the effective transport rate. However, to consider the transport rate as the control parameter leads to some caveats in regards of the independence of the control parameters and the effect of the system size growth.

First, in the Ising model, temperature and coupling are independent parameters but in the gene circuit intrinsic noise and transport rate are not. Transport allows molecules to move in and out the cell and, if varied, modifies the average concentration of proteins and, consequently, the intensity of the intrinsic noise. Yet, if the transport rate is our control parameter then the noise intensity, i.e. the ‘temperature’, must be kept constant. Thus, for being consistent with this scheme we need to modify the rates *α* and *β* when varying the transport rate in our simulations such that the same average concentration for the states *s*_{1} and *s*_{2} is maintained. In order to address this issue, we note that the set of rate equations (2.2) and (2.3) leads to the following deterministic description for the dynamics of the concentration of molecules inside cell *i* and in the extracellular environment:
2.4
2.5

By assuming a concentrated culture condition, that is, *Nr* = 1 (see below), and imposing that in the steady state and (the conditions for the symmetric stable state are the same), we obtain the following dependence of the production rates *α* and *β* as a function of the transport rate for maintaining the same average protein concentration as *D* varies:
2.6and
2.7Consequently and .

Second caveat: in the Ising model, the number of spins is kept constant but in the case of a cellular system, if no dilution protocols are implemented, the number of cells increases and that effectively modifies the transport rate. Note that the ratio of the external volume to the volume occupied by the cell population effectively modifies the transport rate (note that *r* = *V*_{cell}/*V*_{ext} appears in the transport reactions). In the following, when the number of cells in our simulations, *N*, is kept constant, we assume a crowded culture as stated above in order to get a reliable coupling between cells and we fix the ratio *Nr* = *NV*_{cell}/*V*_{ext} = 1. This corresponds to cell concentrations around . Cell densities of that magnitude can be observed in microfluidic traps [40] and can be maintained in time by dilution protocols driven by trap size constraints [41,42]. On the other hand, when the cell population is not kept constant but the colony size grows, the relation *Nr* = *NV*_{cell}/*V*_{ext} = 1 does not hold since the external volume continuously decreases as the volume occupied by the colony, , continually increases. Thus, as the colony grows, *r* increases and, even if *D* is kept constant, as happens in realistic conditions, the transport effectively increases as time progresses.

In our simulations, under growing conditions, individual cells are supposed to increase their volume as: being the volume of cell *i* at time *t*, *τ _{i}* being its cell cycle duration. As for the dynamics of the cell cycle, our

*in silico*experiments assume the cycle duration to be a random variable. Thus, the probability density for a cell to have a cycle of duration

*τ*reads [43]: where is the average cell cycle duration and

*σ*= 10 min is its standard deviation. Note that there is a minimal duration of the cell cycle, . At division time, the cell cycle duration is reset and proteins are distributed binomially between daughter cells.

_{τ}### 2.4. Observables and critical behaviour

As for the observables for characterizing the critical behaviour, by analogy with the Ising model, the local order parameter is defined as the ‘spin’, , where *S _{i}* = + 1 when a cell is in the state

*s*

_{1}and –1 when a cell is in the state

*s*

_{2}. Accordingly, we define the global order parameter,

*M*, as the average ‘magnetization’ [44], , where and stands for the average over time and experiments. A disordered phase corresponds to a situation where cells, at any given time, are equally distributed between

*s*

_{1}and

*s*

_{2}phenotypic states. On the other hand, a fully ordered phase corresponds to a symmetry breaking situation where and all cells are at either the

*s*

_{1}or

*s*

_{2}states. Following the Ising model, we define the ‘susceptibility’ as . The reduced control parameter is defined as ,

*D*

_{c}being the critical effective transport coefficient. Thus, is zero at the critical point and positive/negative in the disordered/ordered phase. Around , i.e. close to the critical point the following scaling relations hold for continuous phase transitions [44]: 2.8and 2.9where and are the rescaled magnetization and susceptibility, respectively. In the thermodynamic limit, the critical exponents,

*β*and

*γ*, are defined as a function of the reduced control parameter as

In practical terms, if the system size is small, the phase transition gets masked and it is difficult to establish precisely the location of the critical point. In these cases, the cumulant intersection method helps to locate the critical point without any previous knowledge of the critical exponents [45]. Thus, the Binder cumulant, does not depend on the system size at the critical point. Moreover, as the system size goes to infinity, the Binder cumulant converges to the following values: and .

## 3. Results

### 3.1. A phase transition as a function of the transport rate

For the sake of characterizing the value of the critical transport rate, we perform stochastic simulations of the rate equations (2.2) and (2.3) using the Gillespie algorithm [46] for different fixed population sizes, from *N* = 50 to *N* = 10^{3} cells, and compute *M* and *χ* as the transport rate, *D*, varies (figure 3; electronic supplementary material). We set the parameter values as stated above and the initial condition corresponds to a situation where half the cells are at state *s*_{1} and the other half at state *s*_{2}. As *D* increases we observe a continuous phase transition from an unordered to an ordered state. As for the former, the genetic dipole in each cell switches as time progresses between *s*_{1} and *s*_{2} and, as a result, *M* is null and no global coordination develops. As for the ordered state, as *D* increases the communication mechanism elicits an emergent behaviour and cells ‘cooperate’ to align their phenotypes, resulting in an ergodicity breaking and a collective organization. The phenotype, either *s*_{1} or *s*_{2}, is selected randomly and depends on the particular realization of the molecular noise. The behaviour of *M* and *χ* reveals a continuous phase transition. The susceptibility presents a peak around the critical effective transport rate and as the population size increases the transition gets more abrupt and the peak diverges. Note that *M* reaches a maximum value that is slightly higher than 1. We attribute this to the fact that the positions of the stable states in stochastic systems differ from the positions of the deterministic stable states, as recently reported [32]. In order to characterize better the phase transition and the values of the critical exponents, we use finite-size scaling theory [44]. We locate the critical point in the thermodynamic limit using the cumulant technique, as shown in figure 3: *D*_{c} = 0.797. Taking into account the value of the degradation rate, corresponds to effective transport rates in the interval (10^{−3}–10^{−2}) min^{−1}. Typical values of the passive transport rate of auto-inducers are around 10 min^{−1} [47], i.e. above *D*_{c}.

Once we have established *D*_{c}, a fit to the finite-size scaling relation reveals the value of the critical exponent *ν* (figure 4): *ν* = 1.7 ± 0.2.

Also, the height of the susceptibility peak depends on the system size as
3.1where *b* is a constant. Fitting *χ*_{max} to equation (3.1), we obtain *γ* = 0.8 ± 0.1 (figure 4). Finally, the magnetization at the critical point as a function of *N* reads
where *c* is a constant. Data regression yields *β* = 0.1 ± 0.4 (figure 4).

In summary, the critical exponents read *ν* = 1.7 ± 0.2, *γ* = 0.8 ± 0.1 and *β* = 0.1 ± 0.4.

In order to test the values of the critical exponents, by following equations (2.8) and (2.9), we collapse *M* and *χ* into single curves using different system sizes *N* (figure 4).

We note that the critical exponents of the Ising model for a dimensionality larger than 3, i.e. as estimated by the mean-field theory, are [12] *ν* = 0.5, *γ* = 1 and *β* = 0.5.

We then conclude that the universality class of the phase transition is different from that of the Ising model. The latter can be explained in terms of the symmetries that are broken when a phase transition develops, i.e. different broken symmetries lead to different critical exponents [12]. In particular, in the Ising model, the symmetry that is broken corresponds to the spin inversion *S _{i}* ↔ −

*S*. In our model, the latter is related to the symmetry that is actually satisfied by equations (2.4) and (2.5): the ‘replacement’ symmetry

_{i}*u*↔

*v*. Note in fact that equations (2.4) and (2.5) do not satisfy the inversion symmetry

*u*↔ −

*u*and

*v*↔ −

*v*. The development of a phase transition in our model ultimately requires the replacement symmetry to be broken and consequently the critical exponents are different.

### 3.2. Collective organization under exponential growth conditions

Once we have shown that the model undergoes a continuous phase transition as a function of the transport rate for a fixed population size, we explore the collective organization of cells under more realistic conditions when all parameters are fixed and cells constantly grow and divide. Thus, by using the same parameter values as above, but starting with a single cell in a volume of 600 µm^{3}, we run stochastic simulations using a fixed transport rate *D* = 2 (above the critical point) and keep constant *α* and *β* during the whole simulation. In this case, as shown in figure 5, the number of cells increases and, upon reaching a particular value, the ‘magnetization’ reveals the appearance of a global coordinated state where cells ‘align’ their phenotypes. The selected collective phenotype depends on the realization of the random processes (figure 5). The critical system size, *N*_{c}, satisfies that . At that point *r* = 1, i.e. the same value as in the simulations with a fixed-size population and the symmetry breaking emerges.

### 3.3. Hypersensitivity and phenotypic memory

In the absence of an external field (i.e. inducer), the model satisfies the symmetry *U* ↔ *V*. Still, due to a spontaneous symmetry breaking phenomenon, one of the phenotypes, either ‘*U*’ or ‘*V*’ expression dominant, is eventually selected and locked at the collective level. The symmetry *U* ↔ *V* can be broken by means of external inducers/repressors that unbalance the protein production of *U* versus *V*. While the latter is not a key for driving a spontaneous symmetry breaking, it is important to study how the system behaves in case *U* and *V* are unbalanced. Importantly, the features of the response of magnetic systems under an external field are notable, hypersensitivity and hysteresis, and constitute the basis of important applications as memory devices. In this section, we investigate the properties of the proposed model under the effect of external inducers (either an activator or a repressor of *U*) and consider the case *α _{u}* =

*α*+

*h*, where

*h*≥ −

*α*. By analogy with the response of the spins to an external magnetic field, we expect the response of the cells to the external inducer to be different whether the transport rate is either below, close to or above the critical point (figure 6).

We run stochastic simulations of a population with *N* = 1000 cells. If *D* < *D*_{c} the protocol for varying the external signal reads as follows. For 0 ≤ *t* < *t*_{1}, the concentration of the external inducer is set to zero. At *t*_{1}, the external inducer/repressor is added and we set *h* = *h*_{1} to a constant value. The response of the cells, , is computed for for the range . As expected, if *D* < *D*_{c} there is only one global phenotypic state for the population for every value of *h*_{1}. For *D* = 0, there is no cell communication and the response is identical to the average of the response of a single cell. At *h*_{1} = 0 nM, the cells randomly jump between the two stable states and the ‘magnetization’ is zero. When *h*_{1} < 0 nM (*h*_{1} > 0 nM), the inducer introduces a bias in the switch and becomes negative (positive). The resulting curve is smooth and monotonous, indicating that the modulation of the inducer yields a smooth modulation of the phenotype. For *D* = 0.7, below the critical point, we observe the same qualitative behaviour as in the case *D* = 0. However, the ‘magnetization’ displays a sigmoidal shape that is more sensitive than the transport-less case. Thus, cell communication enhances the collective response: for *h*_{1} = 0 nM, the slope of the response curve for the transport-less case is
a fourfold increase. For *D* = 1, above the critical point, cells are locked into the same phenotype and present a finite ‘magnetization’ in the absence of inducer as previously discussed (figure 5). In this case, the external inducer is switched on at *t*_{1} = 0, before the cell population chooses one of the two possible symmetry breaking states. Our results show that a bias in the production of *U* as low as 0.01 nM (–0.01 nM), i.e. a 2% perturbation with respect to the basal expression level, is enough to drive the symmetry breaking of the cell population to a positive (negative) state: hypersensitivity. Therefore, as in the case of the Ising model, we observe a discontinuity in the ‘magnetization’ in response to an external signal. Moreover, if we compute the response of the population to an external inducer that varies continuously between *h* = –0.5 nM and *h* = +1 nM and back we observe a hysteresis loop. That is, the collective response of the cell population displays memory effects and is different for the trajectory with increasing *h* from that with decreasing *h*; namely, the cell population ‘remembers' its phenotypic state. The width of the loop, Δ*h*, depends on the rate of variation of *h*, *γ*. In the case shown in figure 6, .

## 4. Discussion

Up to now it was shown to be difficult to span over scales and show that phase transitions can coordinate robustly the behaviour of living matter beyond the single cell scale. Here, we have shown that this mechanism can take place also in multicellular environments. This may be particularly important in the field of developmental biology, where the coordinated behaviour of cells requires a robust collective cell fate establishment. The so-called community effect in *Xenopus* shows that cells in multicellular organisms are able to differentiate together when producing and responding to a signalling factor [48–50]. Yet, regardless of recent advances [51], it is still an open question if this phenomenon can be explained using the critical phenomena formalism. Our work may shed light by identifying a minimal framework for phase transitions in cell groups.

Herein we have used a highly simplified and effective description in order to contain the computational complexity, restrict the number of parameters and make the problem as analytically tractable as possible. More realistic simulations of the circuit would need to include membrane transport rates that, while comparable, should be different for the auto-inducers; the same applies to the degradation rates, the behaviour of promoters and a myriad of other details and processes. However, our intention here is not to provide a precise quantification of a regulatory circuit but to unravel the logic of interactions that can lead to critical phenomena during the organization of living matter. Thus, we argue that regulatory interactions that, on one hand, lead to a dipole-like (i.e. dichotomous) behaviour and that, on the other hand, couple each state with the surrounding cells through an autocatalytic communication mechanism are prone to undergo a phase transition. We stress that the communication mechanism does not necessarily require long-range interactions such as those of QS in prokaryotes or morphogen diffusion in eukaryotes [52]. A phase transition develops in the Ising model as long as the dimensionality is larger than 1 [12]. Consequently, short-range interactions, e.g. ligand–receptor, that couple cells with their nearest neighbours, for example that of Delta-Notch, would also elicit a critical behaviour. The Delta-Notch case is particularly interesting due to its active role in a number of developmental processes patterning primordia and the fact that this interaction may lead to a dipole-like response such as that appearing in the lateral inhibition phenomenon [53]. Moreover, the critical exponents of phase transitions are very sensitive to the dimensionality of the system. Thus, the comparison between modelling and experimental realization of a regulatory circuit leading to a phase transition would in principle allow one to determine the effective range of interactions. In this regard, in the case of the circuit proposed herein, the QS mechanism shows long-range interaction due to fast diffusion of auto-inducers and hence high dimensionality. Here, as is usual when describing diffusive signals, we have assumed global coupling, i.e. ‘infinite’ dimensionality (mean-field approach). While it is plausible that the QS communication does not span the whole-cell colony we do not expect major discrepancies in the universality class since mean-field exponents apply as long as the dimensionality is larger than 3 [12]. Nonetheless, in other systems driven by short-range communication the estimation of the universality class would in principle allow one to establish the interaction range and clarify how many cells are effectively coupled. Even so, we acknowledge the fact that characterizing the universality class involves the controlled variation of parameters and that can be a cumbersome experiment to implement.

As recently reviewed [54], the QS mechanism provides huge possibilities for synthetic biology. As a contribution to that field, our regulatory design shows how cell communication can be used to build a robust, yet tunable, phenotypic memory at the level of a bacterial population. Thus, the response to external signals where the phenotype lags behind the actual value of inducers/repressors could be used in biosensors to record, for example, changes in the environment. In any case, we stress that a precise design of synthetic circuits using several QS pathways requires crosstalk to be minimized between signalling molecules. Here we have neglected the reported Lux/Las crosstalk [39] motivated by previous experiments that revealed a nugatory effect [35]. However, it is possible to include additional regulatory controls, e.g. an AND gate, to minimize cross-reactivity. Moreover, there are QS pathways that rely on different signals and show orthogonal signalling, such as that in *Staphylococcus aureus* [54,55]. Altogether, other alternatives to our design are possible. Yet, we believe that using here QS synthetic biology standards, as Lux and Las, and presenting a simpler, yet functional, circuit is more valuable and enlightening.

Finally, one of the analogies with the Ising model that we have exploited is the role of intrinsic fluctuations as an effective temperature. Such parallelism was introduced in the past to formalize the noisy dynamics of gene networks [56]. In the Ising model, temperature fluctuations allow the system to explore the different configuration states. Likewise, noise due to low copy numbers is here a driving force for cells to explore the phenotypic variability at the collective level. Note that if fluctuations were negligible then the cell population lacks ergodicity and the phase transition would not be possible. In that sense, the intrinsic noise plays a constructive role and the mechanism can be considered as a revisited example of ‘order through fluctuations' in the field of self-organization [57].

In summary, we have introduced a novel framework for understanding how phase transitions can be implemented in multicellular environments and drive a collective decision process. This framework goes beyond the particular regulatory example we have proposed or the cell types since phase transitions do not rely on specific details but on the components that define the core mechanism. Besides the importance of the latter, we also expect our work to be relevant for designing synthetic circuits and biosensors driven by two QS signals.

## Competing interests

We declare we have no competing interests.

## Funding

This work was partially supported by the Dolores T. and William E. Schiesser Faculty Fellowship (J.B.). M.W. was partially supported by an FPU fellowship (Ministry of Economy of Spain) for this research.

## Acknowledgements

We thank Dr Jim Gunton for his helpful comments.

- Received December 21, 2015.
- Accepted May 19, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.