Table Of Content

Spatial self-organization in hybrid models of multicellular adhesion Spatial self-organization in hybrid models of multicellular adhesion Bonforti, A.,1,2,3,4,a) Duran-Nebreda, S.,1,2,3,a) Montan˜ez, R.,1,2,3 and Solé, R. V.1,2,3,5,b) 1)ICREA - Complex Systems Lab, UPF, Dr Aiguadé 88, 08003 Barcelona, Spain 2)Institut de Biologia Evolutiva (CSIC-Universitat Pompeu Fabra), Passeig Mar´ıtim de la Barceloneta 37, 08003 Barcelona, Spain. 3)Department of Experimental and Health Sciences, UPF, 08003 Barcelona, Spain 4)CIDI - Sant Joan de Deu Research Foundation, 08950 Barcelona, Spain 5)Santa Fe Institute, 1399 Hyde Park Road, Santa Fe NM 87501, USA Spatial self-organization emerges in distributed systems exhibiting local interactions when nonlinearities and the appropriate propagation of signals are at work. These kinds of phenomena can be modeled with different frameworks,typicallycellularautomataorreaction-diffusionsystems. Adifferentclassofdynamicalprocesses involvesthecorrelatedmovementofagentsoverspace,whichcanbemediatedthroughchemotacticmovement or minimization of cell-cell interaction energy. A classic example of the latter is given by the formation of 6 spatially segregated assemblies when cells display differential adhesion. Here we consider a new class of 1 dynamical models, involving cell adhesion among two stochastically exchangeable cell states as a minimal 0 model capable of exhibiting well-defined, ordered spatial patterns. Our results suggest that a whole space 2 of pattern-forming rules is hosted by the combination of physical differential adhesion and the value of n probabilitiesmodulatingcellphenotypicswitching,showingthatTuring-likepatternscanbeobtainedwithout a resorting to reaction-diffusion processes. If the model is expanded allowing cells to proliferate and die in an J environmentwherediffusiblenutrientandtoxicwasteareatplay,differentphasesareobserved,characterized 2 by regularly spaced patterns. The analysis of the parameter space reveals that certain phases reach higher 1 population levels than other modes of organization. A detailed exploration of the mean-field theory is also ] presented. Finally we let populations of cells with different adhesion matrices compete for reproduction, B showing that, in our model, structural organization can improve the fitness of a given cell population. The C implications of these results for ecological and evolutionary models of pattern formation and the emergence . of multicellularity are outlined. o bi Keywords: Turing Patterns, Artificial life, Evolutionary Transitions, Cellular Automata, Complexity - q [ I. INTRODUCTION ralmechanismsofpatternformation16–20. Thestructures generated by these processes have a characteristic scale 1 v Theevolutionoflifeshowsanoveralltrendtowardsan whose wavelength depends on the model parameters. 8 increase in size and complexity1,2. One of the determin- Along with this class of pattern-forming mechanisms, 1 ing major innovations that have allowed biological sys- another possible class of models capable of organiz- 9 temstoachieveahighdegreeofcomplexityhasbeenthe ing structures in space is based on cell-cell differential 2 evolution of multicellularity and the emergence of supra- adhesion21,22. Such a mechanism explains the spatial re- 0 . cellular hierarchies beyond single-cell organization3. To- arrangement of different cells belonging to disrupted tis- 1 getherwithmulticellularity,mechanismstomaintainsta- sues when mixed together in vitro23. After a transient, 0 ble phenotypes that underly consistent division of labor clusters involving cells of the same class are often ob- 6 1 had to be developed4,5. served as spatially segregated from other cell types by : The study of the origins of form have a long tradi- meansoftheformationofwelldefinedboundariesorlay- v tion in biology6. Initiated by Turing7 and Rashevsky8, ers. In this case, the underlying mechanism explaining i X numerousattemptstoformalizeamathematicaldescrip- the origin of patterns is that of energy-minimization dy- r tion of pattern formation have been made. As a result, namics, similar to the one used in physics for strongly a spatialinstabilitieswereproposedasapowerfulrationale interacting particle systems. Both reported mechanisms for the creation of spatial order, out of random fluctua- are crucial in the formation of natural self-organized tions, around a homogeneous state in reaction-diffusion structuresindevelopingembryos24,25,andhavebeencon- systems9–13. The main feature of reaction-diffusion sys- nected to the early forms of multicellularity22,26. temsisthepresenceofdiffusion-driveninstabilitiesunder Inthispaperwefocusourattentionontheearlystages certain parametric conditions, by which small perturba- of the transition towards multicellularity, where the ex- tionsinthesystemareamplified,leadingtoorderedspa- plicit connection between fitness, function and structure tial patterns. This family of models has been systemati- has been particularly difficult to elucidate and, thus, is callystudied14,15 andprovidesthebasisforseveralnatu- commonly overlooked. To assess whether the structural organization of multicellular assemblies is related to differential fitness, we have developed an embodied com- putational model where Turing-like structures appear, a)Thesetwoauthorscontributedequallytothiswork stemming from differential adhesion and stochastic phe- b)Correspondingauthor notypicswitching. Fitnessisintrinsicallyobtainedbythe Spatial self-organization in hybrid models of multicellular adhesion 2 introduction of a limiting nutrient and the production of the DA-SPS model, cells are sorted by differential adhe- a toxic waste byproduct, which respectively increase cell sion and can reversibly switch their phenotypes. In the reproduction or death. One of the two cellular states is ECmodelweincludeasimplemetabolismbyaddingnu- abletoprocesswasteatthecostofreducedproliferation. trient and toxic waste, whose concentrations drive cell We observe that different parameter sets produce differ- proliferation and death. In the following sections we ex- ent spatial patterns, and that spatial organization can plain how cellular adhesion, phenotypic switching and have a role in increasing fitness. Finally, we discuss the metabolism are implemented in our models. implications of such results for the transition from unicellular to multicellular organisms and for the evolution of complexity. A. Differential cellular adhesion a b 50 The cell sorting process fundamentally occurs due to the differences in adhesion energy between states. 40 Following Steinberg’s differential adhesion hypothesis (DAH) we assume that the adhesion kinetics are driven 30 by the minimization of adhesion energy between lattice sites, being cells more or less prone to remain to- 20 gether, and avoid or maximize contact with the external medium23,25,27,29. The strength of interactions among 10 different states can be defined by means of an adhesion matrix : 0 0 10 20 30 40 50 J   c kp J(σ0,σ0) J(σ0,σ1) J(σ0,σ2) J =J(σ1,σ0) J(σ1,σ1) J(σ1,σ2) . � � 1 2 J(σ2,σ0) J(σ2,σ1) J(σ2,σ2) Each term J in this matrix describes how favourable 1 kp 1 kq (a,b) � kq � the pairwise interaction between two states is. The matrix is symmetric, i.e. = , and has = J(a,b) J(b,a) J(σ0,σ0) 0 always. To avoid confusion, we will use the notation FIG. 1. (a) Schematic representation of the cell sorting al- σ when we refer to a given state, and the notation S gorithm as described by Steinberg27. Cells of different types n ij to indicate a state occupying a given lattice site (i,j). coexistinaregularsquarelatticealongwithemptymedium, The underlying idea here is that cells will tend to move swapping positions with random neighbouring cells if a par- whenever this allows the system to reach a lower energy. ticular potential H is minimized. (b) The final global con- figuration is a direct consequence of the micro rules imposed It can be shown that the energy function in a given H bytheadhesionmatrix. (c)Markovianprocessmodellingthe position (i,j) can be defined as follows: cellstatetransitionsusedinthispaper. Thisverysimpleap- (cid:88) proach can aptly describe persister cell dynamics and phase = , (1) variation phenomena32. Hij JSkl,Sij Skl∈Γij where Γ is the set defined by the eight nearest neigh- ij boursofacellinposition(i,j)(Moore’sneighbourhood), eachofwhichoccupiesaposition(k,l),andhasadefined II. THE MODEL state S . To calculate the probability that the cell in kl (i,j) will swap with a randomly chosen neighbour, we Our model considers a population of cells living on calculatetheenergyfunctionwhennoswapoccurs. This a two-dimensional square lattice (Fig. 1), along with energyfunction,named ∗,consistsoftwoterms,onein- emptymedium,andfollowingtherulesofacellularPotts H volvingthecellinitsoriginalpositionanditsneighbour- model as described by Steinberg27 (see also28). Within hood set Γ , and another involving the cell’s neighbour, ij thisframework,cellsarediscreteentitiesthatoccupysin- located in i(cid:48),j(cid:48), with its neighbourhood Γi(cid:48),j(cid:48). We then gle lattice positions, have an associated state (σ ) and n virtuallyswapthepositionsofthecellwithitsneighbour moveacrossthelatticetryingtominimizetheirenergetic ini(cid:48),j(cid:48), andcalculatetheenergyfunctionwhenswapoc- potential. Two states correspond to cellular phenotypes, curs ( (cid:48)). The energy difference is then defined as: namely white cells (σ ) and black cells (σ ), while state H 1 2 σ0 represents empty space. ∆ = (cid:48) ∗. (2) In this paper we build and analyse two different ex- H H −H pansions of the basic Potts model: a hybrid differen- When the difference is negative, a decrease in the global tial adhesion-stochastic phenotipic switching (DA-SPS) energyoccursandthestateswillswapposition. Instead, model and an ecology and competition (EC) model. In when∆ >0, thelargerthedifferencethelesstheswap H Spatial self-organization in hybrid models of multicellular adhesion 3 FIG.2. Effectsofstochasticswitchingonpatternformationresultingfromdifferentialadhesion. Inthefirstrowwedisplaythe patternsformedafter105 iterationsofBoltzmannupdaterulesina40×40thoroidallatticewith0.35occupation. Thesecond rowdisplaystheresultingpatternsofthesameconditionswhenswitchingbetweentypesisimplemented(p=q=1,κ=10−3). Byaddingastochasticswitchingrule,onecanobservechangesinthecharacteristicsizeofthestructures,andalsointhespatial arrangement of cell types. Finally, in the third row we show the adhesion matrix for each simulated condition. islikelytohappen,withaprobabilityfollowingtheBoltz- B. Stochastic phenotypic switching mann distribution. If we indicate P(S S ) as the ij kl → probabilitythatourcellmovesfrom(i,j)to(k,l), itcan Cells can perform reversible transitions between their be shown that: states σ and σ , similarly to phase variation and persis- 1 2 tence in natural bacterial populations30–32. Switching is regulated by transition probabilities p and q, 1 P(Sij →Skl)= 1+e∆H/T, (3) P(σ1 →σ2)=κpij P(σ2 →σ1)=κqij where κ is a fixed scaling factor, introduced to regulate the relative speed between adhesion kinetics and pheno- typicswitching. BysimplyaddingSPStoaclassicalDA wheretheparameterT isanoisefactoractingasa‘tem- model, cell sorting properties can change drastically for perature’,essentiallytuningthedegreeofdeterminismof some adhesion matrices (see Fig. 2). It is worth men- our system. The Boltzmann factor e∆H/T acts in such a tioningthatinSPSthetransitionsbetweenstatesarenot way that if ∆ = 0, the probability of swapping is 1/2. H dependentonanymolecularcuenoranycellularmemory Note that cell-cell (and cell-medium) interactions are lo- beyond their current state. cal (Fig. 1), meaning that a cell in (i,j) interacts only with the set Γ of its eight nearest neighbours. ij C. Metabolism Depending on the form of adhesion matrix, different patternscanbeobservedinacellsortingsystemwithtwo Cellularmetabolismisdefinedbytwosimplepathways: celltypes. Unlessotherwiseindicated,inoursimulations the ability of both σ (white) and σ (black) phenotypes wewillapplyasymmetricadhesionmatrix(seeFig. SM- 1 2 totransformnutrientN (constantlyaddedtothelattice) 1), wherecellstendtoattachpreferentiallytoothercells into cellular energy E and waste byproduct W: inthesamestate,andsecondarilytocellsoftheopposite state,whileattachmentwithemptyspaceisnotfavoured. N ρ◦ E +W (4) ij ij ij −→ ρε• Thevaluesoftheadhesionmatrixdeterminethestruc- Nij Eij +Wij, (5) −→ tureofpatternsformedbythecellsortingalgorithm(Fig. and the unique ability of σ cells to degrade waste: 2 2), which can be perturbed by the effects of phenotypic ρ(1−ε)• switching. W . (6) ij −→ ∅ Spatial self-organization in hybrid models of multicellular adhesion 4 Theσ cellscanallocateresourcesforwastedegradation, 2 at the cost of reduced energy production and therefore 1 proliferation, following a linear trade-off (1 ε 0) uninhabited ≥ ≥ consistent with a maximum metabolic load and shared 0.8 resources for protein synthesis. Therefore, the temporal dynamicsofmetabolitesN ,E ,W foragivenposition ij ij ij (i,j) in the lattice are described by: 0.6 q diffusionterm inputrate ∂N (cid:122) (cid:125)(cid:124) (cid:123) (cid:122)(cid:125)(cid:124)(cid:123) 0.4 ij =D 2N + µ N ij ∂t (cid:53) ( η + ξδ + εξδ )N (7) − N (σ1,Sij) (σ2,Sij) ij 0.2 (cid:124)(cid:123)(cid:122)(cid:125) (cid:124) (cid:123)(cid:122) (cid:125) (cid:124) (cid:123)(cid:122) (cid:125) decay absorptionbyσ1 absorptionbyσ2 fully inhabited 0 intakebyσ1 intakebyσ2 decay 0 p0.2 0.4 0.6 0.8 1 ∂E (cid:122) (cid:125)(cid:124) (cid:123) (cid:122) (cid:125)(cid:124) (cid:123) (cid:122) (cid:125)(cid:124) (cid:123) p ij =(ξδ + εξδ )N η E (8) ∂t (σ1,Sij) (σ2,Sij) ij − E ij diffusionterm σ1-produced σ2-produced FIG. 3. Separation between fully inhabited and uninhabited ∂W (cid:122) (cid:125)(cid:124) (cid:123) (cid:122) (cid:125)(cid:124) (cid:123) (cid:122) (cid:125)(cid:124) (cid:123) domainsinthe(p,q)phasespaceinthecomputationalmean- ij =D 2W +(ξδ +εξδ )N ∂t W (cid:53) ij (σ1,Sij) (σ2,Sij) ij field model, for ε = 0.7. See Fig. (SM-4) to see how the separation slope varies at different values of ε. (η +(1 ε)ξδ )W (9) − W − (σ2,Sij) ij (cid:124)(cid:123)(cid:122)(cid:125) (cid:124) (cid:123)(cid:122) (cid:125) decay degradationbyσ2 setofODEsthatconstitutethemetabolismofcellsinour Here, µ is the rate of input of the nutrient resource, system. As a starting point we use the waste differential and DN and DW correspond to the diffusion rate of nu- equation (9). In a well mixed scenario there is no spatial trient and waste respectively. It is worth noting that, structure and all variables are homogeneous. Hence, the being E an intracellular metabolite, it does not diffuse diffusion term and the position subindices cease to be of through the lattice. Variables ηN, ηW and ηE corre- relevance. Therefore, at the steady state we get: spond to the exponential decay parameter of each of the (cid:34) (cid:35) metabolites,whileξ definesthemaximalabsorptionrate. ξN∗ P +εP (1−ε)W∗P ηWW∗ =0, The trade-off parameter ε adjusts the proportion of nu- σ1 σ2 − N∗ σ2 − ξN∗ trient allotted to energy production or waste degrada- tion in σ cells. Taking into account spatial dynamics of where N∗ and W∗ are the equilibrium concentrations of 2 metabolites,bothkindsofcellscandieeitherduetolocal nutrient and waste respectively, and Pσn is the proba- excess of toxic waste or due to lack of internal energy: bility that a cell in the system has state σn. Here we separate the analysis into two solutions: , Wij≥Θ1  {◦ •} −→ ∅ ξN∗ =0 , Eij≥Θ2  (10) {◦ •} −→ ∅ W∗ =N∗(P +εP )/[η /ξ+(1 ε)P ] where , indicates that the process equally affects σ1 σ2 W − σ2 {◦ •} both types of cells. Θ1 is the maximum value of toxic W∗ is equal to zero when N∗ =0 (trivial unstable solu- waste a cell can sustain, Θ2 is the minimum value of in- tion)orwhenthesecondnullclineismet. Todevelopthe ner energy needed for survival, and Θ3 defines the inner mathematical treatment we assume that the two popu- energy threshold needed for a cell to divide, provided an lations of cell states σ and σ are at equilibrium. The 1 2 empty position exists in its neighbourhood: ratio between populations can be then deduced from the persister cell population dynamics: , Eij≥Θ32 , . {◦ •} −→ {◦ •} ∂P σ1 =ξNP +κ(qP pP ) η (W,N)P Motheranddaughtercellshavethesameproperties. En- ∂t σ1 σ2 − σ1 − σ1 σ1 ergy is equally split between the two cells after division. ∂P σ2 =ξεNP +κ(pP qP ) η (W,N)P ∂t σ2 σ1 − σ2 − σ2 σ2 Thethreetermsontheright-handsideoftheseequations III. MEAN-FIELD APPROACH represent reproduction, stochastic switching and death, respectively. At the steady state, given that: To better understand the general properties of the model,wedevelopedamean-field(MF)approachtothe ξNP η (W,N)P =0, σn − σn σn Spatial self-organization in hybrid models of multicellular adhesion 5 the relation between the two populations is: N∗,η andξ havenon-zero,positivevalues,andp,q and n εareconstrainedtotherange[0,1],itcanbeshownthat P =(q/p)P . σ1 σ2 this boundary, if represented in the (p,q) phase space, has a positive linear slope for ε = 0, which decreases SincewehaveP =P +P ,theexpectedprobabilities T σ1 σ2 non-linearly as ε increases, as shown in Fig. (SM-4). for each population at equilibrium are: We tested this MF prediction against the actual sim- q P = (11) ulations of the EC model, finding also in the latter a σ1 (p+q) boundary, now defining a sharp transition between full p P = (12) occupationandsparsestructures. Theslopeisverysimi- σ2 (p+q) larintheECandMFmodels(Fig. 6),indicatingthatthe Being cell death a threshold function, all cells will die if latter properly captures some features of the EC model. W∗ >Θ ,andintheoppositecasenocellwilldie. There- Equations (14) and (16) also show that there exists a 1 fore,whetherthepopulationwillreachfulloccupationor critical εc beyond which the slope is negative: not can be determined by incorporating Θ1 into eq. (1) (1+η /ξ)/(1+N∗/Θ ) (from eq. 14) and transforming the equation into an inequality:  N 1 ε = c N∗(q/p+ε) 1/(1+N∗/Θ ) (from eq. 16) 1 Θ1 > (cid:104) (cid:105) (13) 1 ε+ ηW (p+q) − ξ p Under such conditions no pair p,q in [0,1] can make the inequality true, resulting in a system dominated by Thisexpressiondefinestheregionoftheparameterspace death processes. In the simplified scenario described by in which, even at maximum population, W∗ < Θ and 1 inequality (16), ε is always in the range [0,1], but if we c cells do not die. Reordering the terms, we obtain: consider the non-simplified eq. (14), then ε can assume c (cid:34) (cid:35) a value out of the [0,1] range, when η Θ >ξN∗. Θ (1+η /ξ)/N∗ ε(1+Θ /N∗) N 1 1 N 1 q <p − . (14) (1 Θ η /ξN∗) 1 N − B. Pattern formation in the DA-SPS model This inequality defines the boundary dividing the inhabited from the uninhabited region in the (p,q,ε) phase In the DA-SPS expansion of Potts model, cells can space. In Fig. (SM-4) the boundary for different values perform phenotypic switching through a transition rule ofεisshown. Ifwastedegradationperformedbyσ cells 2 regulated by transition probabilities p and q, and their isfargreaterthanthepassivedecayterm( η W),then W − movement is driven by differential adhesion and the ten- the denominator of eq. (13) becomes: dency to minimize interaction energy between cells. Pa- ηW (p+q) rameter κ is introduced as a scaling factor which reg- 1 ε+ 1 ε − ξ p ≈ − ulates the relative speed at which cell sorting by adhesion and phenotypic switching occur. Using this model, hence eq. (12) gets simplified to: we assess the role in pattern formation of the individual N∗(q/p+ε) transition probabilities (p, q) for a fixed κ. This analysis Θ1 > (1 ε) (15) reveals that, in spite of model simplicity, σ1 and σ2 cells − are able to self-organize in space in periodic structures. and the inequality in eq. (14) becomes simpler: Thesepatternscanrangefromspotstostripestomazes, depending on the relative values of p and q (Fig. 4a). (cid:34) (cid:35) Θ Θ To characterize the phases of the morphospace we ap- 1 1 q <p ε(1+ ) . (16) N∗ − N∗ pliedastandardpercolationalgorithmtothefinalmacro- state of each simulation. Figure (4b) shows that the av- The concentration of nutrient at equilibrium is given by: erage reachable fraction of each of the two states σ and 1 µ σ2 displays a sharp transition. This defines three clear N∗ = . (17) regimes: non-percolating (spots), percolating (fully con- η +ξ(q+εp) N (p+q) nectedmaze)andtransitionregime, wherespotsbecome stripesofincreasinglengthandaremarginallyabletoex- tend to other domains. Interestingly, this transition oc- IV. RESULTS cursforadifferentvalueofthecontrolparameterwithre- specttonon-correlatedpercolationstudies33 (q =0.5in- A. Mean-field model steadofq 0.592). Forthisvaluethestructuresofboth c ≈ cell states percolate, giving rise to the labyrinth phase. Equations (14) and (16) show that the boundary that Although percolation analysis shows a sharp transition, in the MF model separates the inhabited from the unin- the number of domains -i.e. clusters of lattice sites with habited domain depends on p, q, ε and on the other pa- same state- varies smoothly over the phase space (Fig. rametersofthemodelappearingintheequations. AsΘ , 4c). Another interesting feature that can be observed in 1 Spatial self-organization in hybrid models of multicellular adhesion 6 stead, for high log (κ) values, SPS occurs at a faster 10 pacethanthecellsortingprocess,whichbringsaboutan almost random distribution of states. Following the same strategy as before, we ran a set of simulations with periodic boundary conditions, tran- sitionprobabilitiesp=q =0.5, anddifferentvaluesofκ. We then applied a standard Fourier Transform analysis in order to confirm the existence of dominant frequencies in the spatial distribution of lattice states after 106 iterations. Fig. (5b) displays the spatial frequency con- tributioninaradiallyaveragedpowerspectrum(RAPS). ThepresenceofasinglepeakintheRAPSindicatesthat the periodic structures are built by a single dominant frequency without any specific orientation in the spatial domain, i.e. we obtain fixed wavelength structures rem- iniscent of Turing patterns. Furthermore, the peak position and width are subject to the particular value of κ, specifically the wavelength of the pattern decreases as log (κ) increases. The particular mathematical relation 10 between these two variables is displayed in Fig. (5c). C. Phase transitions and fitness landscapes In the ecology and competition (EC) model cells are not only subject to adhesion processes and phenotypic switching, but can absorb the substrate which is con- stantly produced all over the lattice, and transform it intoenergyonly(σ cells)orenergyandwaste(σ cells). 1 2 We compare the EC model with the mean-field (MF) model, to understand which properties can be predicted from the latter and which ones instead emerge from the complexity introduced by spatial organization and inho- mogeneities in the levels of waste and substrate. We run a set of simulations, each with a different pair of p and q values, using periodic boundary conditions and fixing ε=0.7 and k =10−3. The κ parameter is fixed at 10−3, FIG. 4. (a) Pattern formation on a fully occupied lattice in such a way that adhesion processes occur 103 times inwhichtheindividualtransitionprobabilitiesvaryspatially fasterthanphenotypicswitching. Intheinitialtimestep, from 0 to 1, κ = 10−3, T = 10, with a symmetric adhesion 100% of lattice is occupied by randomly distributed σ matrixandafter105 iterations. Differentarrangements-from 1 andσ cells, beingtheratiobetweenthetwophenotypes spotstostripestomazes-ofbothcelltypesareattained,rem- 2 already set at the equilibrium value following (Eq. 11, iniscentoffixedwavelengthstructures. (b)Averagereachable fraction(F)ofcellsofaparticulartype(whiteandblackcir- 12),dependingonpandq values. TheresultsinFig. (6) cles represent σ and σ cells respectively), fixed p=0.5. (c) andinFig. (SM-2)showaclearcorrespondencebetween 1 2 Average domain count (D) -defined as groups of contiguous the analytical MF result and the simulated MF, indicat- cells with the same state- for each pair (p, q). ing that the model we implemented properly reproduces the expected theoretical results. In the EC model, cells can also occupy the region of the phase space which was Fig. (4a, 4c) is that even if the ratio p/q remains con- left empty in the MF model. In fact, while in the MF stant, lower values of the transition parameters generate model simulation all cells die instantaneously when the bigger structures with fewer domains. levelofwastereachesΘ (themaximumlevelofW acell 1 Since in our model κ is a scaling factor for both p and cansustain),intheECmodelthedeathofafewcellscan q, we set out to quantify the effect of this parameter reduce the pressure on the system and allow population in the pattern formation process. Fig. (5a) shows the survivalalsointhoseparametricallyunfavouredhabitats qualitative impact of the tuning of κ in a full-occupation where waste concentration can locally exceed Θ . 1 DA-SPS model. At low log (κ) values, SPS occurs at a Theslopeinthe(p,q)phasespaceseparatinginhabited 10 slower pace than the cell sorting process, which is there- fromuninhabitedregionintheMFmodelvariesdepend- fore able to properly separate cells in two phases. In- ing on the value of ε and of other relevant parameters of Spatial self-organization in hybrid models of multicellular adhesion 7 FIG. 5. (a) The effect of changing the coupling factor between adhesion processes and stochastic phenotype switching. Here p=q=0.5,andlog (κ)istunedcontinuouslyfrom−1to−5. (b)Dominantwavelengthforthesamerangeofκvaluesusing 10 a radially averaged power spectrum. For discrete values of κ (in hues of blue) we show the existence of a moving peak in the frequency domain, with decreasing frequencies -i.e. with increasing wavelength- as κ decreases. (5 replicates of a 250×250 lattice,t=106 algorithmiterations,processedwithMatlab2013bFFT).(c)Averagevalueofthedominantwavelengthplotted against log (κ). We also show a possible fitting of the data. 10 the model (Fig. SM-4). For ε=0, black and white cells structural organization can affect the fitness of cells in areindistinguishable,andtheirbehaviourisindependent the habitat. To do so, we launch a set of simulations in fromthevalueofpandq,ascanbecalculatedintheMF which cells with different adhesion matrices compete for analytical model. As ε increases, σ cells can degrade reproduction. At the beginning of the simulation 10% 2 waste better, but are also less able to elaborate nutrient of the lattice is occupied by cells, which are divided in and can die from lack of energy, unless p and q are set in equal populations differing only by type of adhesion ma- suchawaythatthereisahighprobabilityfortheσ cell trix, which can be aggregate, trabecular, symmetric, null, 2 to switch to a σ cell before inner E value gets below Θ onions, sponge or unicellular (see Fig. SM-1). The ad- 1 2 (the minimum value of E a cell needs for survival). hesion matrix type is transmitted by each cell to its off- Atsteadystatedifferentstructuresemerge, depending spring. Cells are randomly distributed over the lattice on p, q and ε values. The results for ε=0.7 described in regardless of the population they belong to. It is impor- thisparagraphareshowninFig. (6)andinFig. (SM-2). tant to stress that in terms of preferential attachment, For high values of p and q, σ1 and σ2 cells form Turing cellsonlysensethestatesσ1 andσ2 ofneighbouringcells, Patterns in a percolating structure. For high values of independently from population type. p and low values of q, σ2 cells dominate (note that for Foreachsetofparametersε,pandqinthephasespace, ε=0.7theyarestillabletodegradeWsoasnottoreach we assess which type of adhesion matrix brings the re- Θ1). Forlowervaluesofpandhighervaluesofq,σ1 cells lated population to maximum fitness, by comparing the substitute σ2 cells less rapidly, Θ1 can be reached more level of occupation of the lattice for each population. In easily and cells start to die. Fig. (7a) we represent cells of state σ or σ indepen- 1 2 Lastly, we want to assess whether our EC model, inte- dently from their population, while in Fig. (7b) we show grating DA, SPS and metabolism in a habitat with nu- the same cells differentiated by population. At the bot- trient and toxic waste, can present situations in which tom (Fig. 7c, 7d) we represent the change of occupation Spatial self-organization in hybrid models of multicellular adhesion 8 1 1 1 0.98 1 0.8 n n o o ati 0.96 ati0.75 0.6 p0.9 p 0.5 u u c 0.94 c 0.4 oc oc0.25 0.92 0.2 0.8 0 1 1 1 1 0q.7 0.4 0.4 p0.7 0.9 0q.7 0.4 0.4 0p.7 0 0.1 0.1 0.1 0.1 FIG. 6. Occupation levels in the Ecology and competition (EC) model (left) and mean-field (MF) simulations (right), with ε = 0.7 and p and q varying from 0.1 to 1. The occupation level for each pair of values (p,q) is obtained from a separate simulationina150×150thoroidallattice. Theredlineinbothimagesdividesthefullyinhabitedregionfromtheuninhabited or partially inhabited region of the (p,q) phase space. Mean-field simulation reproduces the MF analytical model results, as expected. In the EC model simulation, only few cells die in the partially inhabited region and occupation level remains high becausethespatialinhomogenityofmetabolitesallowscells(forthechosenparameters)tosurviveandovercomethepresence of peaks of toxic waste. The boundary at which the transition happens changes depending on the value of ε, as can be shown by the analytical MF model, the slope becoming smaller at higher values of ε (see Fig. SM-4). As an example we show the patternsofsomelatticesforgivenvaluesofpandq theECmodel(whiteforσ cells,blackforσ cells,grayformedium). An 1 2 alternative representation of the same result is shown in Fig. (SM-2). level in time for each population. In the region of the V. DISCUSSION AND PROSPECTS phase space which is fully inhabited in the MF model, cells can survive independently from their adhesion ma- trixvaluesandneverdie. Forthisreasonthepopulations In this paper we have shown a novel way of construct- relatedtothevariousadhesionmatricesareequallynum- ing periodic arrangements of cell types in the form of bered in this area -with slight differences due only to a hybrid differential adhesion and stochastic switching growth speed before the lattice becomes saturated- and process. This mechanism does not rely on differential randomlydistributedoverspace,withnoemergingstruc- diffusion (normally found between the activator and in- ture. However, in the domain of the phase space which hibitor species in canonical Turing-type systems) yet it was uninhabited in the MF, p and q have such values can create the same kind of structures in a predictable, that do not guarantee survival of cells. Since cell death scalable way. The key ingredients proposed at this level may occur in the EC model for this region of the phase include the differential adhesion hypotheses stemming space, here the values of the adhesion matrix do make from Steinberg’s work (that considers the minimization a difference and some species get selected over others. dynamicsassociatedwithasetofinteractingspinsorad- In particular, the population with ‘trabecular’ adhesion hesion strengths) and genetic switches following Marko- matrix prevails. Moreover, we can observe that in this vian stochastic dynamics, which are the source of cell area cells organize in a maze structure, exhibiting divi- diversity and the basis of some adaptive responses dis- sionoflabourwithinthesamepopulation. Lastly, inthe played by microbial populations34. The switching dy- frontier between the two zones various populations can namics can modify the types of patterns expected from coexist, with a prevalence for ‘onion’ adhesion matrix at the purely energy-driven scenario, thus indicating that high values of q. In Fig. (SM-3) we show the relative potential forms of phenotypic change can lead to addi- occupation levels of each population at varying p and q. tional richness of pattern forming rules. A range of spatially ordered structures is obtained displaying characteristic lengthscales. Beingboth key ingredients present in extant organisms, we consider that this simple mechanism might have been used originally (and might be Spatial self-organization in hybrid models of multicellular adhesion 9 a b 1 1 uni spg 0.8 0.8 oni 0.6 0.6 nul q sym 0.4 0.4 trb agg 0.2 0.2 med 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 p p c d 5 55x104 40 agreg 404.5 14 trab 144 mberofcellsx213 snosuypunnmoliilconmenglsle mberofcellsxNumberofcells213123...555123 atsnosurgypunnamroliibalcoegnmengrgelsgle Nu0 Nu00.5 0 2 4 6 8 10 000 1000 22000 3000 44000S50t0e0ps66000 7000 88000 90001100000 Stepsx103 Stepsx103 FIG. 7. Collage of 100 runs at fixed ε = 0.7 and varying p and q, where populations of cells differing only by their adhesion matricesareputincompetition,toassesswhichadhesionmatrixallowsthecorrespondingpopulationtogainhigherfitness. In a)werepresentthedispositioninthelatticeattheendofthesimulationforthetwocellphenotypesσ (whitecells,represented 1 in yellow), and σ (black cells, represented in blue), regardless of what cell population they belong to. In b) we represent the 2 disposition in the lattice at the end of the simulation for cells of each population, represented in different colors. In the fully inhabitedregionbelowtheboundary,nodeathoccurs,thepopulationsgrowatsimilarratesandsaturatethelatticeinsimilar proportions. In the region of the phase space above the boundary, cell death occurs, allowing the different cell populations to competeforlatticeoccupation. Inthisregion,thetrabecular populationshowsthehighestfitnessoverotherpopulations. The twophasesaredelimitedbytheslopeseparatingtheemptyfromthefulloccupationregioninthemean-fieldanalyticalmodel. Inc)andd)werepresentthechangeofoccupationlevelintimeforeachpopulation,at(p=1,q=0.6)andat(p=0.2,q=0.6) respectively. See Fig. (SM-3) for the occupation levels of each population at varying p and q. reproduced in the future by synthetic means) to create dently from its past and from its neighbour’s state, and regular structures in aggregates and colonies. it is not influenced by the values of the adhesion matrix. In relation to pattern formation dynamics, our hy- The second relevant aspect considered is how these bridadhesionmodeloffersanalternativewaytogenerate forming structures might be of benefit to a developing Turing Patterns, which were up to now directly related cooperative population in presence of nutrient resources with Turing’s RD mechanisms mediated by a diffusible and toxic agents. To do so we developed an ecology and molecule, or with apparently unrelated but mathemati- competition model where a minimal metabolism enables cally equivalent systems such as direct contact-mediated positive or negative interactions between cells. Cells can regulation by means of which cells are affecting each cooperate by metabolizing waste byproducts, yet they other’s internal rates of reactions35. In fact, differently willsufferfromdecreasedgrowthratesathigherpopula- from what was proposed by Babloyantz, in our hybrid tiondensitiesduetosubstrateattrition. IntheECmodel DA-SPS model the molecules on the surface of one cell further pattern-forming processes can be predicted. do not affect the rates of reactions in its neighbours: the Tofurtherasseshowstructuralorganizationcanaffect phenotypic switching process occurs in any cell indepen- the fitness of cells in the habitat, we studied how cells Spatial self-organization in hybrid models of multicellular adhesion 10 with different adhesion matrices compete for reproduc- 15Goodwin, B. C. 1994. How the leopard got its spots. Princeton tion. Interestingly when many populations differing in U.Press,Princeton. terms of adhesion properties compete, in the region of 16Koch, A.J. and Meinhardt, H. 1994. Biological Pattern Forma- the phase space with strong selective pressures only one tion: fromBasicMechanismstoComplexStructuresRev.Mod- ernPhysics66,1481-1507. of the populations survives (Fig. 7). The selected specie 17Lejeune, O., and Tlidi, M. 1999. A model for the explanation consistently develops a periodic multicellular structure of vegetation stripes (tiger bush).JournalofVegetationscience, which is superior to both the unicellular and the un- 10(2),201-208. structured multicellular one, suggesting that higher or- 18Sick,S.,Reinker,S.,Timmer,J.,andSchlake,T.2006.WNTand DKKdeterminehairfolliclespacingthroughareaction-diffusion derpropertiesmightbeofrelevancetotheestablishment mechanism.Science,314(5804),1447-1450. of functionality and cooperation. 19Economou,A.D.,Ohazama,A.,Porntaveetus,T.,Sharpe,P.T., Thissimplecompetitionmodelshowshowminimalin- Kondo,S.,Basson,M.A.andGreen,J.B.2012.Periodicstripe teraction properties pervading the metabolism of multi- formation by a Turing mechanism operating at growth zones in ple species might come to play a central role in forcing the mammalian palate.Naturegenetics,44(3),348-351. 20Raspopovic, J., Marcon, L., Russo, L., and Sharpe, J. 2014. thetransitiontocollectivefitnessandbehaviour,andsets Digit patterning is controlled by a Bmp-Sox9-Wnt Turing net- the groundwork for explicitly evolutionary automata36, workmodulatedbymorphogengradients.Science,345(6196),566- where cells can optimize several genotype dimensions in 570. order to attain more resources. 21Goel, N., Campbell, R. D., Gordon, R., Rosen, R., Martinez, H.,andYcas,M.1970.Self-SortingofIsotropicCells.J.Theor. Biol.28(3): 423-68. 22Newman, S. A. and Baht, R. 2008. Dynamical patterning mod- VI. ACKNOWLEDGEMENTS ules: physico-geneticdeterminantsofmorphologicaldevelopment and evolution.Phys. Biol.5: 015008. We thank the members of the Complex Systems Lab 23Forgacs, G., and Newman, S. A. 2005. Biological physics of the developing embryo.CambridgeU.Press,Cambridge. for useful discussions, and Amad´ıs Pagès for useful hints 24Forgacs,Gabor,etal.Viscoelasticpropertiesoflivingembryonic on code debugging. This work has been supported by tissues: a quantitative study. Biophysical journal 74.5 (1998): the Bot´ın Foundation by Banco Santander through its 2227-2234. Santander Universities Global Division, a MINECO fel- 25Foty,R.A.,andSteinberg,M.S.2005.Thedifferentialadhesion lowship and by the Santa Fe Institute. hypothesis: a direct evaluation.Dev. Biol.278(1): 255-263. 26Newman, Stuart A., and Ramray Bhat. Dynamical patterning 1Bonner, J. T. 1998. The evolution of complexity by means of modules: a”patternlanguage”fordevelopmentandevolutionof natural selection.PrincetonUniversityPress.Princeton. multicellularform.InternationalJournalofDevelopmentalBiol- 2Carroll,S.B.2001.Chance and necessity: the evolution of mor- ogy53.5(2009): 693. phological complexity and diversity.Nature409,1102-1109 27Steinberg, M S. 1975. Adhesion-Guided Multicellular Assem- 3Nedelcu, A.M. and Ruiz-Trillo, I. 2015 (eds.) Evolutionary bly: A Commentary upon the Postulates, Real and Imagined, Transitions to Multicellular Life: Principles and Mechanisms. of the Differential Adhesion Hypothesis, with Special Attention Springer-Verlag,London. toComputerSimulationsofCellSorting.J.Theor.Biol.55(2): 4Bonner, J. T. 2001. First signals: the evolution of multicellular 431-43. development.PrincetonUniversityPress.Princeton. 28Graner,F.,andJ.A.Glazier.Simulationofbiologicalcellsorting 5Newman, S. and Comper, W. D. 1990. Generic physical mech- using a two-dimensional extended Potts model. Physical review anisms of morphogenesis and pattern formation. Development letters69.13(1992): 2013. 110,1-18. 29Hogeweg, P. 2000. Evolving Mechanisms of Morphogenesis: on 6Waddington,C.A.1956.PrinciplesofEmbriology.GeorgeAllen the Interplay between Differential Adhesion and Cell Differenti- Ltd.London. ation.J. Theor. Biol.203: 317-333 7Turing, Alan. 1952. The chemical basis of morphogenesis. Phil. 30Hallet,B.2001.PlayingDrJekyllandMrHyde: combinedmech- Trans.RoyalSoc.LondonB237,37-72. anismsofphasevariationinbacteria.Currentopinioninmicro- 8Rashevsky, N. 1948. On periodicities in metabolizing systems. biology.4(5),570-581. The Bulletin of mathematical biophysics.10(3),159-174. 31BalabanNQ,MerrinJ,ChaitR,KowalikLandLeiblerS.2004. 9Murray,J.D.1981.Apre-patternformationmechanismforan- Bacterial persistence as a phenotypic switch. Science 305:1622- imal coat markings. Journal of Theoretical Biology, 88(1), 161- 1625. 199. 32Darmon, E., and D.R.F. Leach. 2014. Bacterial Genome Insta- 10SoléR.,BascompteJ.andVallsJ.1993.StabilityandComplexity bility.Microbiology and Molecular Biology Reviews: MMBR78 of Spatially Extended Two-species Competition. J. Theor. Biol. (1): 1-39. 159,469-480. 33Malarz,KrzysztofandGalam,Serge2005.Square-latticesiteper- 11Bascompte J. and Solé R. 1995. Rethinking Complexity: Mod- colationatincreasingrangesofneighborbonds.Phys.Rev.E71 elling Spatiotemporal Dynamics in Ecology. Trends Ecol. Evol. 34LewisK2010.Persistercells. Annu.Rev.Microbiol.64: 357-372. 10,361-366. 35Babloyantz, Agnessa. Molecules, dynamics, and life: An intro- 12Theraulaz, G., Bonabeau, E. Nicolis, S. C., Solé, R. et al 2002. duction to self-organization of matter.Wiley-Interscience,1986. Spatial patterns in ant colonies.Proc. Natl. Acad. Sci USA.99: 36Duran-Nebreda, S., Bonforti, A., Montan˜ez, R., Valverde, 9645-9649. S., and Solé , R. 2015. Emergence of proto-organisms from 13JilkineA.andEdelstein-Keshet,L.2011.AComparisonofMath- bistable stochastic differentiation and adhesion. arXiv preprint ematical Models for Polarization of Single Eukaryotic Cells in arXiv:1511.02079. Response to Guided Cues.PLoSComput.Biol.7,e1001121. 14Meinhardt, H., and Gierer, A. 2000. Pattern formation by local self-activation and lateral inhibition.Bioessays,22(8),753-760.