Entropy-driven cell-decision making predicts fluid-to-solid transition in multicellular systems

Cellular decision making allows cells to assume functionally different phenotypes in response to microenvironmental cues, without genetic change. It is an open question, how individual cell decisions influence the dynamics at the tissue level. Here, we study spatio-temporal pattern formation in a population of cells exhibiting phenotypic plasticity, which is a paradigm of cell decision making. We focus on the migration/resting and the migration/proliferation plasticity which underly the epithelial-mesenchymal transition (EMT) and the go or grow dichotomy. We assume that cells change their phenotype in order to minimize their microenvironmental entropy (LEUP: Least microEnvironmental Uncertainty Principle) and study the impact of the LEUP-driven migration/resting and migration/proliferation plasticity on the corresponding multicellular spatio-temporal dynamics with a stochastic cell-based mathematical model for the spatio-temporal dynamics of the cell phenotypes. In the case of the go or rest plasticity, a corresponding mean-field approximation allows to identify a bistable switching mechanism between a diffusive (fluid) and an epithelial (solid) tissue phase which depends on the sensitivity of the phenotypes to the environment. For the go or grow plasticity, we show the possibility of Turing pattern formation for the"solid"tissue phase and its relation with the parameters of the LEUP-driven cell decisions.


Introduction
Cellular decision making allows cells to respond to microenvironmental cues [1]. For instance, if an external ligand from the microenvironment binds to the cell membrane, it may trigger a chain of biochemical reactions inside the cell [2]. After the microenvironmental information has been processed, the cell may respond by changing its properties implying cell differentiation, proliferation, migration, apoptosis etc. Thus, cellular decision making depends on intrinsic signal transduction pathways, the genetic cell network [3], extrinsic stimuli, and molecular noise [4]. Cell-decision making plays a key role in cell fate determination and the maintenance of phenotypic plasticity. In the former, cells irreversibly acquire new fates by following a hierarchical lineage, where pluripotent stem cells in a proper microenvironment (stem cell niche) differentiate into e. g. bone, muscle, epithelial and further specialised cells. On the other hand, phenotypic plasticity allows a reversible adaptation to different microenvironmental stimuli [5]. While cellular decision making is well-studied for single cells, it is an open question how individual cell decisions influence the dynamics at the tissue level.
In this paper, we analyse the implications of phenotypic plasticity on spatio-temporal pattern formation in the cell population and focus on the migration/resting and migration/proliferation plasticity. The migration/resting type of plasticity is observed during the epithelial-mesenchymal transition (EMT/MET) [6][7][8] while the migration/proliferation plasticity characterizes the go-orgrow (GoG) dichotomy found in different types of cancers, such as glioblastoma or breast cancer, but also during normal development [9][10][11][12]. Epithelial cells are immotile, connected by strong adhesion bonds, exhibit regular shapes and show strong apico-basal cell polarity. Epithelial and mesenchymal phenotypes can be interpreted as representing two migration modes, resting and motile, respectively. On the other hand, GoG refers to the mutual exclusion of cell migration and proliferation, leading to the existence of proliferative/immotile and nonproliferative/mobile cell phenotypes. We assume that a Least microEnvironmental Uncertainty Princple (LEUP) dictates the phenotypic switching in both types of phenotypic plastic behaviours. The LEUP principle has been introduced in [13] and assumes that cells change their phenotypes in order to minimize their local microenvironmental entropy. Formulated in the language of statistical physics and Bayesian statistics, cells try to optimize their phenotypic prior to reduce the joint entropy of their own phenotype and microenvironmental variables. Similar ideas using information theoretic measures have been proposed also in seminal works by W. Bialek [14].
The manuscript is structured as follows: In the next section (2), we define a LEUP-driven cell decision model based on the migration/resting plasticity. In Section 3, we develop a stochastic individual-based model (IBM) for moving and resting cells where cells can change their phenotype according to the LEUP-driven cell decision model and show the formation of aggregates of resting cells in corresponding simulations. In section 4, we derive a macroscopic mean-field approximation of our microscopic model which allows to explain the formation of aggregates in the microscopic model. In section 5, we analyze the stability of the steady states and the pattern formation potential of the macroscopic go-or grow model. Finally, in section 6, we discuss the biological implications of our results in terms of multicellular growth and pattern control.

