Spatial self-organization in hybrid models of multicellular adhesion

5.4 MB·

by  Adriano Bonforti

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´e, R. V.1,2,3,5,b)
1)ICREA - Complex Systems Lab, UPF, Dr Aiguad´e 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 dif-
ferential 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 uni-
cellular 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 lat-
tice 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 ma-
trix 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 inhab-
ited 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 adhe-
sion 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 frequen-
cies 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 po-
sition 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 spa-
tially ordered structures is obtained displaying charac-
teristic lengthscales. Beingboth key ingredients present
in extant organisms, we consider that this simple mech-
anism 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`es 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´eR.,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´e 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´e, 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´e , 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.