Table Of ContentMorphogenesis by coupled regulatory networks
Thimo Rohlf and Stefan Bornholdt
Interdisciplinary Center for Bioinformatics, University of Leipzig, Kreuzstraße 7b, D-04102 Leipzig, Germany
(Dated: February 9, 2008)
Based on a recently proposed non-equilibrium mechanism for spatial pattern formation [13] we
studyhowmorphogenesiscan becontrolledbylocally coupleddiscretedynamicalnetworks,similar
to gene regulation networks of cells in a developing multicellular organism. As an example we
study the developmental problem of domain formation and proportion regulation in the presence
of noise and cell flow. We find that networks that solve this task exhibit a hierarchical structure
4 of information processing and are of similar complexity as developmental circuits of living cells.
0 A further focus of this paper is a detailed study of noise-induced dynamics, which is a major
0 ingredient of the control dynamics in the developmental network model. A master equation for
2 domain boundary readjustments is formulated and solved for the continuum limit. Evidence for a
n firstorderphasetransitioninequilibriumdomainsizeatvanishingnoiseisgivenbyfinitesizescaling.
a A second order phase transition at increased cell flow is studied in a mean field approximation.
J Finally, we discuss potential applications.
6
1 PACS:87.18.Hf, 05.70.Ln, 05.50.+q
]
N I. INTRODUCTION also with some nodes in the neighboring cells. Boolean
M networks are minimal models of information processing
in network structures and have been discussed as mod-
Understanding the molecular machinery that regu-
.
o els of gene regulation since the end of the 1960s [8, 9].
lates development of multicellular organisms is among
i The model of Jacksonet al. demonstrated the enormous
b the most fascinating problems of modern science. To-
pattern forming potential of local information process-
- day, a growing experimental record about the regula-
q ing, similar to direct contact induction [10] as known in
tory mechanisms involved in development is accumulat-
[ developmental biology. More recently, Salazar-Ciuadad
ing, in particular in well-studied model-organisms as,
et al. introduced a gene network model based on contin-
1 e.g., Drosophila or Hydra [1, 2]. Still, the genomic de-
v uous dynamics [11, 12] and coupling their networks by
tails known today are not sufficient to derive dynami-
4 direct contact induction. Interestingly, they observe a
cal models of developmental gene regulation processes
2 larger variety of spatial patterns than Turing-type mod-
in full detail. Phenomenological models of developmen-
0 els with diffusive morphogens, and find that patterns
1 tal processes, on the other hand, are well established to-
are less sensitive to initial conditions, with more time-
0 day. Pioneering work in this field was done by Turing,
independent (stationary) patterns. This matches well
4 who in his seminal paper in 1952 [3] considered a purely
0 physico-chemical origin of biological pattern formation. the intuition that networksofregulatorshave the poten-
/ tial for more general dynamical mechanisms than diffu-
o His theory is based on an instability in a system of cou-
siondrivenmodels. Let us here thereforestudy interact-
i pled reaction-diffusion equations. In this type of model,
b ing networksin patternformation and in particular con-
for certain parameter choices, stochastic fluctuations in
- sider information-transfer-based processes. In [13], we
q the initial conditions can lead to self-organization and
introduced a simple stochastic cellular automata model
: maintenance of spatial patterns, e.g. concentration gra-
v of spatial pattern formation based on local information
i dients or periodic patterns. This principle has been suc-
X transfer,thatperformsdenovopatternformationbygen-
cessfully applied to biological morphogenesis in numer-
erating and regulating a domain boundary. In the fol-
r ous applications [4, 5]. However, as current experiments
a lowing we will study how this very general mechanism
makeuswonderabouttheastonishinglyhighcomplexity
could emerge as a result of interacting nodes in coupled
of single regulating genes in development [6], they also
identicalnetworks,similartogeneregulationnetworksin
seem to suggestthat diffusion models will not be able to
interacting cells.
capturealldetailsofdevelopmentalregulation,andpoint
at a complex network of regulating interactions instead. The goalofthis paperis twofold: Firstwewillexplore
The role of information processing in gene regulatory ways to map the earlier cellular automata model [13],
networks during development has entered the focus of onto dynamical networks(Booleannetworks and thresh-
theoretical research only recently. One pioneering study old networks), motivated by the general observation of
was published by Jackson et al. [7] in 1986, who investi- gene regulationnetworks in biological development (sec-
gatedthedynamicsofspatialpatternformationinasys- tionIII).Secondly,wewillgiveamoredetailedaccountof
tem of locally coupled, identical dynamical networks. In thebasicunderlyingmechanism,focusingonthestochas-
thismodel,generegulatorydynamicsisapproximatedby tic dynamics (sectionIV) and on how the model was de-
Booleannetworkswith asubsetofnodes communicating velopedusingageneticalgorithm(Appendix): Insection
not only with nodes in the (intracellular) network, but IV, the statistical mechanics of noise-induced dynamics
2
is studied in detail. Numerical evidence for a first order tion regulation) and to establish the head-foot polarity.
phase transition at vanishing noise rate is given. A mas- Suchregulationofpositioninformationisaquitegeneral
ter equation for boundary readjustments is formulated problem in biological development [16]. An interesting
and solved in the continuum limit. Robustness against problem is how the specific properties of this regulation
cell flow is studied numerically and in a mean field ap- can be achieved by a small network of regulatory genes
proximation. IntheAppendix,thegeneticalgorithmap- andifso,whetherlocalcommunicationbetweenthecells
pliedforidentificationofmodelsolutionsisexplainedand (networks) is sufficient. This basic question is the cen-
statistical properties of different solutions are compared. tral motivation for the present study. In particular, we
considerthesimplifiedproblemofregulatingonedomain
as, for example, the foot region versus the rest of the
II. PROBLEM AND MODEL body. We consider this as a one dimensional problem as
first approximation to the well-defined head-foot-axis in
Hydra.
Let us first motivate and define the morphogenetic
problem that we will use as a test case. Then we un-
dertake a three-step approach to find a genetic network
modelthatsolvesthispatternformationproblem. Inthe head
firststep,wesummarizethepropertiesofthecellularau-
tomata model introduced in [13]. Cellular automata as
dynamical systems discrete in time and state space are
knownto display a wide varietyof complex patterns [14]
1−α
and are capable of solvingcomplex computationaltasks,
including universal computation. We searched for solu-
tions(i.e. ruletableswhichsolvetheproblem)byaidofa
rump
geneticalgorithm(fordetails,seeAppendix). Candidate
solutions have to fulfill four demands: Their update dy-
namics has to generate a spatial pattern which 1) obeys
a predefinedscalingratioα/(1 α),2)is independent of
− Ped
theinitialconditionchosenatrandom,3)isindependent α
of the system size (i.e. the number of cells N ) and 4) is
C foot
stationary (a fixed point). In the second step, this cel- Ped/CnNK2
lular automata rule table is “translated” into (spatially
coupled) Boolean networks, using binary coding of the FIG.1: GeneexpressiondomainsinHydra,herefortheexam-
cellularautomatastates. Thelogicalstructureoftheob- pleofthe“foot”genesPedibin(Ped)andCn-NK2. ThePed-
tained network is reduced to a minimal form, and then, ibin domain shows aquitesharp boundarytowards thebody
in step three, translated into a threshold network. region of the animal. The relative position of the boundary,
given by the ratio α/(1−α), is independent of the absolute
size of theanimal (proportion regulation).
A. A test scenario: Hydra foot formation
Developmental processes exhibit an astonishing ro-
A classical model organism for studies of position de- bustness. This often includes the ability of de novo pat-
pendent gene activation is the fresh water polyp Hy- tern formation, e.g., to regenerate a Hydra even after
dra, which has three distinct body regions - a head with complete dissociationofthe cellensemble in a centrifuge
mouth and tentacles, a body column and a foot region. [17]. Further, they are robust in the face of a steady
The positions of these regions are accurately regulated cell flux: Hydra cells constantly move from the central
along the body axis. body region along the body axis towards the top and
The problem we shall focus on is sketched schemati- bottom, where they differentiate into the respective cell
callyinFig. 1fortheexampleofthe“foot”genesPedibin types according to their position on the head-foot axis.
and CN-NK2 [15]. CN-NK2 is active only in the lowest The global pattern of gene activity is maintained in this
foot region, the expression domain of Pedibin extends a dynamic environment. Let us take these observations as
bit more towards the rump region, where it is turned off a starting point for a detailed study how the interplay
at some point, i.e. the domain shows a sharp boundary. of noise-induced regulatory dynamics and cell flow may
The relative positionof the boundary is almostindepen- stabilize a developmental system.
dent of the animal’s size, i.e. the ratio α/(1 α) (as
−
denoted in Fig. 1) is almost invariant under changes of
body size. As the Pedibin/CN-NK2 system presumably B. One dimensional cellular automata: Definitions
plays an important role in determining the foot region,
this invariance appears to be an essential prerequisite To define a model system that performs the pattern
for maintaining the correct body proportions (propor- formation task of domain self-organization[13], consider
3
on N , i.e. we are looking for a set of solutions where
i(cid:0)1 i i+1 C
the ratio of the domain sizes r := α/(1 α) is invari-
−
ant under changes of the system size. This clearly is a
non-trivial task when only local information transfer is
G1 G2 A(cid:27) (cid:0)(cid:0)(cid:2)
z z AU z - zpp p(cid:2)(cid:13)(cid:2) z zp allowed. The ratio r is a global property of the system,
JRRJJ^1 (cid:16)(cid:17)(cid:26)(cid:19)(cid:19)(cid:16)(cid:17)(cid:26)(cid:19)(cid:26)>(cid:16)(cid:17)(cid:19)(cid:19)7(cid:16)(cid:17)(cid:16)(cid:17)PiQkPQ(cid:16)(cid:17)PQ(cid:16)(cid:17)(cid:16)1(cid:17)3PQPQSoPQSISPQSPQSR(cid:2)(cid:3)3 (cid:0)(cid:1)(cid:17)+(cid:26)=(cid:26)(cid:17)(cid:26)(cid:17)(cid:17)(cid:17)(cid:17)(cid:17) winhteicrhachtiaosnstobeetmweeregnetfhreomcelpl’usrsetlaytelos.cal (next neighbor)
R2 R4
C. Translation into spatially coupled Boolean
networks
FIG. 2: Diagram showing the interaction structure of the
minimal network needed to solve the asymmetric expression Onecannowtakeastepfurthertowardsbiologicalsys-
task. Forthesakeofclarity,intracellularinteractionsbetween tems, by transferring the dynamics we found for a cellu-
the two genes G1 and G2 are shown only for cell i, and like- larautomatachainontocells ina line that communicate
wiseoutgoingintercellularsignalsfromthetwogenestwothe
with each other, similar to biological cells. Identifying
neighborcells i−1andi+1wereleftout. Thetranscription
states with cell types and assuming that the model cells
factors produced by gene G1 and G2 in cell i−1 couple to
haveanetworkofregulatorsinside,eachofthemcapable
thereceptorsystemsR1 andR2,respectively,whereas incell
i+1thetranscriptionfactorsproducedbythesegenescouple toreproducetherulesofacellularautomaton,weobtain
tothereceptorsystems R3 and R4 (biased signaling). Incell amodelmimickingbasicpropertiesofabiologicalgenetic
i, the receptors release factors which regulate the activity of network in development.
G1 and G2. Cellular automata rule tables can easily be translated
intological(Boolean)networks,e.g.,forn=3,twointer-
nalnodescanbeusedforbinarycodingofthecellstates.
a one-dimensional cellular automaton with parallel up- One then has
date [14]. N cells are arranged on a one-dimensional
C
lattice, and each cell is labeled uniquely with an index (σi(t),σi(t))=f (σi−1(t 1),σi (t 1),σi+1(t 1)) (3)
i 0,1,...,N 1 . Eachcellcantakenpossiblestates 1 2 1 1,2 − 1,2 − 1,2 −
C
∈{ − }
σi 0,1,..,n . The state σi(t) of cell i is a function with f : 0,1 6 0,1 2. The so obtained rule tables
∈ { } 1
of its own state σi(t 1) and of its neighbor’s states are, by ap{plica}tio7→n o{f Bo}olean logic, transformed into a
−
σi−1(t 1) and σi+1(t 1) at time t 1, i.e. minimized conjunctive normal form, which only makes
− − −
use of the the three logical operators NOT, AND and
σi(t)=f[σi−1(t−1),σi(t−1),σi+1(t−1)] (1) ORwithaminimalnumberofANDoperations. Thisisa
ratherrealisticassumptionfor generegulatorynetworks,
withf : 0,1,...,n 3 0,1,...,n (acellularautomaton
{ } 7→{ } as the AND operation is more difficult to realize on the
withneighborhood 3). Atthesystemboundaries,forsim-
basisofinteractionsbetweentranscriptionfactors. Other
plicity we choose a discrete analogue of zero flux bound-
logical functions as, e.g., XOR, are even harder to real-
ary conditions, i.e. we set σ−1 = σNC+1 = const. = 0. ize biochemically [18]. The structure of the constructed
Other choices, e.g. asymmetric boundaries with cell up-
networkandits biologicalinterpretationisshowninFig.
date depending only on the inner neighbor cell, lead to
2. Note that the up-down symmetry of the body axis is
similar results. The state evolution of course strongly
broken locally by a spatially asymmetric receptor distri-
depends on the choice of f: for a three-state cellular au-
bution [19].
tomaton (n = 3), there are 327 7.626 1012 possible
≈ ·
update rules, each of which has a unique set of dynam-
ical attractors. As we will show in the results section,
D. Translation into spatially coupled threshold
n=3 is the minimal number of states necessaryto solve
networks
the pattern formation problem formulated above.
Now we can formulate the problem we intend to solve
asfollows: Findaset offunctions(updaterules)which, Perhaps the simplest model for transcriptional reg-
Tgivuepndaatreanstdeopms eivnoitlivaeFlsvtehcetosrys~σte=m(’sσ0d,y.n..a,σmNiccs−1t)o,awfiitxheidn uBloaotiloeannnneettwwoorrkkss,awrheertehlroegshicoalldfunnecttwioonrskas,reamsoudbesleedt boyf
point attractor with the property: weightedsumsofthenodes’inputstatesplusathreshold
h [20, 21]. They have proven to be valuable tools to ad-
~σ∗ := σi =2 if i<[α·NC] (2) dressquestionsassociatedtothedynamicsandevolution
σ =2 if i [α N ] of gene regulatory networks [22, 23, 24, 25, 26].
i C
(cid:26) 6 ≥ ·
Any Boolean network which has been reduced to its
where [] is the Gauss bracket. The scaling parameter minimizedform(here,theconjunctivenormalform),can
·
α may take any value 0 < α < 1. For simplicity, it is be coded as a dynamical threshold network, with mini-
fixed here to α = 0.3. Notice that α does not depend mally three hierachies of information processing (“input
4
layer”: signals from the neighbor cells at time t 1, conditionsbecomes arbitrarilysmallwith increasingsys-
−
“hidden layer”: logical processing of the signals, “out- tem size. Hence, the pattern self-organization in this
put layer”: states of the two “pattern genes” in cell i at system exhibits considerable robustness against fluctu-
time t). The genes’states now may take values σ = 1, ations in the initial conditions. The main mechanism
i
andlikewise forthe interactionweights onehas clij =±±1 leadingtostabilizationatα∞ =0.281isamodulationof
foractivatingandinhibitingregulation,respectively,and thetravelingvelocityoftherightphaseboundaryinFig.
cl =0 if gene i does not receive an input from gene j in 2suchthattheboundaryonaveragemovesslightlymore
ij
cell l. The dynamics then is defined as thanonecelltothe leftperupdate step,whereasthe left
boundary moves one cell to the right exactly every third
σji(t)=sign(fj(t−1)) (4) update step. the modulation of the right boundary can
beseenastheresultofinteractingphaseboundariesrem-
with
iniscentofparticleinteractions. This pictureof“particle
2 i+1 computation” is a useful concept also in various other
f (t)= cl σl +h (5) contexts [27].
j kj k j
k=1l=i−1
X X
for the “hidden” genes, where σl,k 1,2 is the state
k ∈ { }
of the kth output gene in cell l (there are no couplings
between the genes in the “hidden” layer). The threshold
h is given by
j
2 i+1
h = cl 2, (6)
j | kj|−
k=1l=i−1
X X
which implements a logical OR operation. For the “out-
put” (pattern forming) genes, one simply has
kk
in
f (t)= σ kk , (7)
k l− in
l=1
X
i.e. the weights are all set to one, and the (negative)
threshold equals the number of inputs kk that gene k
in
receives from the “hidden” genes (logical AND).
FIG.3: AtypicaldynamicalrunfortheAutomataasdefined
in Table I, here for a system of size NC = 250 cells (deter-
ministic dynamics, no noise), starting from a random initial
III. RESULTS FOR DETERMINISTIC
configuration. Timeisrunningonthey-axisfromtoptobot-
DYNAMICS
tom. Cells with state σi =0aredepictedin blackcolor, cells
with σi=1 in red and cells with σi =2 in blue.
A. Cellular Automata Model
The first major outcome of the cellular automata
model is that a number of n = 3 different states is nec- B. Interaction topology of the minimal network
essaryandsufficientfor this classof systems to solvethe
given pattern formation task. The update table of the
In this section let us translate the rule table found in
fittest solution found during optimization runs, which
the previous section into a (gene) regulatory network.
solves the problem independent of system size for about
We will see that this network has biologically realistic
98 percent of (randomly chosen) initial conditions (i.e.
properties regarding the number of genes necessary for
has fitness Φ = 0.98), is shown in Table I. Fig. 3 shows
informationprocessingandthe complexityofinteraction
the typical update dynamics of this solution. The finite
structure,makingitwellconceivablethatsimilar“devel-
sizescalingoftheself-organizedrelativedomainsizeαas
opmental modules” exist in biological systems.
afunction ofthe the numberofcellsN isshownin Fig.
C
4. Inthe limitoflargesystemsizes, αconvergestowards
1. Boolean representation
α∞ =0.281 0.001. (8)
± The ruletable ofTable 1 firstis translatedinto binary
The variance of α vanishes with a power of N , i.e. the coding, i.e 0 00, 1 01 and 2 10, this corre-
C
→ → →
relative size of fluctuations induced by different initial sponds to two “genes” G and G one of which (G )
1 2 1
5
0.3 Gi2(cid:0)1 L
i+1
:G2
0.25 i(cid:0)1
:G1 L
i(cid:0)1 i+1
:G2 L :G1
0.2 0.1
) )℄ 0.01 Gi2 Gi2(cid:0)1
C C
(cid:11)(N 0.15 (cid:11)(N 00.0.000011 Gi1+1 L :Gi1 L
ar[ 1e-05 Gi2+1 :Gi1+1
0.1 V 1e-06 :Gi1(cid:0)1 :Gi2(cid:0)1
1e-07
1e-08 Gi1 L J Gi1(t) :Gi1 L J Gi2(t)
0.05 10 1000 100000 Gi2 Gi1+1
NC Gi1 Gi1(cid:0)1
0
10 100 1000 10000 100000 Gi2 L :Gi2 L
NC Gi2+1 Gi2+1
FIG. 4: Finite size scaling of the self-organized relative do- :Gi2(cid:0)1 L Gi1
main size αas a function of thetotal numberof cells NC. In
the limit of large system sizes, α converges towards a fixed Gi1+1 Gi2 L
value α∞ = 0.281±0.001 (as denoted by the straight line Gi1+1
fit). The inset shows the finite size scaling of the variance
i+1
Var[α(NC)]; the straight line in this log-log plot has slope G2
−1.3 and indicates that fluctuations vanish with a power of
FIG.5: Booleanrepresentationoftheminimalnetwork,min-
thesystem size.
imized conjunctive normal form. Gb with a ∈ {1,2} and
a
b ∈ {i−1,i,i+1} denotes gene a in cell number b. The in-
index σi−1 σi σi+1 σi index σi−1 σi σi+1 σi puts in the left branches of the trees are given by the genes’
0 0 0 0 0 14 1 1 2 2 statesat timet−1. ¬denotesNOT,⊙denoteslogical AND
1 0 0 1 2 15 1 2 0 0 and ⊕ logical OR.
2 0 0 2 1 16 1 2 1 0
3 0 1 0 0 17 1 2 2 1
4 0 1 1 2 18 2 0 0 0 involve a higher number of logical sub-processing steps,
5 0 1 2 2 19 2 0 1 0 i.e. a higher number of genes, hence we will not discuss
6 0 2 0 1 20 2 0 2 0 them here.
7 0 2 1 2 21 2 1 0 1 Considering the huge number of possible input con-
8 0 2 2 2 22 2 1 1 2 figurationswhichthe outputs theoreticallycould depend
9 1 0 0 0 23 2 1 2 2 on, the complexity of the resulting network is surpris-
10 1 0 1 1 24 2 2 0 1 ingly low. As shown in Fig. 5, the output state of gene
11 1 0 2 1 25 2 2 1 2
G only depends on five different input configurations
1
12 1 1 0 0 26 2 2 2 2
of at maximum four different inputs, gene number two
13 1 1 1 1
onsix different input configurationsof at maximum four
different inputs. This indicates that the spatial infor-
TABLEI:Ruletableofmostsuccessfulbestcellularautomata
mation flowing into that network is strongly reduced by
solution found during genetic algorithm search. In the left
column,theruletableindexisshown,runningfrom0to26,in internal information processing (only a small number of
themiddlecolumnthethreeinputstatesattimetareshown, input states leads to output “1”), as expected for the
the right column shows the corresponding output states at simple stationary target pattern. Nevertheless, this in-
time t+1. formation processing is sufficient to solve the non-trivial
task of domain size scaling.
is active only in a domain at the left side of the cell
chain. The so obtained Boolean update table is reduced 2. Threshold network implementation
toitsminimizedconjunctivenormalform,usingaQuine-
McCluskey algorithm [28]. For the construction of the Alternatively, the Boolean representation of the pat-
network topology we use the conjunctive normal form, tern formation system can be translated into a three-
as it is a somewhat biologically plausible solution with a layered threshold network (Fig. 6). The states of the
minimalnumberoflogicalANDoperations. Inprinciple, genes G and G at time t in a cell i and its two neigh-
1 2
other network topologies,e.g. with more levels of hierar- borcellsi 1andi+1serveasinputsof11information-
−
chy,arepossibleandbiologicallyplausible,however,they processing genes (“hidden” layer). The state of these
6
gene number spatial pattern
G i1−1 G 1i+1 1 111111111111111111111111111111
2 111111111110111111111111111111
3 111111111100000000000000000000
G i2−1 G 2i+1 4 111111111110000000000000000000
5 111111111110111111111111111111
6 111111111011111111111111111111
7 111111111110000000000000000000
8 100000000111111111111111111111
9 000000000111111111111111111111
10 111111111111111111111111111111
(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)G(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) 1(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)i (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)G(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) 2i(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) GG1112 110110110110110110110110110110101100100100100100100100100100100100100100100100100100100100
FIG. 6: Threshold network realization of the pattern forma- TABLE II: Spatial gene activity patterns of the network
tion system. Solid line arrows denote links with wij = +1, shown in Fig. 6. (stationary patterns after the system has
dashed arrows denote links with wij =−1. The inputs from reached equilibrium), NC = 30. The numbers identify the
the genes in the neighbor cells (Gi−1,Gi−1, Gi+1 and Gi+1) “hidden”genesin Fig. 6from theleft handsidetotheright.
1 2 1 1
are processed by a layer of “hidden genes” (hollow circles) States were transformed by σ→(σ+1)/2.
with different thresholds h (as these “genes” perform a log-
ical OR operation, one has h = kin −2). The threshold of
gene G1 is h1 = −5, for gene G2 one has h2 = −6 (logical
AND).
genes then defines the state of G and G in cell i at
1 2
time t + 2 (output layer). Additionally, there is some
feedback from G and G to the information process-
1 2
inglayer,asexpectedforthe dependenceoncell-internal
dynamicsalreadypresentinthecellularautomataimple-
mentationofthe model. The resulting stationaryspatial
patternsofgeneG (the“domaingene”),G (activeonly
1 2
at the domain boundary) and the “hidden” genes after
systemrelaxationareshowninTable2,forasystemsize
of N = 30 cells. From this seemingly redundant pat-
C
terninequilibriumone hardlyguessesthe higherlevelof
geneticinformationprocessingduringthe self-organizing FIG. 7: Quasi-particles, started by stochastic update errors,
phase. The network we construct here, regarded as a lead to control of the boundary position under noise (left
“developmental module” defining the head-foot polarity panel). The Γ particle (top right) leads to readjustment of
throughspatially asymmetric gene expression,has a size the boundary two cells to the left, the ∆ particle (bottom
similar to biologicalmodules (compare, for example, the right) leads to readjustment of the boundary one cell to the
right.
segment polarity network in Drosophila [29]) as well as
similar complexity (average connectivity K¯ 3).
≈
ppercell,leadingtoanaverageerrorrater =pN . In-
e C
IV. DYNAMICS UNDER NOISE AND CELL terestingly,this stochastic noisestarts moving“particle”
FLOW excitations in the cellular automaton which, as a result
indeed stabilize the developmental structure of the sys-
In the following we will study dynamics and robust- tem. To prepare for the details of these effects, define
ness of the model with respect to noise. Two kinds of first how we measure the boundary position properly in
perturbationsfrequently occur: Stochastic update errors the presence of noise. Let us use a statistical method to
and external forces cased by a directed cell flow due to measure the boundary position in order to get conclu-
cell proliferations. Both types of perturbations are very sive results also for high p: Starting at i = 0, we put
common during animal development, e.g., in Hydra cells a “measuring frame” of size w over cell i and the next
continuouslymovefromthecentralbodyregionalongthe w 1 cells, move this frame to the right and, for each i,
−
body axis towards the top and bottom, and differentiate measure the fraction z of cells with state σ = 2 within
into the respective cell types along the way according to the frame. The algorithm stops when z drops below 1/2
their position on the head-foot axis. and the boundary position is defined to be i+w/2.
Letusdefine stochasticupdate errorswithprobability It is easy to see that, for not too high p, there are
7
onlytwodifferentquasi-particles(i.e. stateperturbations
0.34
sw=200
moving through the homogeneous phases), as shown in
sw=1600
Fig. 7. In the following, these particles are called Γ and sw=12800
0.33
∆. TheΓparticleisstartedintheσ phasebyastochas- sw=51200
2
tic error σ = 2 σ = 2 at some i < αN ), moves to
i i C
→ 6 0.32
therightand,whenreachingthedomainboundary,read-
)
justs it two cells to the left of its original position. The re
(
∆particleisstartedintheσ0 phasebyastochasticerror (cid:3)(cid:11) 0.31
σ =0 σ =0 at some i>αN and movesto the left.
i i C
→ 6
Interaction with the domain boundary readjusts it one 0.3
cell to the right. Thus we find that the average position
∗
α of the boundary is given by the rate equation 0.29
∗ ∗
2α r =(1 α )r , (9)
e e
− 0.28
∗ 1e-07 1e-05 0.001 0.1
i.e. α = 1/3. Interestingly, for not too high error rates
∗ re
r , α is independent from r and thus from p. If we
e e
consider the average boundary position α∗ as a system- FIG. 9: Average domain boundary position α∗ as a function
specific order parameter which is controlled by the two oftheerrorratere,sampledoverupdatewindowsofdifferent
quasi-particles, then comparing the solution of Eqn. (9) lengthssw(ensemblestatistics,400differentinitialconditions
for each data point). The abscissa is logarithmic. With in-
tothe equilibriumpositioninthenoiselesscaseindicates
that the system u∗ndergoes a first order phase transition cdreeradsientgersmwi,ntishteictrdaynnsaitmioincsfrtoomα∗th=e1so/l3utuinodneαr∗dneoti=se0is.2s8h1ifutend-
with respect to α at p=r =0.
e towards re =0. The two straight lines definea lower bound-
∗ ∗
aryα andaupperboundaryα ,asexplainedinthetext.
low up
0.4
0.35 0.1
0.3 0.01 rreloeuwp((sss(cid:0)www1))
0.25
)
(cid:3)(re 0.2 s)w 0.001
(cid:11) (
ns
a
0.15 mean(cid:12)eldappr. NNCC ==110000 trre 0.0001
0.1 NC =400
mean(cid:12)eldappr. NC =400
0.05 NC =1600 1e-05
mean(cid:12)eldappr. NC =1600
0
1e-07 1e-05 0.001 0.1 10 1e-06
100 1000 10000 100000
re
sw
∗
FIG. 8: Average boundary position α as a function of the
FIG.10: Finitesizescaling oftheupperandlowertransition
error rate re for system sizes NC = 100, NC = 400, and pointsrup andrlow,i.e. thepointswhereα∗ crossesα∗ and
NarCe a=ver1a6g0e0d. ovTehre20a0bsdciiffsseareinstliongiatiraitlhcmonicd.itiNonusmweriticha2l ·d1a0ta6 α∗up,respeectivelye(Fig. 7), asa function of thesamplinlogwwin-
updateseach. Thedashedcurvesshowthemeanfieldapprox- dowlengthsw. Bothreup andrelow vanish∝s−w1,asindicated
bythe line with slope −1 in thislog-log-plot.
imationgivenbyEqn.(5),thestraightdashedlinesmarkthe
∗
unperturbedsolutionα =0.281andthesolutionundernoise,
∗
α =1/3.
twotransitionpoints rup and rlow areobtained. We find
e e
that rup c s−1 and rlow c s−1 with c > c
Figs. 9 and 10 show noise dependence and finite e ≈ up w e ≈ low w up low
as expected (Fig. 10), which implies that the difference
size scaling of the transition. In case of a first order
∆rtrans(s ):=rup rlow scales as
phase transition at re = 0, we would expect a shift of e w e − e
itsheprtorpanorsittiioonnalptooints−r1etraanss(wsewl)l atoswaarddisverrege=nce0 owfhtihche ∆retrans(sw)=(cup−clow)s−w1, (10)
slope at the transitiown point when s is increased, i.e. hence, because ∆α∗(rtrans) = const. = α∗ α∗ , in-
dα∗/dr (rtrans) when s .w deeddα∗/dr (rtrans)deivergeswhenthesamuppli−ngwloiwndow
The eshieft of r→tra∞ns(s ) iswm→ost∞easily measured by size goes to ienfienity.
e w
∗ ∗ ∗
defining a lower and a upper boundary α and α , The solution α =1/3 is stable only for 0<r 1/2.
∗ low up e ≤
respectively (Fig. 9); when α crosses these boundaries, AsshowninFig.7,theinteractionofaΓparticlewiththe
8
boundary needs only one update time step, whereas the
0.08
boundaryreadjustmentfollowinga∆particleinteraction NC =400
takes three update time steps. Hence, we conclude that 0.07 analyti
result
NC =800
the termonthe righthandside ofEqn.(10),whichgives
0.06 analyti
result
the flow rate of ∆ particles at the boundary, for largere NC =1600
will saturate at 1/3, leading to 0.05 analyti
result
(cid:11)) 100000
2α∗r = 1 (11) p( 0.04
e
3
0.03 t
with the solution 50000
0.02
1
α∗ = r−1+Θ(N ) (12)
6 e C 0.01 0
0.25 0.4 0.55
forre >1/2. Hence,thereisacrossoverfromthesolution 0
α∗ =1/3 to another solution vanishing with r−1 around 0.25 0.3 0.35 0.4 0.45 0.5 0.55
e
r = 1/2. The finite size scaling term Θ(N ) can be es- (cid:11)
e C
timated from the following consideration: for p 1, the FIG.11: Forthesystemwithstochasticupdateerrors,fluctu-
→
averagedomainsizecreatedby“purechance”isgivenby ationsoftheboundarypositionαaroundtheaverageposition
α∗ = N−1 NC (1/3)n n (3/4)N−1. If the measur- α∗ =1/3 are Gaussian distributed. The figure compares the
ingwindCowhans=s0izew,w·eo≈btainΘ(NC ) (3/4)wN−1. numericallyobtainedstationaryprobabilitydistributionwith
P C ≈ C theanalyticresultofEqn.(18)forthreedifferentsystemsizes.
To summarize, we find that the self-organized boundary
All data are gained for re = 0.1 and averaged over 100 dif-
position is given by ferent initial conditions with 2·106 updates each. The inset
shows a typical timeseries of theboundary position.
0.281 0.001 if r =0
e
∗ ±
α = 1/3 if 0<r 1/2 (13)
e
(1/6)re−1+Θ(NC) if re >≤1/2 The stationary solution of this equation is given by
with a first-order phase transition at re = 0 and a 3 3
crossoveraround r =1/2. f(x)= exp x2 , (17)
Now let us conseider the fluctuations of α around α∗ r2π (cid:20)−2 (cid:21)
given by the master equation i.e. inthe longtime limitt ,the probabilitydensity
→∞ ∗
fortheboundarypositionαis aGaussianwithmeanα :
pτ(α) = 2αr pτ−1(α+2δ)+(1 α)r pτ−1(α δ)
e e
− −
+ (N r )pτ−1(α) 2αr pτ−1(α)
C e e
− (1−−α)repτ−1(α) − (14) p(α,NC)= 32NπC exp −3N2C (α−α∗)2 . (18)
r (cid:20) (cid:21)
with δ = 1/N . Eqn. (15) determines the probability
C
pτ(α) to find the boundary at position α at update time FromEqn. (18)we see that the variance of α vanishes
∼
step τ, given its position at time τ 1. This equation 1/NC and the relative boundary position becomes sharp
canbesimplifiedasweareinterested−onlyinthe station- in the limit of large system sizes. Fig. 11 shows that
ary probability distribution of α. It is easy to see that this continuum approximation for NC 400 provides
≥
the errorrate r just provides a time scale for relaxation verygoodcorrespondencewiththe numericallyobtained
e
towards the stationary distribution and has no effect on probability distributions.
thestationarydistributionitself. Therefore,wemaycon- The stochastic nature of boundary stabilization under
sider the limit r rmax := N , divide through r and noise is also reflected by the probability distribution of
neglectthelastteh→reeteermsontCherighthandsideofeEqn. waiting times t for boundary readjustments due to par-
(14) (which become zero in this limit). We obtain ticle interactions: the particle production is a Poisson
process with the parameter λ=r and the waiting time
e
pτ(α)=2αpτ−1(α+2δ)+(1 α)pτ−1(α δ). (15) distribution is given by
− −
To study this equation, we consider the continuum limit p (t)=r exp( r t) (19)
wait e e
−
N . Let us introduce the scaling variables x =
f(α(Cx−,→t)α∗=∞)√NNCCp,τ(tα=NCτ)/.NICnsaenrtdintghethpesreobdaebfiinliittyiodnesnisnittyo wthiethwaanitianvgetriamgee dwiastitriinbguttiiomnsefhotrid=iffree−re2n.tFerirgo.r1r2atsehsorwes.
Eqn.(16)andignoringallsubdominantpowers (1/N ), In a biological organism, a pattern has to be robust
C
we obtain a Fokker-Planckequation: O not only with respect to dynamical noise, but also with
respect, e.g., to “mechanical” perturbations. In Hydra,
∂f(x,t) ∂2 ∂ e.g., there is a steady flow of cells directed towards the
= +3 x f(x,t). (16)
∂t ∂x2 ∂x animal’s headandfoot, due to continuedproliferationof
(cid:18) (cid:19)
9
0.01 1
re=0:01 re=0:1,left
analyti
result re=0:01,left
re=0:005 re=0:001,left
analyti
result 0.8 re=0:1,right
re=0:001 re=0:01,right
analyti
result re=0:001,right
) ) 0.6
p(tw 0.001 (cid:3)(cid:11)(rf
0.4
0.2
0.0001 0
0 200 400 600 800 1000 1e-07 1e-06 1e-05 0.0001 0.001 0.01 0.1 1
tw rf
FIG. 12: Probability distribution p(t) of waiting times for FIG.13: Averagedomainsizeα∗ asafunctionofthethecell
boundary readjustements in the model with stochastic up-
flow rate rf for three different error rates re; numerical data
date errors for three different values of re, semi-log plot. As (crossesandpoints)weresampledover10differentinitialcon-
expected for a Poisson process, p(t)is an exponential.
ditionsand1e6updatesforeach datapoint. “Left”indicates
cell flow directed to the left system boundary, “right” to the
rightsystemboundary,respectively. Thedashedlinesarethe
stem cells [30]; the stationary pattern of gene activity is corresponding solutions of Eqn.(20).
maintained is spite of this cell flow. Let us now study
the robustness of the model with respect to this type of
perturbation. Let us consider a constant cell flow with limitr 1theboundarypositionα∗detectedinnumer-
f
→
rate rf, which is directed towards the left or the right ical simulations deviates from the mean field prediction,
system boundary. In Eqn. (9), we now get an additional due to a boundary effect at the left system boundary
drift term rf on the left hand side: (stochastic production of finite lifetime stationary oscil-
lators,leadingtointermittentflowsofΓparticlesthrough
∗ ∗
2α re rf =re(1 α ), (20) the system).
± −
To summarize this part, we see that in the model
with the solution stochastic errors in dynamical updates for r > r in-
e f
deed stabilize the global pattern against the mechanical
α∗ = 31 1− rrfe if re ≥rf (21) stress of directed cell flow.
( (cid:16) 0 (cid:17) if re <rf
for the case of cell flow directed towards the left system V. SUMMARY AND DISCUSSION
∗
boundary (plus sign in eqn. (20)). One observesthat α
undergoes a second order phase transition at the critical
In this paper we considered the dynamics of pattern
valuercrit =r . Belowrcrit,thedomainsizeα∗vanishes,
e f e formation motivated by animal morphogenesis and the
andabovercrit itgrowsuntilitreachesthevalueα∗ =
e max largelyobservedparticipationofcomplexgeneregulation
1/3ofthesystemwithoutcellflow. Forcellflowdirected
networks in their coordination and control. We there-
towards the right system boundary (minus sign in eqn.
fore chose a simple developmental problem to study toy
(20)), we obtain
models of interacting networks that control pattern for-
mation and morphogenesis in this setting. In particular,
α∗ = 13 1+ rrfe if rf ≤2re (22) thegoalwastoexplorehownetworkscanofferadditional
( (cid:16) 1 (cid:17) if rf >2re. mechanismsbeyondthestandarddiffusionbasedprocess
of the Turing instability. Our results suggest that main
Inthiscase,thecriticalcellflowrateisgivenbyr =2r , functionsofmorphogenesiscanbeperformedbydynami-
f e
for cell flow rates larger than this value the σ -domain calnetworkswithoutrelyingondiffusivebiochemicalsig-
2
∗
extends overthe whole system, i.e. α =1. Fig. 13 com- nals, but using local signalingbetween neighboringcells.
paresthe results ofnumericalsimulations withthe mean This includes solving the problem of generating global
field approximation of Eqn. (20). In numerical simula- position information from purely local interactions. But
tions,cellflowisrealizedbyapplicationofthetranslation alsoitgoesbeyonddiffusionbasedmodels asitoffersso-
operatorΘσ :=σ toallcellswith0 i<N 1ev- lutions to developmental problems that are difficult for
i i+1 C
eryrf−1 timestepsandleavingσNC−1un≤changed. −Incase such models and avoids their inherent problem of fine
of cell flow directed to the right system boundary, in the tuned model parameters. It may be not too far fetched
10
to speculate that some developmental processes as, e.g., morphogenesis, the rate of cell divisions is controlled by
the establishment of position information, may rely on celldifferentiationandcell-to-cellsignaling. Dynamicsin
this type of internal information processing rather than both models, however, is deterministic. An extension of
oninterpretationof,e.g.,chemicalgradients. Asonlyin- our model as outlined above may open up for interest-
gredient to the present model, symmetry has to be bro- ing studies how stochastic signaling events could control
ken locally by a spatially asymmetric receptor distribu- and stabilize a global expressionpattern and cell flow as
tion, e.g., as proposed in [19], to provide the necessary an integrated system. Other possible extensions of the
input information. Also hybrid approaches are conceiv- model concern the dimensionality: In two or three di-
able. For example such broken symmetry could be pro- mensionsothermechanismsofsymmetrybreakingmight
vided by a gradient, even without fine-tuning (contrary be present, possibly leading to new, interesting dynami-
to standardgradient-basedmodels) since only the direc- cal effects.
tion of the gradient has to be read out. AJavaappletsimulationofthemodelcanbe foundat
Thenetworkmodelderivedhereperformsaccuratereg- [34].
ulation of position information and robust de novo pat-
tern formation from random conditions, with a mecha-
nism basedon localinformationtransfer ratherthan the VI. ACKNOWLEDGEMENTS
Turing instability. Non-local information is transmitted
through soliton-like quasi-particles instead of long-range We thank T.C.G. Bosch, T.W. Holstein, and U. Tech-
gradients. Two realizations as discrete dynamical net- nau for pointing us to current questions in Hydra devel-
works, Boolean networks and threshold networks, have opment. T.R. acknowledges financial support from the
been developed. The resulting networks have size and Studienstiftung des deutschen Volkes (German National
complexitycomparabletodevelopmentalgeneregulation Merit Foundation).
modules as observed in animals, e.g., Drosophila [29] or
Hydra [6]. The threshold networks (as models for tran-
scriptional regulation networks) process position infor- APPENDIX A: GENETIC ALGORITHM
mation in a hierarchical manner; in the present study, SEARCHES
hierarchy levels were limited to three, but realizations
with more levels of hierarchy, i.e. more “pre-processing” Let us briefly recapitulate here how the rule table of
ofinformationarealsopossible. Similar hierarchicaland the model has been obtained by the aid of a genetic al-
modular organization are typical signatures of gene reg- gorithm.
ulatory networks in organisms [18].
Robustness of the model is studied in detail for two
types of perturbations, stochastic update errors (noise) 1. Definition of the GA
and directed cell flow. A first order phase transition is
observed for vanishing noise and a second order phase
In order to find a set of update rules that solve the
transition at increasing cell flow. Fluctuations of the F
problem as formulated in section II, cellular automata
noise-controlledboundary position were studied numeri-
have been evolved using genetic algorithms [35]. Ge-
callyforfinitesizesystemsandanalyticallyinthecontin-
neticalgorithmsarepopulation-basedsearchalgorithms,
uum limit. We find that the relative size of fluctuations
whichare inspiredby the interplay ofrandommutations
vanisheswith1/N ,whichmeansthattheboundarypo-
C and selection as observed in biological evolution [36].
sition becomes sharp in the limit of large system sizes.
StartingfromarandomlygeneratedpopulationofP rule
Dynamics under cellflowis studied indetailnumerically
tables f , the algorithm optimizes possible solutions by
n
and analytically by a mean field approximation. A ba-
evaluating the fitness function
sicobservationisthatnoise-inducedperturbationsactas
quasi-particles that stabilize the pattern against the di- 1 Tu [α·NC]−1
rimseacatinseedcsiofzonerdcoeorroadfedcreoplmlhflaaoisnwe.terxAatnetnsaidticionrnigtitocovawelracretdlhlseflawovwhaonrliaesthesiy,nstgthedemroe-, Φ(fn)= (Tu−T)·NC t=XTu−T Xi=0 δσin(t),2
depending on the direction of cell flow, respectively. NC−1
Several extensions of this model are conceivable. In + {1−δσin(t),2}.(A1)
the present model, the cell flow rate rf is considered as i=[Xα·NC]
a free parameter, the global pattern, however, can be
The optimization algorithm then is defined as follows:
controlled easily by an appropriate choice of the error
rate re. This may suggest to extend the model by in- 1. Generate a random initial population =
troduction of some kind of dynamical coupling between f ,...,f of rule tables. F
1 P
r and r , treating r as a function of r . Interestingly, { }
e f f e
similar approaches have been studied by Hogeweg [31] 2. RandomlyassignsystemsizesNmin Nn Nmax
C ≤ c ≤ C
and Furusawa and Kaneko [32, 33]: In both models of to all rule tables.