A minimal model of LEUP-driven cell migration plasticity
Here, we define a mathematical model for LEUP-based switching between moving and resting phenotypes. We assume that cells can sense their microenvironment and can change their own phenotype X i on a domain L ⊂ R 2 accordingly, where X i = 0, 1 corresponds to the resting and migrating state, respectively. We define the microenvironment of cell i as the number N 0 i of cells having phenotype X i = 0 and by the number N 1 i of cells having phenotype X i = 1. By assuming a maximum cell capacity N of the microenvironment , we can define where N φ i are the free spaces/slots. The total number of cells N T is defined as the sum of cells having phenotype (X i = 1) and cells having phenotype (X i = 0). So, Cells will select their phenotype in a Bayesian fashion, i.e.
where P (Y i | X i ) can be interpreted as the probability that the cell perceives all other cells in its surroundings, and P (X i ) is the prior probability distribution of the cell's phenotypes. The former models the microenvironmental sensing the intrinsic programming of cellular phenotype. However, sensing other cells and evaluating P (Y i | X i ) entails an energy cost. It is reasonable to assume that the cell will try to optimize its prior P (X i ) for the sake of energetic frugality.
Here we define Y i as the extrinsic random variable of the i-th cell, where Y i = N X i i | N T , X i = 0, 1 represents the number of cells of phenotype X i given the total local number of cells N T . The conditional probability of having N 1 i number of cells present in the microenvironment follows a binomial distribution (B) (for the derivation details see S.I.): where p 1 is the probability of N 1 i number of cells having phenotype X i = 1 out of N T cells: Using statistical mechanical arguments [13], the problem of finding the probability of the cell's phenotype is equivalent to finding the prior that minimizes the entropy of a cell's phenotype and its surroundings. Let S (X i , Y i ) be the entropy of the cell-surroundings, S (X i ) the internal entropy of the cell, and S (Y i | X i ) the entropy of the information sensed by the cell. Then, the entropies are connected by the relation where δ δP (X i ) is the functional derivative,S (Y i | X i ) is the expected ensemble statistics [13], and λ and β are Lagrange multipliers. Taking into account the relations among entropies, eqn. (6) yields where Z is the normalization factor Please note that the parameter β has a biological interpretation, since it quantifies the intensity with which a cell senses and complies to the microenvironment. To estimate the entropy we use the binomial distribution of the resting and migrating cells. Given N T cells in the environment and assuming the frequency of cells at a moving state p 1 = N 1 i N T , the probability for N 1 i moving cells is and the entropy of the binomial distribution in the last equation is When there are no cells in the microenvironment the entropy S(N 1 i , p 1 ) is 0. To facilitate the evaluation of the entropy distribution of the microenvironment, we have used a Gaussian distribution, which approximates the binomial distribution (for N T > 5), i.e.
for the microenvironment [15]. Now we can evaluate the equilibrium distribution P (X i = 0) as with the cell densities ρ 0,1 i := Fig. (12), we show the dependency of P (X i = 0) on resting density ρ 0 and β). Accordingly, the probability of a i-th cell's phenotype P (X i = 1) reads Please note that, the difference between the microenvironmental entropy of the i-th cell is  The microenvironment is defined by the numbers of two phenotypes only. (A2) We have assumed a binomial distribution for the occurrence of the two phenotypes. (A3) Cells are making decisions at a fast time scale. This justifies to assume an equilibrium distribution for the different phenotypes.
In this section, we have defined a minimal LEUP-driven cell decision making model and derived the corresponding phenotypic steady states. Based on this, we subsequently develop a cell-based model to understand the resulting multicelluar spatiotemporal dynamics.  [9]. The switching between motile (X i = 1) and resting (X i = 0) phenotype is shown, where the transtion probabilities are defined as W 10 and W 01 . The proliferation rate is defined by r.
equations. Langevin equations are well-suited to model cell-cell interactions and cell migration [16,17].
Our Langevin's equation is defined on a domain L ⊂ R 2 with periodic boundary conditions. We define an interaction radius ∈ R around the i-th cell at position x i ∈ L. The expected interaction volume is V ∝ 2 . The time evolution of our model is defined by the following rules: (R1) Cells change their phenotype by sensing their microenvironment within the interaction radius according to LEUP.
(R3) Once cells become migratory they move with a constant speedv. The above are translated to the following Langevin's equations where v i = (cos θ i , sin θ i ) T is the direction of movement of cell i. For the temporal evolution of the probability p i of the motile state, we use the BGK (Bhatnagar-Gross-Krook) operator technique [18]. We assume that the probability p i evolves weakly out of its equilibrium probability p eq i , which is the LEUP steady state probability P (X i = 1) (see eq. (13)). In turn, the parameter τ is the relaxation time towards the corresponding probability distribution p eq i . Here noise is assumed to have a zero-mean, white noise term, which has the statistical properties ξ θ i (t) = 0 and ξ θ where t 1 and t 2 are two time points, D θ ∈ R + is the angular diffusion coefficient, δ(t) is the Dirac delta, and δ ij is the Kronecker delta. Parameterv is the constant speed of every motile cell.
We simulate the Langevin model, on a two-dimensional domain, with a varying interaction radius , mean cell density and sensitivity β. We assume an initial state that is approximately homogeneous in space. For sufficiently large densities and sensitivities, we observe aggregation patterns, as shown in Figure (2). To quantify cell clustering, we calculate the radial distribution function for different sensitivities and different interaction radii as shown in Figures (3a) and (3b). The key observations from the simulation study are: 1. There is a critical threshold for the parameter β > 0, where clustering of resting cells occurs (Fig. 3a). This is rather expected since a sufficient sensing of the microenvironment is required. Fig. 3b. Very high corresponds to a large number of sensed cells. This leads to also equal steady state probabilities, as the entropy difference in the microenvironment, associated with the single cell states becomes negligible. On the other hand, the lower bound of is expected, since enough sampling size of cellular microenvironment is required to induce aggregation patterns.

