Morphogenesis 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 first-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.