For an intermediate interaction radius , we observe cell clusters as shown in
3. In Fig. 3d, we show the phase diagram of the system, when parameter β and the average density are varied. In particular, we observe that there is a parametric regime where clustering behavior is emerging. Interestingly, high densities do not always imply patterning and require higher values of β to support cluster emergence (see Fig. 3c). Finally, as expected, low densities also reduce the area of the patterning regime for low β.

Mean-field approximation
In this section, we derive a continuous approximation of the aforementioned discrete model using a mean-field approach. The goal is to shed light of the pattern formation mechanisms, i.e. cell clustering, as observed in the microscopic simulations. The main idea of the mean-field approximation is to replace the description of many-particle interactions by a single particle description based on an average or effective interaction. Thereby, any multi-particle problem can be replaced by an effective description, that can be stated in the form of an ordinary (ODE) or partial differential equation (PDE). In order to proceed, we will first treat the switch dynamics and the migration process separately.

Mean-field description of the phenotypic switch dynamics
Let's consider cell i at position x i with states X i = 0, 1. The cell changes its state in dependence on the local microenvironment according to the LEUP with rates W 01 , W 10 and the probabilities to be in either states are denoted by p 0,1 . We assume that the system is always close to the steady state, so the master equation reads and p 0,1 are given by eq. (7). By rearranging the terms, we obtain for the switching rates For simplicity, we set W 01 = p 1 where the transition probability towards motile phenotypes equals to the moving steady state probability p 1 . This coincides with the "detailed balance condition". The steady state probabilities p 0,1 depend on the number of cells N 0,1 i in the There is an optimal sensing radius for a given mean density, so that interaction radii that are too large or too low do not lead to aggregation, indicated by an almost flat radial distribution, see purple and blue lines. (c) Radial distribution function for different mean densities, at β = 20 and = 4. When the mean density becomes too large, no clustering can be observed. (d) Phase space diagram for cluster pattern formation. Here we define clustering when max r g(r) > 1.09 (see text for explaination). The interaction radius is fixed at = 4. All data points in (a -d) correspond to the mean of 20 independent simulations. respective phenotypes in the microenvironment. Cell i at position x i senses Here, we sum over all cells j, and ξ x j (t) − x i (t), t is a Boolean stochastic variable that serves as the sensing function of a cell at x i . It depends on the distance between cells and time, and ξ = 1 indicates that the cell j at x j is sensed by cell i at position x i . To match the IBM, we assume that with the Heaviside step function Θ(x), so that all cells in a ball of radius around x i are sensed.
To proceed, we apply a mean-field approximation to calculate the expected switching rate . Note that we dropped the dependence on space and time for better readability. Let φ i (X, x, t) denote the probability of finding cell i in a small volume dV around x at time t with the phenotype X = 0, 1. Formally we have where the average is the ensemble average, and the total densities of resting/migrating cells are For the expected density of sensed cells in the microenvironment around x we obtain where V ∝ d , and where we have already used eq. eq. (19). Consequently, we obtain the approximate switching rateW and the switch dynamics is (dropping dependencies on space and time for simplicity) Summing over all cells, we obtain for the total density of resting cells Note that this is a non-local, non-linear set of PDEs, which are difficult to treat analytically. However, we can further simplify our analysis by making another approximation, assuming that the sensing radius is much smaller than the total domain size. Then, we can replace the non-local sensing function ξ x j (t) − x i (t), t by a local delta distribution In this case, the expected number of sensed cells in the microenvironment around x simply becomes the local density of the respective phenotype Finally, the rateW 10 reduces tō and the reversed transition probabilityW 01 is where the volume is defined as V = d .

Mean-field description of the cell migration process
In this section, we derive the macroscopic equation in two spatial dimensions for the motile cell population. As before we can write the corresponding stochastic Langevin's equations as d dt This process can be considered as a special kind of active Brownian motion. In this kind of Langevin's equation, the stochastic force creates variations of orientation. According to [19] we can derive the corresponding Fokker-Planck equation for migrating cells using adiabatic elimination and averaging it over different noise realizations obtaining the following diffusion equation Here we define D 1 as diffusion coefficient which isv 4 2D θ .

Coupling migration and switching dynamics
Combining our results in the previous sections, we can easily formulate a system of PDEs where E(ρ 0 , ρ 1 ) is the phenotypic exchange term and ν is the corresponding timescale ratio of the switching and diffusion processes, i.e. ν = τ S τ D . To ensure the numerical consistency of the above system, we assume that resting cells diffuse in a very slow manner i.e. D 0 D 1 , which results in the following reaction-diffusion system of equations

Spatio-temporal dynamics of the LEUP-driven migration/proliferation plasticity
In this section, we study the mean-field approximation of the aforementioned stochastic plasticity dynamics for different regimes of the sensitivity β and the interaction radius . In particular, we initially study the system dynamics, when cells have a very large interaction radius ( 1) and then a finite one. Finally, we study also the special case β = 0, i.e. cells decide independently of their microenvironment.

Large interaction radius case: existence of a critical sensitivity
Here, we focus on the very large interaction radius system ( 1 =⇒ V 1). Although, this parameter regime is not biologically relevant, it is very instructive since it allows us to derive analytical estimates for our systems dynamics. Here, we generalize the macroscopic system by adding proliferation dynamics (logistic growth). Our phenotypic switch dynamics are recapitulated by setting the proliferation rate to zero, i.e. r = 0. The full system reads In turn, we conduct a non-dimensionalization of eqns. (40) to identify the variables. Moreover, this helps us to gain a knowledge about the relationships between the different model parameters.
By assuming that the system size is fixed at L the non-dimensional quantities read In the limit V 1, the eqn. (40) can also be written as To understand the behaviour of the system at long times, we conduct a fixed point analysis.
Initially we assume a well-stirred system, i.e. no spatial interactions. Then, eqs. (42) can be written as coupled non-linear ODEs which have three fixed points The above imply a pitchfork bifurcation for theβ = β 8V parameter, i.e. there exists a critical valueβ c that introduces the bistable state. From eqs.(43), we can easily deduce that the critical sensitivity value This is an acceptable approximation even for the finite interaction radius system (see next section). In the following, we analyze if the system is able to produce spatial patterns. Applying linear stability analysis for the spatially resolved system eqs. (42), we can deduce that no pattern formation is possible (for details check next section). Any perturbations to the homogeneous state lead always to a spatially homogeneous steady state. When r = 0 this result is consistent with our findings in discrete IBM simulations where very large interaction radii do not confer any clustering, as shown in Fig. 3b.

The finite interaction radius case: emergence of a bistable switch between "fluid" and "solid" tissue phases and pattern formation
Now we turn to the full system for intermediate interaction radius also (assuming proliferation). This implies a finite interaction volume V and for analytical feasibility we are interested in the Gaussian approximation of the switching probabilities and their corresponding mean-field terms. The full system of PDEs assuming also proliferation reads Since closed expression of steady states are not analytically feasible, we obtain the bifurcation diagram numerically, see Fig (4a). We observe the existence of a supercritical pitchfork bifurcation and the existence of a critical β c . For β ≥ β c , the systems depart from balanced state (ρ 0 = ρ 1 = 1 2 ) to the coexistence of a "fluid" (most cells migrate) and a "solid" (most cells are resting) phase. The switch is controlled by the perturbation on the ratio of migratory and resting cells. Interestingly, we can compare the analytic estimate of β c = V , from the V 1 case, with the one calculated for the finite V system. Fig. (4b) shows that the infinite system approximation provides an upper bound for β c , which is not too far from the real value of finite systems.
The existence of a critical β c with respect to the two phases is evident in the IBM simulations as well. In particular if we quantify the ratio of stationary over motile cells, we observe a similar behavior of the critical sensitivity for increasing V , since it decreases (see Fig.13 in SI ). However, we cannot conduct a strict quantitative comparison since our bifurcation analysis does not involve any diffusion as opposed to the IBM simulation. In turn, we apply linear stability analysis to identify parameter regimes that promote diffusion driven pattern formation or Turing instability. To analyze the Turing instability [20] we have to find the system's steady state, (i.e. when diffusion is not present in the systems of eqn.(40)). It has been shown that d = D 1 D 0 1 is a necessary condition for the emergence of spatially heterogeneous solution i.e. patterns. Now, we can write the system of PDEs (i.e eqn.(40)) in a generalized matrix form where ρ is defined as (ρ 0 − ρ * 0 , ρ 1 − ρ * 1 ) T and R is the Jacobian at ρ * = (ρ * 0 , ρ * 1 ). Using the Turing conditions of instability [20], we found patterns when N is finite in zero-flux boundary conditions. We have checked all the Turing instability conditions: Interestingly, only when the density of resting cells is larger than that of the moving cells, patterns are formed under Turing instability conditions. In order to investigate the system's potential to exhibit pattern formation, we check the range of validity of the Turing instability criteria (47). Diffusion-driven instability conditions are satisfied for a large portion of the parameter space. In turn for these parameters, we calculate the critical diffusion coefficient d c .
For values d > d c we are able to observe patterns, The existence of d c is associated with the existence of a critical wavenumber k c [20].
In Fig. (5) we identify the parameter regime that allows us to observe patterns and then corresponding d c and k c .
By fixing the initial conditions to a cosine wave, we observe in one dimension the existence of regular spikes as shown in Fig (6). Please note that the exact pattern is sensitive to the initial conditions. In turn, we simulate our system in 2D and for the same parameters and we observe patterning in the form of dots. Finally, we investigate if the type of patterns changes for variations in parameters β and r (see Fig. (7)). If d > d c then we need a large domain or size of the system to observe the patterns. In the 2D case, we have observed the discoidal patterns of resting phenotype which resemble the 1D case. So, the radius of the circles of the patterns are increased if we fix the domain size. Moreover, we observe that in both dimensions  Figure 8. Two-dimensional simulation of the mean-field approximation for finite interaction radius system. The patterns are obtained for parameters that correspond to the first point (TS1) in the Turing space, as displayed in Fig.(7) (i.e. 1D and 2D) the critical spatial frequency increases with decreasing r . In comparison with the discrete IBM simulations, we have to state that the simulation clusters can be identified by the discoidal mean-field patterns. Under this statement, we find that: 1. There exists a critical β that allows for the emergence of patterns as in the IBM simulations.
2. For very low and very high interaction radius , we observe no patterns, which is consistent with our IBM results (see Fig.3b).

Microenvironment independent phenotypic switching leads to uncontrolled growth dynamics
Finally we investigate the case β = 0 where cells do not sense their microenvironment and so we end up with the following system, We can do a fixed point analysis similar to equation (42). Assuming no spatial dynamics, we can find two fixed points i.e.
For further details see Fig. 9 in S.I. We can clearly show that the fixed point (0,0) is unstable and 1 2 , 1 2 is a saddle point. We can also write the coupled PDE equations in a single PDE as well, if we consider ρ 0 = ρ 1 =ρ. Then we obtain Now we can clearly see that the eqn.(51) is similar to the Fisher-Kolmogorov equation [20]. It is known that the Fisher-Kolmogorov equation does not exhibit any pattern formation instabilities. The latter observation is consistent with our discrete IBM simulations, where no clustering is observed.

Discussion
Recently, the Least microEnvironmental Uncertainty Principle (LEUP) has been proposed as an organization principle for cell decision-making in multicellular systems. In this paper, we apply this principle to shed light on the effects of phenotypic plasticity for tissue dynamics with a mathematical model. We focus on two types of plasticity: a go or rest, and a go or grow phenotypic dichotomy, which play key roles in important processes in biological development and pathological situations as cancer invasion. We assume that in any given spatial mesoscopic sample the presence of cell phenotypes follows a binomial distribution. Using this assumption, we are able to calculate the microenvironmental entropy and the LEUP-driven probability distribution of each phenotype. On the basis of this distribution, we defined an appropriate microscopic stochastic model for the spatio-temporal dynamics of both phenotypes, and in turn derived the corresponding mean-field description resulting in a system of coupled reactiondiffusion equations. The main results of our study are: (i) in the case of go or rest plasticity, there exists a supercritical pitchfork bifurcation that defines a switch between a "fluid" and a "solid" phase, and (ii) in the case of go or grow plasticity, for the "solid" phase, we can derive conditions for the emergence of Turing patterns. Interestingly, the existence of Turing patterns requires a critical LEUP sensitivity to the microenvironment and a minimum interaction range.
Our model assumes binary transitions between two discrete phenotypes, which is an extreme simplification of biological reality. Certainly, one could include a continuous spectrum of motile/proliferative phenotypes that is expected to imply even richer spatio-temporal dynamics. Indeed, assuming a continuous state (velocity) space would potentially lead to further interesting bifurcations, such as the metastable EMT state as found in [21,22], or complex spatiotemporal patterns as indicated in [23]. But the very abstracted case chosen here is still instructive to indicate how entropy-driven phenotypic decisions can lead to particular types of spatio-temporal pattern formation.
Interestingly, our LEUP-driven IBM is an extension of a Vicsek-type model [24], formulated in the context of self-propelled particles [25]. Our model exhibits a novel collective behavior when compared to the past published results from Vicsek-type of models. In particular, to our knowledge it is the first time to produce with such models Turing patterns, i.e. dynamics clusters of non-motile cells of specific characteristic wavelength. Typically, in Viscek-type models we may observe moving clusters of swirling cells (e.g. the milling Viscek model) but never static ones.
At this point, let us focus on the biological assumptions and implications of our study. The molecular regulatory mechanisms involved in EMT or GoG remain largely unknown, where the latter can be viewed as an EMT with proliferation constrained to the epithelial/resting phase. Here, we assume that the phenotypic regulation of both mechanisms is based on the minimization of microenvironmental entropy in physiological tissues, which allows us to predict the multicellular spatiotemporal dynamics. This assumption is supported by the fact that healthy physiological processes, biological development or processes like wound healing, where EMT/MET or GoG are present, typically lead to an ordered (low entropy) tissue from a disordered initial condition. On the other hand, deregulation of EMT and GoG have been already identified as pivotal elements in invading cancers [9,26], where genetic and phenotypic heterogeneity, characterized by high entropy, is a key characteristic. Assuming a LEUPdriven migration/proliferation phenotypic regulation allows us to understand how cells control multicellular dynamics in terms of growth and patterning.
Implications in multicellular growth control: The central finding of our study is associated with the bifurcation diagram of the LEUP sensitivity parameter β in Fig. 4a. As stated before, the parameter β quantifies how prone cells are in sensing and responding to their microenvironmental stimuli. When β = 0 cells migrate and proliferate independent of their microenvironmental cues, which corresponds to one of the cancer hallmarks [27]. In this case, the systems grows uncontrollably, resembling a cancerous tissue [27]. The resulting Fisher-Kolmogorov macroscopic behavior has been prototypically used to model invading tumours [28]. On the other hand, by adding a death process for any motility state in the GoG model, we can recapitulate the Allee effect (bistability between extinction and growth) as found in Boettger et al [29]. A cell population at the "fluid" state will go extinct (motile cells do not proliferate but still have a probability to die,) where as systems in the "solid" state will always grow until carrying capacity. Therefore, in the bistability regime cell sensing properties lead to multicellular growth control.
Implications in multicellular pattern control: Increased cell sensing β represents the physiological tissue dynamics, since it allows the system to control its behavior. By tuning β and the ratio of motile/resting cells, the system exhibits a bistable behavior between a "fluid"/mesenchymaltype and a "solid"/epithelial-type tissue phase. This kind of tissue level switch is of utmost importance in physiological processes such as wound healing or embryogenesis [30]. After a tissue injury, the healing is characterized by a "fluid" diffusive expansion of fibroblast cells, that adopt a migratory phenotype via EMT [31]. After covering the wound, the "solid" phase emerges as cells stop the migration program and proliferate to finalize tissue repair. In the abscence of proliferation, the bistable switch from the "solid" to the "fluid" state could potentially explain the jamming phase transition observed in epithelial colonies, under EGF modulation and Rab5a knock-out [32]. Typically, the "solid" multicellular phase is prone to the emergence of pattern formation, which frequently occurs in physiological epithelial tissues [31,33]. When EMT is combined with a Notch-Delta cell-cell communication then epithelial/immotile cell clusters emerge [22], as observed in our IBM and mean-field simulations. Adding proliferation in the GoG model, the type and the size of such emerging patterns require a tight regulation of microenvironmental sensing β and proliferation rate r, as indicated by Fig. 11a, since the critical wavelength k c depends on the ratio β r . Such Turing patterns are in agreement with previous GoG studies [34], where Turing patterns where emerging.
In conclusion, our study shows how individual LEUP-driven cell decisions the dynamics at the tissue level and how knowledge of collective cell decision-making can be used to control of growth and pattern. Diffusion coefficient of resting cells D 1 Diffusion coefficient of migratory cells

S.1. Calculation of microenvironmental entropy
We assume that a maximum number of N cells is present inside the cell's microenvironment, where l is the radius of the microenvironment (Fig.1a). Since there is the possibility of free slots, if there are less than N cells in the cell neighbourhood, we have started by assuming a trinomial distribution for the i-th cell is P N 1 i , N φ i where the number N 0 i of cells having phenotype X i = 0 and by the number N 1 i of cells having phenotype X i = 1. In addition, a number N φ i of empty slots are included inside the microenvironment. The joint probability is defined by where p and θ are the probabilities of having a number N 1 i of cells with phenotype (X i = 1) out of N cells and having a number N φ i of free slots. The conditional probability of having a number N 1 i of cells present in the microenvironment given a number N φ i of free slots is Now, we use the form of binomial distribution to calculate the entropy of the Binomial distribution. Please note that P (N 1 i | N φ i ) can be written as P (N 1 i | N T ) due to Binomial distribution. where According to LEUP we have to evaluate the microenvironmental entropy of S (Y i | X i = 0) and S (Y i | X i = 1) to calculate the probability of the internal states (X i = 0) and (X i = 1) From the Gaussian approximation we can write the entropy difference as

S.2. Calculation of cell proliferation rate
In our go-or-grow model only resting cells are allowed to proliferate but not the moving cells. The growth rate depends in a specific way on the number of moving and resting cells in the microenvironment. In particular, we assume that the per-capita growth rate for resting cells is a linearly decreasing function of N 0 i and N 1 i , and is also decreasing with the number of migratory cells and a constant per-capita death rate d 1 . Accordingly, where q r = 1 and r = h 1 − d 1 .   Throughout the simulations, we kept the total density at 0.2 and the total number of cells was 500.