- Email: [email protected]
Receptor-based models with hysteresis for pattern formation in hydra
Mathematical Biosciences 199 (2006) 97–119 www.elsevier.com/locate/mbs
Receptor-based models with hysteresis for pattern formation in hydra Anna Marciniak-Czochra Institute of Applied Mathematics, University of Heidelberg, Im Neuenheimer Feld 294, 69120 Heidelberg, Germany Received 10 August 2005; received in revised form 27 August 2005; accepted 5 October 2005
Abstract In this paper, we propose a new receptor-based model for pattern formation and regulation in a freshwater polyp, namely hydra. The model is defined in the form of a system of reaction–diffusion equations with zero-flux boundary conditions coupled with a system of ordinary differential equations. The production of diffusible biochemical molecules has a hysteretic dependence on the density of these molecules and is modeled by additional ordinary differential equations. We study the hysteresis-driven mechanism of pattern formation and we demonstrate the advantages and constraints of its ability to explain different aspects of pattern formation and regulation in hydra. The properties of the model demonstrate a range of stationary and oscillatory spatially heterogeneous patterns, arising from multiple spatially homogeneous steady states and switches in the production rates. 2005 Elsevier Inc. All rights reserved. Keywords: Pattern formation; Hydra; Receptor-based model; Hysteresis
1. Introduction One of the crucial issues in developmental biology is to understand how coordinated systems of positional information are established during an organism’s development and how cells in the organism respond to the associated signaling cues, processes which ultimately result in the subdivided and patterned tissues of multicellular organisms [20,30,33]. E-mail address: [email protected] 0025-5564/$ - see front matter 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.mbs.2005.10.004
98
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
Experiments suggest that during development cells respond to local positional cues that are dynamically regulated. The hypothesis is that cells differentiate according to positional information [33]. The question is how this information is supplied to the cells. There exists a number of models for pattern formation and regulation based on the idea that positional information is supplied to cells by a diffusible biochemical morphogen [3,16]. Both regulatory and signaling molecules (ligands) act by binding and activating receptor molecules which are located in the cell membrane (or, with lipophilic ligands, in the cytoplasm) [20,18,11]. Thus in this paper we hypothesize that the positional value of the cell can be determined by the density of bound receptors which do not diffuse [18,19]. In this paper the test organism for mathematical modeling is a fresh-water polyp, the hydra, a simple organism which can be treated as a model for axis formation in higher organisms [20,30]. Hydra is an example of an organism in which the development of theoretical models has had strong influence on experimental design, starting with the positional information ideas of Wolpert [33] and the activator–inhibitor model of Gierer and Meinhardt [3,14]. Experiments performed on hydra provide a good ground for testing the abilities and limitations of such mathematical models.
2. A receptor-based model of pattern formation in hydra Hydra, a fresh-water polyp, has a tubular body about 5 mm long with a whorl of tentacles surrounding the mouth at the upper end and a disk-shaped organ for adhesion at the lower end. The longitudinal pattern is subdivided into a head, a gastric region, a budding zone (where new animals are generated by a process of natural cloning), a stalk and a foot. Experiments performed on hydra suggest that the function of cells is determined by their location. Following Mu¨ller [18], we assume that the positional information is determined by the density of bound receptors. Our models are based on the idea that epithelial cells secrete ligands (regulatory biochemical) which diffuse locally within the interstitial space (the space between two layers of cells lining the inner and outer body surface respectively) and bind to free receptors on the cell surface. A number of different receptor-based models for the evolution of spatial patterns have already been proposed, for example for the formation of cAMP patterns in Dictyostelium slime mold cells in [13] and [15]. The models considered here in this paper are simple compared with these aforementioned rather detailed and complicated schemes. This is due to the absence of a detailed biochemical knowledge of pattern control in hydra. Using simple and general models we try to clarify the ideas of the receptor-based mechanism of pattern control and explore their implications. The models we study are motivated by the model of Sherratt et al. and its receptor-based approach [27]. The objective is to build a receptor-based model without imposing an initial gradient as was done in the model of Sherratt et al. We investigate what kind of mechanisms regulating the dynamics of cell surface receptors and diffusing biochemical molecules can explain the observed results of experiments. We validate our models with the results of two crucial experiments: • Cutting experiments: After a cut, cells of the gastric region of the animal can reconstitute the missing structure, i.e. the head or foot, according to their position along the body axis. In hydra, regeneration occurs mainly by the re-patterning of existing tissues and is an example
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
99
of morphallaxis [33]. The lack of growth requirements for regeneration is shown in heavily irradiated hydra. No cell divisions occur in these animals, but they can still regenerate normally. Thus the mechanism of pattern formation in hydra seems to be independent of growth. Overlapping cut levels show that the same cells can form either the gastric region, or the head, or the foot, according to their position along the body axis (see Fig. 1). The most significant is the experiment which shows that after a transverse cut both parts of the animal can regenerate [18]. • Grafting experiments: Pieces of tissue are grafted from one animal to another and the outcome of such experiments depends on the change of the position along the body axes. Transplantation of the tissue from parts of the body column close to the head induces head formation while transplantation of the tissue from the position close to the foot results in foot formation [17] (see Fig. 2). In [12] we proposed a model describing receptor-ligand binding. The model takes into consideration the density of free receptors, of bound receptors and of ligands, as well as of a second diffusible substance functioning as an enzyme. We assume that new ligands and new free receptors are produced on the cell surface via a combination of recycling (dissociation of bound receptors) and de novo production within the cell. The ligands then bind reversibly to free receptors and bound receptors are internalized into the cell. Bound receptors also may dissociate. Both ligands and free receptors undergo natural decay. In addition, cells secrete an enzyme which diffuses along the body column and degrades the ligands. We consider a one-dimensional epithelial sheet of length L and denote the concentration of free receptors, bound receptors, ligands and enzyme densities by rf, rb, l and e respectively. For
Fig. 1. Cutting experiment. The hydra regenerates after a transverse cut of cells of the gastric region (from both upper and lower half of the body column). The experiment shows the availability of positional information cues. A hydra, as shown in the middle, is cut in two ways. In one experiment (left hand-side) the lower body column is removed, in the second experiment (right hand-side) the upper part is removed. The cut levels are not identical but somewhat different to show that one and the same group of cells (marked in grey) can form a foot (left hand-side), or a head (right handside) or a gastric segment (original state in the middle). The function of the cells depends on the position along the body column. [Courtesy of W. Mu¨ller].
100
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
Fig. 2. Grafting experiment. Determination of relative positional information values by transplantation. Pieces of tissue are grafted from one animal to another and one of three outcomes is observed. (1) If the tissue is transplanted from the upper position along the body column to the lower position then a new head is formed. (2) If the former and new position is the same then the piece is integrated and nothing is observed. (3) If the tissue is grafted to the upper position a new foot is formed. We can see that the disparities between the positional value of the transplant and the surrounding host tissue result in the head formation or the foot formation respectively. [Courtesy of W. Mu¨ller].
simplicity we assume that all binding processes are governed by the law of mass action without saturation effects. The model has the following form: o ð1Þ rf ¼ lf rf þ pr ðrb Þ brf l þ drb ; ot o ð2Þ rb ¼ lb rb þ brf l drb ; ot o o2 ð3Þ l ¼ d l 2 l ll l brf l þ pl ðrb Þ þ drb be le; ot ox o o2 ð4Þ e ¼ d e 2 e le e þ pe ðl; rb Þ ot ox with zero-flux boundary conditions for l and e. dl and de are diffusion coefficients for the ligands 2 and enzyme. By rescaling the system, we set dl to 1c and de to dc2 , where c ¼ Ld l and d 2 ¼ dd el and consider the model on the domain [0, 1]. The rates of decay of free receptors, bound receptors,
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
101
ligands and enzyme are denoted by lf, lb, ll and le respectively. d is the rate of dissociation, b is the rate of binding of ligands and free receptors and be is the rate of binding of ligands and the enzyme. By pr, pl and pe we denote production rates of new free receptors, ligands and enzyme. In [12], we considered different nonlinear production terms of free receptors and diffusing biological m1 rb , substances, using as an example Michaelis–Menten kinetics, which resulted in pr ðrb Þ ¼ 1þr b m3 rb m2 r b pl ðrb Þ ¼ 1þrb and pe ðrb Þ ¼ 1þrb . We considered also the possibility of production of the enzyme which would depend on the density of ligands (it could be regulated via some other receptors m3 l m1 rb not considered in our model). This resulted in the kinetic terms pr ¼ 1þr , pl = m2rb2 and pe ¼ 1þl . b In [12] we showed that a diffusion-driven (Turing-type) mechanism for pattern formation could be involved in a receptor-based model. This uses the special features of diffusion which results in the destabilization of the spatially homogeneous steady state and the formation of spatially heterogeneous patterns. In models involving this mechanism patterns can arise spontaneously (see the survey by Murray [16]). We note that the above system with diffusion-driven instability is able to simulate the cutting experiments. However, the model fails to reproduce the grafting experiments (compare Section 5.2). The reason for this is that the final heterogeneous pattern is stable to spatial perturbations such as those arising in models of transplantation. Transplantation of the tissue from one part of the body column of a host to another part of the body column of a donor corresponds to the changing of values of all the variables in that part of the domain. New values create a local maximum or local minimum depending on the relationship between the host and donor positions of the transplant. In the models with diffusion-driven instability the gradient pattern is stable for a given range of diffusion coefficients and a local perturbation of this state cannot lead to the creation of another type of spatial structure observed in the grafting experiments. Simulations of the model show, for example, that if two peaks are present in the initial conditions (which corresponds to grafting a new head to the lower part of the body column), the second peak will disappear after a short time and a gradient pattern will appear. Stable solutions with 2 peaks can only be obtained by increasing the size of the domain at the same time. This deficiency of models with the Turing-type mechanism of pattern formation encourages us to look for an alternative mechanism which could also explain the results of grafting experiments. 3. Modeling hysteresis As stated above, the results of grafting experiments on hydra suggest that the model should involve a memory-based relation. Thus far we have modeled the functions pl and pe as nonlinear functions depending on rb (and l in the case of pe) and then determined which nonlinearities the system allows a Turing-type pattern formation. It seems appropriate to consider pl and pe to be determined by a system of ordinary differential equations, which is a simplified model of the metabolic processes of the cells. As in many other cases the complexity of the biological process makes it difficult to give a rigorous derivation of these equations. We can only use heuristic reaction terms to depict and investigate the patterning process. Existing qualitative theory of reaction–diffusion systems is mainly focused on stability conditions for homogeneous steady states, see, for example [16,28]. It is very important to study the dynamics of systems functioning far from such equilibria. Here, we consider the mechanism of pattern formation which results from the existence of multiple steady states. Switching between
102
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
different states arises in many biological processes [9,26]. One such mechanism is hysteresis, which seems to be important in the modeling of biological development, since according to the experimental observations, inductive signals are present only during a certain time interval of development. The signal triggers changes in the cell’s nucleus and evokes differentiation which does not revert when the signal is stopped. Hysteresis results from the existence of multiple steady states (see Fig. 3). Suppose we have three steady states, two outer states which are stable, and one inner state which is unstable. The system may jump from the first stable state to the second when certain parameters are changed. If these parameters are changed back to their previous values, the system remains in the second state until the parameters reach another critical value. It is known that hysteresis kinetics may lead to spatial patterns in chemical and biological systems, such as Liesegang rings formed by precipitating colloids [21], bacterial growth patterns [6,7] and DNA melting/reannealing. Models based on a hysteresis-type bistability were also used in chemistry to describe enzymatic reactions [1,25,10]. Prominent and well-documented examples of hysteresis-driven mechanisms in biology are cell-cycle transitions [31], for example in Xenopus laevis egg extracts [26,22] and cell division cycles of fission yeast [23]. A mathematical model of cell-cycle progression predicts that irreverspl
5
4.5
C
4
stable 3.5
3
2.5
D
2
unstable
1.5
1
0.5
0
0
B
stable
A 0.5
1
l1
1.5
2
2.5
3
3.5
4
l
l2
Fig. 3. Nullclines of gl – hysteretic dependence of pl on l. The two outer branches are stable, the inner branch is unstable. If the system evolves along the stable branch AB, the density of ligands l increases and it reaches the value l2 where the system becomes unstable and it jumps to the other stable branch CD. Decreasing the value of l below l2 does not cause any change and the system remains on the upper stable branch till the value of l approaches l1 from above. The l value needs to drop below l1 to shift the system back to the lower stable branch AB.
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
103
ible transitions into and out of mitosis are driven by hysteresis in a molecular control system [26,32]. The first experimental demonstration of hysteresis in the eucaryotic cell engine has been reported recently [2]. 3.1. Additional equations modeling the production terms We propose to modify the model (1)–(4) and consider the production of diffusing substances described by new variables pl and pe, determined by additional equations, o p ¼ gl ðrb ; l; pl Þ; ot l ð5Þ o p ¼ ge ðrb ; l; e; pe Þ. ot e Here all the variables depend both on time t and space x, that is pl = pl(t, x) and so forth. However the evolution in time of pl and pe is prescribed in each point x by the same rule, given by the kinetics gl and ge. In this sense we refer to system (5) as a system of ordinary differential equations. Spatial changes in the dynamics of pl and pe are caused by coupling these dynamics with the dynamics of the ligands and enzyme, which diffuse spatially. Biological signal transduction mechanisms are complicated cascades of different biochemical reactions which are not, as yet, understood in complete detail. Therefore in this paper, we consider particular reaction kinetics which allow the testing of the possibilities and restrictions of our model (1)–(5). Such functions should be biologically realistic and fulfill some basic assumptions. We assume that: • Production of ligands and enzyme in a steady state has a hysteretic dependence on the amount of ligands and the enzyme present (see Fig. 4). • Bound receptors trigger the production of ligands and enzyme. • Production is controlled by its own size. To model the hysteretic dependence of pl on l and pe on e we take the following kinetic functions: pl m2 lrb þ ; 2 2 1 þ pl ð1 þ rl pl bl pl Þð1 þ al rb Þ pe m3 le . þ ge ¼ de 2 2 1 þ pe ð1 þ re pe be pe Þð1 þ ae rb Þ gl ¼ dl
ð6Þ ð7Þ
Fig. 4 shows the nullclines of gl crossing the nullclines of the kinetic system of Eq. (3). We show this in the coordinate system l, pl. The plot is drawn for rb = 1, but it remains qualitatively the same for all rb. We have chosen the parameters to obtain three solutions of the subsystem o p ¼ 0, oto l ¼ 0. ot l To model the production of the enzyme we use a similar function. It exhibits a hysteretic dependence between the density of the enzyme e and the production of the enzyme pe. We assume that the production of the enzyme also depends on the density of ligands and is a function of the product le. It means that production of new enzyme depends on the number of ligands degraded by the
104
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119 l 14
12
10
8
6
s 4
u
2
0
s 0
1
2
3
4
5
6
7
8
9
10 p l
Fig. 4. Nullclines oto pl ¼ 0 and oto l ¼ 0 in the absence of diffusion. The curves cross in three points. This corresponds to the existence of three steady states of the system, two of them stable and one unstable.
enzyme. Numerical simulations show that the dependence on l is necessary for pattern formation in this model. We also assume a dependence on rb, but based on our simulations we conclude that it is not as crucial. The complete model has the form: o rf ¼ lf rf brf l þ drb þ pr ðrb Þ; ot o rb ¼ lb rb þ brf l drb ; ot o 1 o2 l ll l brf l þ drb þ pl be le; l¼ ot c ox2 o d 2 o2 e le e þ p e ; e¼ ot c ox2 o pl m2 lrb pl ¼ dl þ ; 2 2 ot 1 þ pl ð1 þ rl pl bl pl Þð1 þ al rb Þ o pe m3 le pe ¼ de ; þ 2 2 ot 1 þ pe ð1 þ re pe be pe Þð1 þ ae rb Þ
ð8Þ
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
105
rb where pr(rb) = m1 or pr ðrb Þ ¼ m1 1þr . The particular form of pr is chosen for simplicity. We also b tried to model pr in a way similar to pl and pe, which makes the model more complicated and does not give qualitatively new results.
3.2. Properties of the system Using the framework of invariant rectangles (see [28]) we can show that solutions of system (8) remain positive for positive initial conditions. The time-behavior of the solutions depends on the parameters. First of all we assume that all the parameters are positive. To guarantee global existence of the solutions we must assume that rl, bl, rl and be are chosen in such way that polynomials 1 þ re p2e be pe and 1 þ rl p2l bl pl have no real roots. Otherwise, the solutions might grow indefinitely in finite time (blow-up). Accordingly, we require the following restrictions on the parameters: b2l 4rl < 0; b2e 4re < 0.
ð9Þ
We investigate only those parameter values which allow a hysteretic relation in Eqs. (6) and (7). It can be verified that for such parameter values conditions (9) are indeed fulfilled. The existence of global solutions for the parameter values satisfying (9) results from a standard argument based on the theory of bounded invariant rectangles and a priori estimates [28]. Methods outlined in [28, Chapter 14] can be used without major modifications. Since the model has a complicated structure, we investigate the behavior of solutions numerically. However, in a simplified case, a little more can be accomplished, as is explained in the following section.
4. Spatially inhomogeneous stationary solutions in a simplified case We consider a model with hysteresis but without an enzyme, i.e. o rf ¼ lf rf brf l þ drb þ pr ðrb Þ; ot o rb ¼ lb rb þ brf l drb ; ot o 1 o2 l ll l brf l þ drb þ pl ; l¼ ot c ox2 o pl m2 lrb p ¼ dl þ ot l 1 þ p2l ð1 þ rl p2l bl pl Þð1 þ al rb Þ rb with pr(rb) = m1 or pr ðrb Þ ¼ m1 1þr . b Stationary solutions of system (10)–(13) satisfy
ð10Þ ð11Þ ð12Þ ð13Þ
0 ¼ lf rf brf l þ drb þ pr ðrb Þ;
ð14Þ
0 ¼ lb rb þ brf l drb ;
ð15Þ
106
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
1 o2 l ll l brf l þ drb þ pl ; c ox2 pl m2 lrb þ 0 ¼ dl . 1 þ p2l ð1 þ rl p2l bl pl Þð1 þ al rb Þ 0¼
ð16Þ ð17Þ
If the subsystem (14)–(16) has multiple solutions, then we may find many stationary solutions of system (10)–(13) (compare [4]). By solving the first two equations (14) and (15), we reduce system (14)–(17) to the system of two equations, 1 0 ¼ Du þ f ðu; vÞ; c 0 ¼ gðu; vÞ;
ð18Þ
with u := l and v := pl and f ðu; vÞ ¼ v
bðlb þ dÞm1 u bdm1 u ll u þ ; blb u þ lf lb lf d blb u þ lf lb þ lf d
gðu; vÞ ¼ dl
v m2 m1 bu2 . þ 1 þ v2 ðbl u þ l l l dÞð1 þ rl v2 b vÞ 1 þ al bm1 u b f b f l blb uþlf lb lf d
ð19Þ
Now we assume that there exist two solutions, h0(u) and h1(u), such that g[u, hi(u)] = 0, i = 0, 1, and that there exist stationary solutions, (u0, v0) and (u1, v1), of the kinetic system, i.e. f(u, v) = 0 and g(u, v) = 0, which are stable. Stationary solutions of system (18) satisfy 1 0 ¼ Du þ f ðu; vÞ; c vðxÞ ¼ h/ðxÞ ðuðxÞÞ
ð20Þ
for / : R ! f0; 1g. We are looking for solutions of (20) satisfying v ¼ h0 ðuÞ; v ¼ h1 ðuÞ;
x 6 x; x P x.
Due to the regularity of the equations, ux is continuous [24]. Therefore, considering the zero-flux boundary conditions, ux(0) = ux(1) = 0, the function u must satisfy the following: Z u Z uð0Þ 1 u2x þ f ðs; h0 ðsÞÞds ¼ const ¼ f ðs; h0 ðsÞÞds; x < x; c 2 0 0 ð21Þ Z uð1Þ Z u 1 u2x f ðs; h1 ðsÞÞds ¼ const ¼ f ðs; h1 ðsÞÞds; x > x; þ c 2 0 0 where [0, 1] is the support of u. At x ¼ x, we obtain Z uð1Þ Z uð0Þ f ðs; h0 ðsÞÞds ¼ f ðs; h1 ðsÞÞds; u
u
where u ¼ uðxÞ.
ð22Þ
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
107
Since we are looking for solutions connecting u0 and u1 we may take u(0) = u0 and u(1) = u1. Therefore, there exists a stationary spatially inhomogeneous solution if and only if there exists u*, such that the equation, Z u0 Z u1 AðuÞ :¼ f ðs; h0 ðsÞÞds f ðs; h1 ðsÞÞds ¼ 0 ð23Þ u
u
is satisfied. Next we consider the kinetic system (19). Solutions of the algebraic equations g(u, v) = 0 and f(u, v) = 0 are presented in Fig. 5. We find that there are two positive stable stationary solutions of the kinetic system and that there are two stable branches, h0(u) and h1(u), of solutions of the equation g(u, v) = 0. We look for the value of u* such that (23) can be satisfied. Since we cannot explicitly obtain the functions hi but we know that they are monotone (see Fig. 5), we change the variables in the integrals in (23) and consider the equation, Z v1 Z v0 0 f ðH 0 ðtÞ; tÞH 0 ðtÞdt f ðH 1 ðtÞ; tÞH 01 ðtÞdt; ð24Þ AðuÞ :¼ h0 ðuÞ
v
h1 ðuÞ
8
7
6
5
v+
(u1,v1 )
4
h1(u) 3
T1
2
1
T2
h0(u) v– 0 (u0 ,v0 )
0
u–
1
2
u+
3
4
5
6
u
Fig. 5. Solutions of the algebraic equations g(u, v) = 0 and f(u, v) = 0 (u := l, v := pl) for the parameter set: b = 1.0, d = 0.9, ll = 1.2, lf = 0.1, lb = 0.01, m1 = 2, m2 = 1, al = 9, bl = 0.8, rl = 0.2 and dl = 1. (u0, v0) and (u1, v1) are stable steady states of the kinetic system, u is u-coordinate of the lower turning point (T1), u+ is u-coordinate of the upper turning point (T2), v is v-coordinate of the point of the curve g(u, v) = 0 that lies on the lower branch (h0(u)), v+ is v-coordinate of the point of the curve g(u+, v) = 0 that lies on the upper branch (h1(u)).
108
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
where Hi is the inverse function of hi, i.e., Hi hi = Identity, i = 0, 1. We explicitly find Hi using the computer algebra program Maple and calculate A for a given u. u* lies in the interval [u, u+] such that for u 2 [u, u+] there exist both h0(u) and h1(u). Thus we can calculate A(u) and A(u+) and if they are of different signs, then there exists u* such that A(u*) = 0. Numerical calculations using Maple show that A(u) = 0 exactly at u = u. Numerical simulations performed for system (10)–(13) are in agreement with the analysis. For small enough diffusion coefficients, d l ¼ 1c (which is equivalent to the domain size large enough), we observe the formation of spatially inhomogeneous solutions (fronts), which are stationary in time (see Fig. 6). The spatial profile of the solution depends on the scaling parameter c and it is
Fig. 6. The solution of system (10)–(13) for parameter set II (see Table 1) and diffusion coefficient d = 0.01. Simulations are performed using initial conditions: rf(x) = 10 for x 2 [0, 1], rb(x) = 10 for x 2 [0, 1], l(x) = 10 for x 2 [0, 1] and pl(x) = 1 for x 2 [0,0.3], pl(x) = 0.1 for x 2 (0.3,1].
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
109
Fig. 7. The solution of system (10)–(13) for different diffusion coefficients: (a) d = 0.1, (b) d = 0.0001. We observe that the profile for the bigger diffusion coefficient is smoother.
2
steeper when c is larger (compare Fig. 7(a) and (b)). We recall that c ¼ Ld l and the model is defined on the interval [0, 1] with diffusion coefficients d l ¼ 1c. This means that an increasing c corresponds to an increase in the domain size or a decrease in the diffusion coefficient of the initial problem. m1 rb leads to qualitatively the same results. Remark 4.1. Considering the function pr ¼ 1þr b
Numerical solutions of model (8) with zero-flux boundary conditions show similar behavior. They have the form of stationary or moving fronts depending on the parameters. However, for the complete system we cannot perform a similar analysis. Even assuming that one of the diffusion coefficients is equal to zero, we obtain a sub-system of ordinary differential equations which is difficult to reduce to a single equation.
5. Results of simulations for the complete model with hysteresis 5.1. Basic settings We performed simulations for model (8) with zero-flux boundary conditions for l and e. As an initial condition we chose the homogeneous state, rf = 1, rb = 1, l = 1, e = 1, pl = 1 and pe = 1 and then perturbed the density of ligands by setting l = 10 at one end of the domain on an interval of the length 1/10th or 2/10th of the domain. Simulations show that for model (8) we can have a gradient-like solution for the density of bound receptors, which is stationary in time (i.e. a standing wave, see Figs. 8 and 9(b)) or a spatio-temporal solution oscillating in time (see Fig. 9(a)). Parameters involved in the equations describing the production dynamics are chosen in such a way as to allow a hysteretic dependence of pl on l. Of course, this is only one particular choice. Furthermore, the numerical values used do not represent a systematic review of the parameter space. Rather, it is our purpose here to establish, numerically, the capabilities of the models to display the desired types of behavior. The diffusion coefficient for the ligand should be smaller than for the enzyme.
110
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
rb Fig. 8. Time evolution of rb for the parameter set I and pr ¼ m1 1þr . b
rb Fig. 9. Time evolution of rb for the parameter set II and: (a) constant pr = m1, (b) pr ¼ m1 1þr . b
Table 1 Parameters values used in simulations of the model (8) Set
b
be
d
lf
lb
ll
le
m1 = m2 = m3
al = ae
bl = be
rl = re
dl = de
I II
2.0 1.0
0.1 0.8
0.5 0.9
0.1 0.1
0.1 0.1
0.1 0.01
0.3 0.3
1.0 1.0
9.0 9.0
0.8 0.8
0.2 0.2
1.0 1.0
The parameter values used in the simulations are listed in Table 1. For parameter set I we observe the formation of a stationary pattern (standing wave) (see Fig. 8). We performed further simulations for different parameter sets, usually changing the d, b and be coefficients. For parameter set II, we observe two different behaviors depending on the form of the function pr. A constant pr results in a spatio-temporal pattern, while the Michaelis–Menten function for pr results in a stationary pattern (see Fig. 9(a) and (b)).
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
111
5.2. Simulations of cutting and grafting experiments We now explain the mathematical formulation of the models corresponding to the experiments performed on hydra in laboratory. Two basic types of such experiments, i.e. cutting and grafting experiments, involve changes to the initial state of the hydra. Thus, we consider the models studied before for normal development with different initial conditions. We formulate these conditions and explain how they are related to the experimental setting. 5.2.1. Cutting experiments Assume that v(x) is a vector of densities corresponding to the case of normal development and we cut a piece of the body column in x1 and x2. The initial conditions for the cutting experiment are u0 ðxÞ ¼ vðx1 þ xðx2 x1 ÞÞ
for x 2 ½0; 1; 0 6 x1 < x2 6 1.
ð25Þ
Since we assume that x = 0 corresponds to the upper end of the body column (head) and x = 1 to the lower end (foot), the case x1 = 0, x2 < 1 simulates a surgical removal of the foot, while the case x1 > 0, x2 = 1 – removal of the head. 2 Let us recall that for the dimensionless models the domain size equals 1 and c ¼ Ld l is a scaling coefficient depending on the domain size of the original problem. The cutting experiment changes the domain size of the problem. The ‘new’ domain is smaller. Thus, to obtain a model correspond2 ing to the cutting experiments, one has to also rescale the coefficient c. Since c ¼ Ld l , we obtain 2 cnew = c(x2 x1) . Since regeneration in hydra is morphallactic (i.e. does not involve the domain growth), we model the cutting experiments with a constant domain size (constant c). The function v(x) is determined by the values obtained from simulations of normal development (gradient-like solutions). We have numerically solved the model equations for initial conditions corresponding to various cutting experiments and found that the model with hysteresis cannot explain all cutting experiments. If we simulate cutting of the head end, we can observe reorganization. But this happens only in the case when the density of ligands is high enough for the size of production (pl) to stay in the basin of attraction of the higher steady state. Similarly if we simulate cutting of the foot, we observe the formation of the gradient-like solution if the lower end values of the initial state are close enough to the lower steady state. However, it seems to be meaningful that there is experimental evidence that after cutting a hydra in half, both parts simultaneously regenerate [18]. This means that to model this experiment we should obtain a gradient-like pattern formation for the initial conditions corresponding to both parts of the initial domain (as is the case for the models with diffusion-driven instability described in [12]). It is an open problem whether such behavior can be obtained using a model of the form (8) for any set of parameters and production kinetics. A similar observation holds for the receptor-based model proposed by Sherratt et al. [27]. Results of the simulations described by these authors only concern cutting and regeneration of the head or foot. Simulations of the cutting experiments described in [27] concern cutting of the upper or lower 1/4 of the body column. This means that new simulations are performed for the initial conditions corresponding to 3/4 of the gradient-like solution. For such initial conditions model (8) also reproduces the results of cutting experiments.
112
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
Fig. 10. The solutions of the model for initial condition corresponding to different transplantations of a head.
5.2.2. Grafting experiments Now we formulate the initial conditions describing grafting experiments. Assuming again that v(x) is a solution corresponding to the case of normal development, we put vðxÞ for x 2 ½0; x1 [ ½x2 ; 1; ð26Þ u0 ðxÞ ¼ vðx x1 þ xgraft Þ for x 2 ½x1 ; x2 . This means that, instead of the portion of tissue with coordinates [x1, x2], we transplant the portion of tissue of the same length from the position starting in xgraft. Corresponding simulations can be interpreted as predicting the outcome of grafting experiments. Fig. 10(a) and (b) show the spatio-temporal evolution of the density of bound receptors from initial conditions corresponding to the transplantation of tissue from the upper part of the body column to the lower one. We observe the formation of a second peak corresponding to a second head. The formation and persistence of the second maximum is the result of the bistability of the reaction terms. For such a model we observe various stable solutions, which interpolate between the two stable steady states.
6. Stationary homogeneous steady states In [12] we showed that a receptor-based model (1)–(4) with diffusion-driven instabilities can fully reproduce the outcome of the cutting experiments. Thus, one of the questions arising is whether it is also possible to obtain, using model (8), diffusion-driven instabilities, and how the corresponding parameter values are related to the parameter space for which we obtained the wave solutions. rb . We consider the case with the production of free receptors given by the function pr ¼ m1 1þr b We also discuss results for the constant production pr = m1.
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
113
Solving the algebraic equations for the spatially homogeneous steady state (for details see the Appendix), we obtain the two equations (33) and (34) given in the Appendix. It is interesting to find the steady states for the values of parameters used in simulations for which we obtained spatial patterns i.e. parameter values from the parameter sets I and II (see Table 1). In the case pe = 0, only one real positive value of rb satisfies Eq. (33). Now we assume that pe 5 0. We solve (33) and (34) for the parameter set II. Eq. (33) has a solution which is illustrated in Fig. 11 above using the light-color line. The solution of Eq. (34) is shown in Fig. 11 above using the dark-color line. In summary, for this set of parameters, there is only one positive real solution u1 ¼ ð0:486; 8:456; 17:383; 0:876; 13:208; 0:262Þ. There is also one real solution with u2 ¼ ð2:538; 6:041; 2:379; 0; 0:627; 0ÞÞ. For the complete system we have also the trivial pe = 0 ð solution ð u3 ¼ ð0; 0; 0; 0; 0; 0ÞÞ.
20 Eq. 28 Eq. 29
rb
15
10
5
0 0
2
4
6
8
6
8
pe 20 Eq. 28 Eq. 29
rb
15
10
5
0 0
2
4
pe
Fig. 11. Solutions of (34) crossing the solution of (33) for the parameter set II and a Michaelis–Menten production of m1 r b free receptors pr ¼ 1þr – on the picture above, and constant pr = m1 – on the picture below. We observe that there is b only one common point.
114
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119 Eq. 28 16
Eq. 29
rb
12
8
4
0 0
2
4
6
8
pe
Fig. 12. Solutions of (34) crossing the solution of (33) for the parameter set I and a production of free receptors pr = m1. We observe that there are three common points.
We perform a linear stability analysis of system (27) and find that the trivial solution, u3 , is stable, and that the solution u1 is also stable. The solution u2 is unstable. We investigate whether the steady 2state u1 can be destabilized by diffusion. Thus, we study the eigenvalues of the matrix A D lcm , where A is a linearization of the kinetic system at the state u1 , D is a matrix of diffusion coefficients and l2n;k are wavenumbers obtained from the eigenproblem for the Laplacian. We find that there exist values of the diffusion coefficients and of c for which we can select an unstable mode lm. The simulations performed for such values show that the steady state is indeed destabilized. However, the solution tends to the trivial equilibrium and no spatially inhomogeneous spatial pattern is obtained. We stress here that for the values of diffusion coefficients (and c) for which we obtained a spatially inhomogeneous pattern, the steady state u1 is always stable. We obtained qualitatively similar results for the constant production function pr = m1. For the parameter set I, there exist three positive steady states of (27) (see Fig. 12). The picture looks qualitatively the same as for the constant production of bound receptors. In both cases, two of these steady states are unstable and one is stable and remains stable also for the system with diffusion. We did not find numerically any parameter set such that there is a diffusion-driven instability in system (33), (34) leading to the formation of spatially inhomogeneous patterns.
7. Discussion – A Turing-type model versus a model with hysteresis In the framework of reaction–diffusion systems, there are essentially two ways in which a system of initially identical cells can start to differentiate: • There is a critical number of cells (related to the size of the domain), above which the spatially homogeneous attractor loses stability which leads to ‘spontaneous’ spatial patterning. This is the case for the models with a Turing-type instability.
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
115
• There is an external inducing signal or forcing term which drives the system into a new, spatially inhomogeneous state. Such a signal originates from a transplant formed by a group of already differentiated cells. In such a case the signal must be strong enough to trigger differentiation. This corresponds to a sufficiently large initial perturbation of the homogeneous steady state. This type of initialization of the pattern-forming mechanism is involved in the model with hysteresis. In [12] and in this work we approach the question of mechanisms responsible for the formation of new structures and regeneration of missing ones using the specific example of a fresh-water polyp, the hydra. We consider both mechanisms of pattern formation i.e. diffusion-driven instability and a hysteresis-driven mechanism and demonstrate their advantages and constraints in explaining the different aspects of pattern formation and regulation. The model proposed in this work and the previous Turing-type model of pattern formation in hydra [12] are based on the idea that patterns are controlled by specific cell-surface receptors, which transmit to the cells signals responsible for their differentiation. Our studies show that by taking into account the receptor reactions on the cell surface, activator–inhibitor models of pattern formation can be considerably improved in the sense that they cease to require the artificial hypothesis of the existence of two diffusible morphogens and the autocatalytic production of the molecules. In the previous Turing-type model of pattern formation in hydra [12], it was shown that an arbitrary, ‘random’ initial condition leads to a gradient-like solution. In our model organism, the hydra, such a model can explain the results of cutting experiments. For a certain range of domain size values, we observe the formation of a gradient-like pattern from any initial conditions corresponding to the distributions of molecules on the different parts of a body column of a complete organism. This mimics well the regeneration of the complete animal from an arbitrary piece of the body column which is sufficiently large. However, models with diffusion-driven instabilities (and also the activator–inhibitor models) fail to reproduce the grafting experiments. The present work shows that our model with multiple steady states and hysteretic behavior can explain the results of grafting experiments, because in such a model multiple spatially inhomogeneous solutions may coexist. The solution depends on the initial condition in a manner that is similar to how the results of the grafting experiments depend on the graft position and not on the size of the animal. We have shown that stationary fronts already arise in a simplified model which does not include any equations describing the dynamics of the enzyme or its production. For the complete model we have found numerically similar stationary fronts. However, the dynamics of the complete system is more complicated and we have also observed spatio-temporal patterns oscillating in time. Numerical simulations show that the model with hysteresis cannot reproduce the results of all types of cutting experiments. Only regeneration of a limited interval of a cut head or foot can be simulated (as in [27]). Experiments in which the animal is cut in half and both parts regenerate (compare Fig. 1) are not reproduced by the model with hysteresis (only one of the two halves of the animal is regenerated). However, such experiments are reproduced by the receptor-based system which exhibits a Turing-type instability [12]. The question arises as to whether these two kinds of experiments can be explained using the same mechanism or whether it could be a combination of previously considered mechanisms.
116
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
We checked that in model (8) diffusion-driven instabilities do not lead to the formation of a gradient-like pattern. However, the behavior of the system for different parameter sets and different kinetics describing the production of ligands and enzyme is still an open problem. As in many other cases the complexity of the biological process makes it difficult to give a rigorous derivation of the model equations. In this work we have proposed heuristic kinetics for the production of new ligands and enzyme molecules. These kinetics aim to summarize the interactions of a complicated cascade of biochemical intracellular and extracellular reactions. The main aim is to check whether a hysteretic mechanism arising from the existence of multiple steady states and switches in the dynamics is capable of explaining the observed experimental results which so far have not been able to be explained by any existing models. It would be interesting to confirm experimentally the existence of a hysteretic dependence of ligand production and their concentration in the regulatory system of hydra, in a way analogous to the work of Tyson and Novak [32] in the context of the cell cycle. This would manifest itself in rapid changes of the production rate under a continuous change of the concentration of the molecules. Further biological results are needed to gain insight into the molecular nature of cell communication and of positional information. To understand morphogenesis in hydra it is necessary to bridge the gap between experimental observations at the cellular level and at the genetic and biochemical levels. So far, about 20 ligand-activated signaling pathways (among others Wnt, TGF-b, Hedgehog and Notch) have been identified in hydra, in which at least one component is likely to play a developmental role [29]. In some cases, there is definitive experimental evidence for a developmental role as a head activator for example. In other cases, such a role is assumed based on the results obtained for other organisms. Recently, the Wnt signaling cascade has become the focus of much research attention and fragments of candidate receptor-genes were found in the cDNA library [5]. The question of how the Wnt regulatory network may be responsible for spatial patterning is the subject of our further work. We aim to build more detailed models based on the experimental results on the Wnt signaling cascades in hydra cells performed in the lab of Prof. Holstein (Zoology Institute, University of Heidelberg). Studies with heuristic kinetics suggest possible a pattern formation scenario and show that the existence of even one diffusing substance in the complicated signaling network can lead to the formation of different spatial structures. On the other hand, the constraints of the models suggest that communication over a distance may involve cell-to-cell signaling mechanisms different from simple diffusion. This is why there arises the need to explain how the intracellular and extracellular interactions are translated into macroscopic form and pattern. The complicated structures of the cell cytoplasm, of the cell nucleus and of the extracellular matrix make it evident that modeling of transport processes must be handled carefully, in many cases reflecting the porous structure of the medium [8]. This current work shows that it is justified to look for complicated nonlinear dynamics in the signaling pathways leading to synthesis of new ligands. The models proposed in this paper have shown how the nonlinearities of the signaling pathways give rise to the final spatial patterning. They may be interesting for biologists who are looking for a possible scenario of pattern formation in hydra, and also for mathematicians looking for the mechanisms involved in reaction–diffusion systems which would lead to the formation of spatially inhomogeneous solutions consistent with experiments.
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
117
Acknowledgment I would like to thank Prof. Willi Ja¨ger and Prof. Werner Mu¨ller for supporting me and Prof. Marek Kimmel for valuable comments on the manuscript. This work was supported by the grant Graduiertenkolleg 484: ‘Signalling systems and gene expression in model systems of developmental biology’. Appendix Stationary homogeneous steady states of the model (8) with a Michaelis–Menten production of rb satisfy the system, the free receptors, pr ¼ m1 1þr b rb 0 ¼ lf rf brf l þ drb þ m1 ; 1 þ rb 0 ¼ lb rb þ brf l drb ; 0¼
1 o2 l ll l brf l þ drb þ pl be le; c ox2
d 2 o2 e le e þ p e ; c ox2 pl m2 lrb þ 0 ¼ dl ; 1 þ p2l ð1 þ rl p2l bl pl Þð1 þ al rb Þ pe m3 le þ . 0 ¼ de 2 2 1 þ pe ð1 þ re pe be pe Þð1 þ ae rb Þ
ð27Þ
0¼
We solve the 1st, 2nd and 4th equation with respect to rf, l and e. Assuming that rb 5 0, we obtain p e¼ e; ð28Þ le rb ðlb þ lb rb m1 Þ ; ð29Þ rf ¼ lf ð1 þ rb Þ l¼
ðlb þ dÞlf ð1 þ rb Þ . bðlb þ lb rb m1 Þ
ð30Þ l d ðb2 þr l2 l2 b l l b Þ
For rb = 0 we can calculate that the steady states are (0, 0, 0, 0, 0, 0) and ð0; 0; e e em3ðbe 2eþll 2 l2eÞ e l e , e l e lbel ; 0; lbe le l Þ, the second of which is non-positive. Substituting (28)–(30) in (27) we reduce system (27) to the following system of three equations: ll ðlb þ dÞlf ð1 þ rb Þ be ðlb þ dÞlf ð1 þ rb Þpe ; þ drb þ bðlb þ lb rb m1 Þ bðlb þ lb rb m1 Þle m2 rb ðlb þ dÞlf ð1 þ rb Þ dl p l 0¼ ; 2 ð1 þ pl Þ bðlb þ lb rb m1 Þð1 þ al rb Þð1 þ rl p2l bl pl Þ m3 ðlb þ dÞlf ð1 þ rb Þpe de p e . 0¼ 2 ð1 þ pe Þ bðlb þ lb rb m1 Þle ð1 þ ae rb Þð1 þ re p2e be pe Þ
0 ¼ pl rb ðlb þ dÞ þ
ð31Þ
118
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
We solve the first of the above equations with respect to pl and obtain pl ¼ ðbe lf pe lb be lf pe d be lf pe lb rb be lf pe drb þ rb ble l2b þ r2b ble l2b rb ble lb m1 ll lf le lb ll lf le d ll lf le lb rb ll lf le drb Þ=bðlb þ lb rb m1 Þle .
ð32Þ
Finally we obtain two equations 0 ¼ dl ðbe lf pe lb belf pe d be lf pe lb rb belf pe drb þ rb ble l2b þ r2b ble l2b rb ble lb m1 ll lf le lb ll lf le d ll lf le lb rb ll lf le drb Þ=b=ðlb þ lb rb m1 Þ=le =ð1 þ ðbe lf pe lb belf pe d belf pe lb rb belf pe drb þ rb ble l2b 2
þ r2b ble l2b rb ble lb m1 ll lf le lb ll lf le d ll lf le lb rb ll lf le drb Þ =b2 =ðlb þ lb rb m1 Þ2 =l2e Þ m2 rb ðlb þ dÞlf ð1 þ rb Þ=b=ðlb þ lb rb m1 Þ =ð1 þ al rb Þ=ð1 þ rl ðbe lf pe lb be lf pe d be lf pe lb rb be lf pe drb þ rb ble l2b þ r2b ble l2b rb ble lb m1 ll lf le lb ll lf le d ll lf le lb rb 2
2
ll lf le drb Þ =b2 =ðlb þ lb rb m1 Þ =l2e bl ðbelf pe lb be lf pe d be lf pe lb rb be lf pe drb þ rb ble l2b þ r2b ble l2b rb ble lb m1 ll lf le lb ll lf le d ll lf le lb rb ll lf le drb Þ=b=ðlb þ lb rb m1 Þ=le Þ; 0¼
m3 ðlb þ dÞlf ð1 þ rb Þpe de pe þ . 2 ð1 þ pe Þ bðlb lb rb þ m1 Þle ð1 þ ae rb Þð1 þ re p2e be pe Þ
ð33Þ ð34Þ
References [1] A. Babloyantz, M. Sanglier, Chemical instabilities of ‘‘All-or-None’’ type in b-galactosidase induction and active transport, FEBS Lett. 23 (1972) 346. [2] F.R. Cross, V. Archambault, M. Miller, M. Klovstad, Testing a mathematical model of the yeast cell cycle, Mol. Biol. Cell 13 (2002) 52. [3] A. Gierer, H. Meinhardt, A theory of biological pattern formation, Kybernetik 12 (1972) 30. [4] S. Heinze, B. Schweizer, H. Schwetlick, Existence of front solutions in degenerate reaction diffusion systems, SFB preprint, University of Heidelberg, 2004. [5] B. Hobmayer, F. Rentzsch, K. Kuhn, C.M. Happel, C.C. Laue, P. Snyder, U. Rothba¨cher, T.W. Holstein, Wnt signaling and axis formation in the diploblastic metazoan Hydra, Nature 407 (2000) 186. [6] F. Hoppensteadt, W. Ja¨ger, Pattern formation by bacteria, in: S. Levin (Ed.), Lecture Notes in Biomathematics: Biological Growth and Spread, Springer, Heidelberg, 1980. [7] F. Hoppensteadt, W. Ja¨ger, C. Po¨ppe, A hysteresis model for bacterial growth patterns, in: S. Levin (Ed.), Lecture Notes in Biomathematics: Modelling of Patterns in Space and Time, Springer, Heidelberg, 1983. [8] U. Hornung, W. Ja¨ger, A. Mikelic, Reactive transport through an array of cells with semi-permeable membranes, Mathematical Modeling and Numerical Analysis 28 (1994) 59. [9] J. Keener, J. Sneyd, Mathematical Physiology, Springer, New York, 1998. [10] Ch.T. Klein, Hysteresis-driven structure formation in biochemical networks, J. Theor. Biol. 194 (1998) 263. [11] D.A. Lauffenburger, J.J. Linderman, Receptors. Models for Binding, Traficking, and Signaling, Oxford University, New York, 1993.
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
119
[12] A. Marciniak-Czochra, Receptor-based models with diffusion-driven instability for pattern formation in hydra, J. Biol. Syst. 11 (2003) 293. [13] J.-L. Martiel, A. Goldbeter, A model based on receptor desensitization for cyclic AMP signaling in Disctyostelium cells, Biophys. J. 57 (1987) 807. [14] H. Meinhardt, A model for pattern formation of hypostome, tentacles and foot in hydra: how to form structures close to each other, how to form them at a distance, Dev. Biol. 157 (1993) 321. [15] P.B. Monk, H.G. Othmer, Cyclic AMP oscillations in suspensions of Disctyostelium Discoideum, Philos. Trans. R. Soc. Lond. B323 (1989) 185. [16] J.D. Murray, Mathematical Biology, Springer, Berlin, 1993. [17] W.A. Mu¨ller, Ectopic head and foot formation in hydra. Diacylglycerol induced increase in potential value and assistance of the head in foot formation, Differentiation 42 (1990) 131. [18] W.A. Mu¨ller, Pattern control in hydra: basic experiments and concepts, in: H.G. Othmer, P.K. Maini, J.D. Murray (Eds.), Experimental and Theoretical Advances in Biological Pattern Formation, Plenum, New York, 1993. [19] W.A. Mu¨ller, Competition of factors and cellular resources as a principle of pattern formation in hydra, Dev. Biol. 167 (1995) 175. [20] W.A. Mu¨ller, Developmental Biology, Springer, New York, 1997. [21] S.C. Mu¨ller, G. Venzl, Pattern formation in precipitation processes, in: S. Levin (Ed.), Lecture Notes in Biomathematics: Modelling of Patterns in Space and Time, Springer, Heidelberg, 1983. [22] B. Novak, J. Tyson, Numerical analysis of a comprehensive model of M-phase control in Xenopus oocyte extracts and intact embryos, J. Cell. Sci. 106 (1993) 1153. [23] B. Novak, Z. Pataki, Mathematical model of cell division cycle of fission yeast, Chaos 11 (2001) 277. [24] F. Rothe, Global Solutions of Reaction–Diffusion Systems, Springer, Berlin, 1984. [25] F.F. Seelig, B. Denzle, Hysteresis without autocatalysis in enzyme systems with substrate inhibition, FEBS Lett. 24 (1972) 283. [26] W. Sha, J. Moore, K. Chen, A. Lassaletta, Y. Chung-Seon, J. Tyson, Hysteresis drives cell-cycle transitions in Xenopus Laevis egg extracts, PNAS 100 (2003) 975. [27] J.A. Sherratt, P.K. Maini, W. Ja¨ger, W. Mu¨ller, A receptor-based model for pattern formation in hydra, Forma 10 (1995) 77. [28] J. Smoller, Shock-Waves and Reaction–Diffusion Equations, Springer, New York, 1994. [29] R. Steele, Developmental signaling in hydra: what does it take to build a ‘‘simple’’ animal? Developmental Biology 248 (2002) 199. [30] U. Technau, B. Hobmayer, F. Rentzsch, T. Holstein, Molecular and cellular analysis of de novo pattern formation in Hydra, in: A. Deutsch (Ed.), Function and Regulation of Cellular Systems: Experiments and Models, Birkha¨user-Verlag, Basel, 2002. [31] C.D. Thron, Bistable biochemical switching and the control of the events of the cell cycle, Oncogene 15 (1997) 317. [32] J.J. Tyson, B. Novak, Regulation of the eukaryotic cell-cycle: molecular antagonism, hysteresis, and irreversible transitions, J. Theor. Biol. 210 (2001) 249. [33] L. Wolpert, Positional information and the spatial pattern of cellular differentiation, J. Theor. Biol. 25 (1969) 1.
Receptor-based models with hysteresis for pattern formation in hydra Anna Marciniak-Czochra Institute of Applied Mathematics, University of Heidelberg, Im Neuenheimer Feld 294, 69120 Heidelberg, Germany Received 10 August 2005; received in revised form 27 August 2005; accepted 5 October 2005
Abstract In this paper, we propose a new receptor-based model for pattern formation and regulation in a freshwater polyp, namely hydra. The model is defined in the form of a system of reaction–diffusion equations with zero-flux boundary conditions coupled with a system of ordinary differential equations. The production of diffusible biochemical molecules has a hysteretic dependence on the density of these molecules and is modeled by additional ordinary differential equations. We study the hysteresis-driven mechanism of pattern formation and we demonstrate the advantages and constraints of its ability to explain different aspects of pattern formation and regulation in hydra. The properties of the model demonstrate a range of stationary and oscillatory spatially heterogeneous patterns, arising from multiple spatially homogeneous steady states and switches in the production rates. 2005 Elsevier Inc. All rights reserved. Keywords: Pattern formation; Hydra; Receptor-based model; Hysteresis
1. Introduction One of the crucial issues in developmental biology is to understand how coordinated systems of positional information are established during an organism’s development and how cells in the organism respond to the associated signaling cues, processes which ultimately result in the subdivided and patterned tissues of multicellular organisms [20,30,33]. E-mail address: [email protected] 0025-5564/$ - see front matter 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.mbs.2005.10.004
98
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
Experiments suggest that during development cells respond to local positional cues that are dynamically regulated. The hypothesis is that cells differentiate according to positional information [33]. The question is how this information is supplied to the cells. There exists a number of models for pattern formation and regulation based on the idea that positional information is supplied to cells by a diffusible biochemical morphogen [3,16]. Both regulatory and signaling molecules (ligands) act by binding and activating receptor molecules which are located in the cell membrane (or, with lipophilic ligands, in the cytoplasm) [20,18,11]. Thus in this paper we hypothesize that the positional value of the cell can be determined by the density of bound receptors which do not diffuse [18,19]. In this paper the test organism for mathematical modeling is a fresh-water polyp, the hydra, a simple organism which can be treated as a model for axis formation in higher organisms [20,30]. Hydra is an example of an organism in which the development of theoretical models has had strong influence on experimental design, starting with the positional information ideas of Wolpert [33] and the activator–inhibitor model of Gierer and Meinhardt [3,14]. Experiments performed on hydra provide a good ground for testing the abilities and limitations of such mathematical models.
2. A receptor-based model of pattern formation in hydra Hydra, a fresh-water polyp, has a tubular body about 5 mm long with a whorl of tentacles surrounding the mouth at the upper end and a disk-shaped organ for adhesion at the lower end. The longitudinal pattern is subdivided into a head, a gastric region, a budding zone (where new animals are generated by a process of natural cloning), a stalk and a foot. Experiments performed on hydra suggest that the function of cells is determined by their location. Following Mu¨ller [18], we assume that the positional information is determined by the density of bound receptors. Our models are based on the idea that epithelial cells secrete ligands (regulatory biochemical) which diffuse locally within the interstitial space (the space between two layers of cells lining the inner and outer body surface respectively) and bind to free receptors on the cell surface. A number of different receptor-based models for the evolution of spatial patterns have already been proposed, for example for the formation of cAMP patterns in Dictyostelium slime mold cells in [13] and [15]. The models considered here in this paper are simple compared with these aforementioned rather detailed and complicated schemes. This is due to the absence of a detailed biochemical knowledge of pattern control in hydra. Using simple and general models we try to clarify the ideas of the receptor-based mechanism of pattern control and explore their implications. The models we study are motivated by the model of Sherratt et al. and its receptor-based approach [27]. The objective is to build a receptor-based model without imposing an initial gradient as was done in the model of Sherratt et al. We investigate what kind of mechanisms regulating the dynamics of cell surface receptors and diffusing biochemical molecules can explain the observed results of experiments. We validate our models with the results of two crucial experiments: • Cutting experiments: After a cut, cells of the gastric region of the animal can reconstitute the missing structure, i.e. the head or foot, according to their position along the body axis. In hydra, regeneration occurs mainly by the re-patterning of existing tissues and is an example
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
99
of morphallaxis [33]. The lack of growth requirements for regeneration is shown in heavily irradiated hydra. No cell divisions occur in these animals, but they can still regenerate normally. Thus the mechanism of pattern formation in hydra seems to be independent of growth. Overlapping cut levels show that the same cells can form either the gastric region, or the head, or the foot, according to their position along the body axis (see Fig. 1). The most significant is the experiment which shows that after a transverse cut both parts of the animal can regenerate [18]. • Grafting experiments: Pieces of tissue are grafted from one animal to another and the outcome of such experiments depends on the change of the position along the body axes. Transplantation of the tissue from parts of the body column close to the head induces head formation while transplantation of the tissue from the position close to the foot results in foot formation [17] (see Fig. 2). In [12] we proposed a model describing receptor-ligand binding. The model takes into consideration the density of free receptors, of bound receptors and of ligands, as well as of a second diffusible substance functioning as an enzyme. We assume that new ligands and new free receptors are produced on the cell surface via a combination of recycling (dissociation of bound receptors) and de novo production within the cell. The ligands then bind reversibly to free receptors and bound receptors are internalized into the cell. Bound receptors also may dissociate. Both ligands and free receptors undergo natural decay. In addition, cells secrete an enzyme which diffuses along the body column and degrades the ligands. We consider a one-dimensional epithelial sheet of length L and denote the concentration of free receptors, bound receptors, ligands and enzyme densities by rf, rb, l and e respectively. For
Fig. 1. Cutting experiment. The hydra regenerates after a transverse cut of cells of the gastric region (from both upper and lower half of the body column). The experiment shows the availability of positional information cues. A hydra, as shown in the middle, is cut in two ways. In one experiment (left hand-side) the lower body column is removed, in the second experiment (right hand-side) the upper part is removed. The cut levels are not identical but somewhat different to show that one and the same group of cells (marked in grey) can form a foot (left hand-side), or a head (right handside) or a gastric segment (original state in the middle). The function of the cells depends on the position along the body column. [Courtesy of W. Mu¨ller].
100
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
Fig. 2. Grafting experiment. Determination of relative positional information values by transplantation. Pieces of tissue are grafted from one animal to another and one of three outcomes is observed. (1) If the tissue is transplanted from the upper position along the body column to the lower position then a new head is formed. (2) If the former and new position is the same then the piece is integrated and nothing is observed. (3) If the tissue is grafted to the upper position a new foot is formed. We can see that the disparities between the positional value of the transplant and the surrounding host tissue result in the head formation or the foot formation respectively. [Courtesy of W. Mu¨ller].
simplicity we assume that all binding processes are governed by the law of mass action without saturation effects. The model has the following form: o ð1Þ rf ¼ lf rf þ pr ðrb Þ brf l þ drb ; ot o ð2Þ rb ¼ lb rb þ brf l drb ; ot o o2 ð3Þ l ¼ d l 2 l ll l brf l þ pl ðrb Þ þ drb be le; ot ox o o2 ð4Þ e ¼ d e 2 e le e þ pe ðl; rb Þ ot ox with zero-flux boundary conditions for l and e. dl and de are diffusion coefficients for the ligands 2 and enzyme. By rescaling the system, we set dl to 1c and de to dc2 , where c ¼ Ld l and d 2 ¼ dd el and consider the model on the domain [0, 1]. The rates of decay of free receptors, bound receptors,
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
101
ligands and enzyme are denoted by lf, lb, ll and le respectively. d is the rate of dissociation, b is the rate of binding of ligands and free receptors and be is the rate of binding of ligands and the enzyme. By pr, pl and pe we denote production rates of new free receptors, ligands and enzyme. In [12], we considered different nonlinear production terms of free receptors and diffusing biological m1 rb , substances, using as an example Michaelis–Menten kinetics, which resulted in pr ðrb Þ ¼ 1þr b m3 rb m2 r b pl ðrb Þ ¼ 1þrb and pe ðrb Þ ¼ 1þrb . We considered also the possibility of production of the enzyme which would depend on the density of ligands (it could be regulated via some other receptors m3 l m1 rb not considered in our model). This resulted in the kinetic terms pr ¼ 1þr , pl = m2rb2 and pe ¼ 1þl . b In [12] we showed that a diffusion-driven (Turing-type) mechanism for pattern formation could be involved in a receptor-based model. This uses the special features of diffusion which results in the destabilization of the spatially homogeneous steady state and the formation of spatially heterogeneous patterns. In models involving this mechanism patterns can arise spontaneously (see the survey by Murray [16]). We note that the above system with diffusion-driven instability is able to simulate the cutting experiments. However, the model fails to reproduce the grafting experiments (compare Section 5.2). The reason for this is that the final heterogeneous pattern is stable to spatial perturbations such as those arising in models of transplantation. Transplantation of the tissue from one part of the body column of a host to another part of the body column of a donor corresponds to the changing of values of all the variables in that part of the domain. New values create a local maximum or local minimum depending on the relationship between the host and donor positions of the transplant. In the models with diffusion-driven instability the gradient pattern is stable for a given range of diffusion coefficients and a local perturbation of this state cannot lead to the creation of another type of spatial structure observed in the grafting experiments. Simulations of the model show, for example, that if two peaks are present in the initial conditions (which corresponds to grafting a new head to the lower part of the body column), the second peak will disappear after a short time and a gradient pattern will appear. Stable solutions with 2 peaks can only be obtained by increasing the size of the domain at the same time. This deficiency of models with the Turing-type mechanism of pattern formation encourages us to look for an alternative mechanism which could also explain the results of grafting experiments. 3. Modeling hysteresis As stated above, the results of grafting experiments on hydra suggest that the model should involve a memory-based relation. Thus far we have modeled the functions pl and pe as nonlinear functions depending on rb (and l in the case of pe) and then determined which nonlinearities the system allows a Turing-type pattern formation. It seems appropriate to consider pl and pe to be determined by a system of ordinary differential equations, which is a simplified model of the metabolic processes of the cells. As in many other cases the complexity of the biological process makes it difficult to give a rigorous derivation of these equations. We can only use heuristic reaction terms to depict and investigate the patterning process. Existing qualitative theory of reaction–diffusion systems is mainly focused on stability conditions for homogeneous steady states, see, for example [16,28]. It is very important to study the dynamics of systems functioning far from such equilibria. Here, we consider the mechanism of pattern formation which results from the existence of multiple steady states. Switching between
102
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
different states arises in many biological processes [9,26]. One such mechanism is hysteresis, which seems to be important in the modeling of biological development, since according to the experimental observations, inductive signals are present only during a certain time interval of development. The signal triggers changes in the cell’s nucleus and evokes differentiation which does not revert when the signal is stopped. Hysteresis results from the existence of multiple steady states (see Fig. 3). Suppose we have three steady states, two outer states which are stable, and one inner state which is unstable. The system may jump from the first stable state to the second when certain parameters are changed. If these parameters are changed back to their previous values, the system remains in the second state until the parameters reach another critical value. It is known that hysteresis kinetics may lead to spatial patterns in chemical and biological systems, such as Liesegang rings formed by precipitating colloids [21], bacterial growth patterns [6,7] and DNA melting/reannealing. Models based on a hysteresis-type bistability were also used in chemistry to describe enzymatic reactions [1,25,10]. Prominent and well-documented examples of hysteresis-driven mechanisms in biology are cell-cycle transitions [31], for example in Xenopus laevis egg extracts [26,22] and cell division cycles of fission yeast [23]. A mathematical model of cell-cycle progression predicts that irreverspl
5
4.5
C
4
stable 3.5
3
2.5
D
2
unstable
1.5
1
0.5
0
0
B
stable
A 0.5
1
l1
1.5
2
2.5
3
3.5
4
l
l2
Fig. 3. Nullclines of gl – hysteretic dependence of pl on l. The two outer branches are stable, the inner branch is unstable. If the system evolves along the stable branch AB, the density of ligands l increases and it reaches the value l2 where the system becomes unstable and it jumps to the other stable branch CD. Decreasing the value of l below l2 does not cause any change and the system remains on the upper stable branch till the value of l approaches l1 from above. The l value needs to drop below l1 to shift the system back to the lower stable branch AB.
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
103
ible transitions into and out of mitosis are driven by hysteresis in a molecular control system [26,32]. The first experimental demonstration of hysteresis in the eucaryotic cell engine has been reported recently [2]. 3.1. Additional equations modeling the production terms We propose to modify the model (1)–(4) and consider the production of diffusing substances described by new variables pl and pe, determined by additional equations, o p ¼ gl ðrb ; l; pl Þ; ot l ð5Þ o p ¼ ge ðrb ; l; e; pe Þ. ot e Here all the variables depend both on time t and space x, that is pl = pl(t, x) and so forth. However the evolution in time of pl and pe is prescribed in each point x by the same rule, given by the kinetics gl and ge. In this sense we refer to system (5) as a system of ordinary differential equations. Spatial changes in the dynamics of pl and pe are caused by coupling these dynamics with the dynamics of the ligands and enzyme, which diffuse spatially. Biological signal transduction mechanisms are complicated cascades of different biochemical reactions which are not, as yet, understood in complete detail. Therefore in this paper, we consider particular reaction kinetics which allow the testing of the possibilities and restrictions of our model (1)–(5). Such functions should be biologically realistic and fulfill some basic assumptions. We assume that: • Production of ligands and enzyme in a steady state has a hysteretic dependence on the amount of ligands and the enzyme present (see Fig. 4). • Bound receptors trigger the production of ligands and enzyme. • Production is controlled by its own size. To model the hysteretic dependence of pl on l and pe on e we take the following kinetic functions: pl m2 lrb þ ; 2 2 1 þ pl ð1 þ rl pl bl pl Þð1 þ al rb Þ pe m3 le . þ ge ¼ de 2 2 1 þ pe ð1 þ re pe be pe Þð1 þ ae rb Þ gl ¼ dl
ð6Þ ð7Þ
Fig. 4 shows the nullclines of gl crossing the nullclines of the kinetic system of Eq. (3). We show this in the coordinate system l, pl. The plot is drawn for rb = 1, but it remains qualitatively the same for all rb. We have chosen the parameters to obtain three solutions of the subsystem o p ¼ 0, oto l ¼ 0. ot l To model the production of the enzyme we use a similar function. It exhibits a hysteretic dependence between the density of the enzyme e and the production of the enzyme pe. We assume that the production of the enzyme also depends on the density of ligands and is a function of the product le. It means that production of new enzyme depends on the number of ligands degraded by the
104
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119 l 14
12
10
8
6
s 4
u
2
0
s 0
1
2
3
4
5
6
7
8
9
10 p l
Fig. 4. Nullclines oto pl ¼ 0 and oto l ¼ 0 in the absence of diffusion. The curves cross in three points. This corresponds to the existence of three steady states of the system, two of them stable and one unstable.
enzyme. Numerical simulations show that the dependence on l is necessary for pattern formation in this model. We also assume a dependence on rb, but based on our simulations we conclude that it is not as crucial. The complete model has the form: o rf ¼ lf rf brf l þ drb þ pr ðrb Þ; ot o rb ¼ lb rb þ brf l drb ; ot o 1 o2 l ll l brf l þ drb þ pl be le; l¼ ot c ox2 o d 2 o2 e le e þ p e ; e¼ ot c ox2 o pl m2 lrb pl ¼ dl þ ; 2 2 ot 1 þ pl ð1 þ rl pl bl pl Þð1 þ al rb Þ o pe m3 le pe ¼ de ; þ 2 2 ot 1 þ pe ð1 þ re pe be pe Þð1 þ ae rb Þ
ð8Þ
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
105
rb where pr(rb) = m1 or pr ðrb Þ ¼ m1 1þr . The particular form of pr is chosen for simplicity. We also b tried to model pr in a way similar to pl and pe, which makes the model more complicated and does not give qualitatively new results.
3.2. Properties of the system Using the framework of invariant rectangles (see [28]) we can show that solutions of system (8) remain positive for positive initial conditions. The time-behavior of the solutions depends on the parameters. First of all we assume that all the parameters are positive. To guarantee global existence of the solutions we must assume that rl, bl, rl and be are chosen in such way that polynomials 1 þ re p2e be pe and 1 þ rl p2l bl pl have no real roots. Otherwise, the solutions might grow indefinitely in finite time (blow-up). Accordingly, we require the following restrictions on the parameters: b2l 4rl < 0; b2e 4re < 0.
ð9Þ
We investigate only those parameter values which allow a hysteretic relation in Eqs. (6) and (7). It can be verified that for such parameter values conditions (9) are indeed fulfilled. The existence of global solutions for the parameter values satisfying (9) results from a standard argument based on the theory of bounded invariant rectangles and a priori estimates [28]. Methods outlined in [28, Chapter 14] can be used without major modifications. Since the model has a complicated structure, we investigate the behavior of solutions numerically. However, in a simplified case, a little more can be accomplished, as is explained in the following section.
4. Spatially inhomogeneous stationary solutions in a simplified case We consider a model with hysteresis but without an enzyme, i.e. o rf ¼ lf rf brf l þ drb þ pr ðrb Þ; ot o rb ¼ lb rb þ brf l drb ; ot o 1 o2 l ll l brf l þ drb þ pl ; l¼ ot c ox2 o pl m2 lrb p ¼ dl þ ot l 1 þ p2l ð1 þ rl p2l bl pl Þð1 þ al rb Þ rb with pr(rb) = m1 or pr ðrb Þ ¼ m1 1þr . b Stationary solutions of system (10)–(13) satisfy
ð10Þ ð11Þ ð12Þ ð13Þ
0 ¼ lf rf brf l þ drb þ pr ðrb Þ;
ð14Þ
0 ¼ lb rb þ brf l drb ;
ð15Þ
106
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
1 o2 l ll l brf l þ drb þ pl ; c ox2 pl m2 lrb þ 0 ¼ dl . 1 þ p2l ð1 þ rl p2l bl pl Þð1 þ al rb Þ 0¼
ð16Þ ð17Þ
If the subsystem (14)–(16) has multiple solutions, then we may find many stationary solutions of system (10)–(13) (compare [4]). By solving the first two equations (14) and (15), we reduce system (14)–(17) to the system of two equations, 1 0 ¼ Du þ f ðu; vÞ; c 0 ¼ gðu; vÞ;
ð18Þ
with u := l and v := pl and f ðu; vÞ ¼ v
bðlb þ dÞm1 u bdm1 u ll u þ ; blb u þ lf lb lf d blb u þ lf lb þ lf d
gðu; vÞ ¼ dl
v m2 m1 bu2 . þ 1 þ v2 ðbl u þ l l l dÞð1 þ rl v2 b vÞ 1 þ al bm1 u b f b f l blb uþlf lb lf d
ð19Þ
Now we assume that there exist two solutions, h0(u) and h1(u), such that g[u, hi(u)] = 0, i = 0, 1, and that there exist stationary solutions, (u0, v0) and (u1, v1), of the kinetic system, i.e. f(u, v) = 0 and g(u, v) = 0, which are stable. Stationary solutions of system (18) satisfy 1 0 ¼ Du þ f ðu; vÞ; c vðxÞ ¼ h/ðxÞ ðuðxÞÞ
ð20Þ
for / : R ! f0; 1g. We are looking for solutions of (20) satisfying v ¼ h0 ðuÞ; v ¼ h1 ðuÞ;
x 6 x; x P x.
Due to the regularity of the equations, ux is continuous [24]. Therefore, considering the zero-flux boundary conditions, ux(0) = ux(1) = 0, the function u must satisfy the following: Z u Z uð0Þ 1 u2x þ f ðs; h0 ðsÞÞds ¼ const ¼ f ðs; h0 ðsÞÞds; x < x; c 2 0 0 ð21Þ Z uð1Þ Z u 1 u2x f ðs; h1 ðsÞÞds ¼ const ¼ f ðs; h1 ðsÞÞds; x > x; þ c 2 0 0 where [0, 1] is the support of u. At x ¼ x, we obtain Z uð1Þ Z uð0Þ f ðs; h0 ðsÞÞds ¼ f ðs; h1 ðsÞÞds; u
u
where u ¼ uðxÞ.
ð22Þ
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
107
Since we are looking for solutions connecting u0 and u1 we may take u(0) = u0 and u(1) = u1. Therefore, there exists a stationary spatially inhomogeneous solution if and only if there exists u*, such that the equation, Z u0 Z u1 AðuÞ :¼ f ðs; h0 ðsÞÞds f ðs; h1 ðsÞÞds ¼ 0 ð23Þ u
u
is satisfied. Next we consider the kinetic system (19). Solutions of the algebraic equations g(u, v) = 0 and f(u, v) = 0 are presented in Fig. 5. We find that there are two positive stable stationary solutions of the kinetic system and that there are two stable branches, h0(u) and h1(u), of solutions of the equation g(u, v) = 0. We look for the value of u* such that (23) can be satisfied. Since we cannot explicitly obtain the functions hi but we know that they are monotone (see Fig. 5), we change the variables in the integrals in (23) and consider the equation, Z v1 Z v0 0 f ðH 0 ðtÞ; tÞH 0 ðtÞdt f ðH 1 ðtÞ; tÞH 01 ðtÞdt; ð24Þ AðuÞ :¼ h0 ðuÞ
v
h1 ðuÞ
8
7
6
5
v+
(u1,v1 )
4
h1(u) 3
T1
2
1
T2
h0(u) v– 0 (u0 ,v0 )
0
u–
1
2
u+
3
4
5
6
u
Fig. 5. Solutions of the algebraic equations g(u, v) = 0 and f(u, v) = 0 (u := l, v := pl) for the parameter set: b = 1.0, d = 0.9, ll = 1.2, lf = 0.1, lb = 0.01, m1 = 2, m2 = 1, al = 9, bl = 0.8, rl = 0.2 and dl = 1. (u0, v0) and (u1, v1) are stable steady states of the kinetic system, u is u-coordinate of the lower turning point (T1), u+ is u-coordinate of the upper turning point (T2), v is v-coordinate of the point of the curve g(u, v) = 0 that lies on the lower branch (h0(u)), v+ is v-coordinate of the point of the curve g(u+, v) = 0 that lies on the upper branch (h1(u)).
108
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
where Hi is the inverse function of hi, i.e., Hi hi = Identity, i = 0, 1. We explicitly find Hi using the computer algebra program Maple and calculate A for a given u. u* lies in the interval [u, u+] such that for u 2 [u, u+] there exist both h0(u) and h1(u). Thus we can calculate A(u) and A(u+) and if they are of different signs, then there exists u* such that A(u*) = 0. Numerical calculations using Maple show that A(u) = 0 exactly at u = u. Numerical simulations performed for system (10)–(13) are in agreement with the analysis. For small enough diffusion coefficients, d l ¼ 1c (which is equivalent to the domain size large enough), we observe the formation of spatially inhomogeneous solutions (fronts), which are stationary in time (see Fig. 6). The spatial profile of the solution depends on the scaling parameter c and it is
Fig. 6. The solution of system (10)–(13) for parameter set II (see Table 1) and diffusion coefficient d = 0.01. Simulations are performed using initial conditions: rf(x) = 10 for x 2 [0, 1], rb(x) = 10 for x 2 [0, 1], l(x) = 10 for x 2 [0, 1] and pl(x) = 1 for x 2 [0,0.3], pl(x) = 0.1 for x 2 (0.3,1].
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
109
Fig. 7. The solution of system (10)–(13) for different diffusion coefficients: (a) d = 0.1, (b) d = 0.0001. We observe that the profile for the bigger diffusion coefficient is smoother.
2
steeper when c is larger (compare Fig. 7(a) and (b)). We recall that c ¼ Ld l and the model is defined on the interval [0, 1] with diffusion coefficients d l ¼ 1c. This means that an increasing c corresponds to an increase in the domain size or a decrease in the diffusion coefficient of the initial problem. m1 rb leads to qualitatively the same results. Remark 4.1. Considering the function pr ¼ 1þr b
Numerical solutions of model (8) with zero-flux boundary conditions show similar behavior. They have the form of stationary or moving fronts depending on the parameters. However, for the complete system we cannot perform a similar analysis. Even assuming that one of the diffusion coefficients is equal to zero, we obtain a sub-system of ordinary differential equations which is difficult to reduce to a single equation.
5. Results of simulations for the complete model with hysteresis 5.1. Basic settings We performed simulations for model (8) with zero-flux boundary conditions for l and e. As an initial condition we chose the homogeneous state, rf = 1, rb = 1, l = 1, e = 1, pl = 1 and pe = 1 and then perturbed the density of ligands by setting l = 10 at one end of the domain on an interval of the length 1/10th or 2/10th of the domain. Simulations show that for model (8) we can have a gradient-like solution for the density of bound receptors, which is stationary in time (i.e. a standing wave, see Figs. 8 and 9(b)) or a spatio-temporal solution oscillating in time (see Fig. 9(a)). Parameters involved in the equations describing the production dynamics are chosen in such a way as to allow a hysteretic dependence of pl on l. Of course, this is only one particular choice. Furthermore, the numerical values used do not represent a systematic review of the parameter space. Rather, it is our purpose here to establish, numerically, the capabilities of the models to display the desired types of behavior. The diffusion coefficient for the ligand should be smaller than for the enzyme.
110
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
rb Fig. 8. Time evolution of rb for the parameter set I and pr ¼ m1 1þr . b
rb Fig. 9. Time evolution of rb for the parameter set II and: (a) constant pr = m1, (b) pr ¼ m1 1þr . b
Table 1 Parameters values used in simulations of the model (8) Set
b
be
d
lf
lb
ll
le
m1 = m2 = m3
al = ae
bl = be
rl = re
dl = de
I II
2.0 1.0
0.1 0.8
0.5 0.9
0.1 0.1
0.1 0.1
0.1 0.01
0.3 0.3
1.0 1.0
9.0 9.0
0.8 0.8
0.2 0.2
1.0 1.0
The parameter values used in the simulations are listed in Table 1. For parameter set I we observe the formation of a stationary pattern (standing wave) (see Fig. 8). We performed further simulations for different parameter sets, usually changing the d, b and be coefficients. For parameter set II, we observe two different behaviors depending on the form of the function pr. A constant pr results in a spatio-temporal pattern, while the Michaelis–Menten function for pr results in a stationary pattern (see Fig. 9(a) and (b)).
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
111
5.2. Simulations of cutting and grafting experiments We now explain the mathematical formulation of the models corresponding to the experiments performed on hydra in laboratory. Two basic types of such experiments, i.e. cutting and grafting experiments, involve changes to the initial state of the hydra. Thus, we consider the models studied before for normal development with different initial conditions. We formulate these conditions and explain how they are related to the experimental setting. 5.2.1. Cutting experiments Assume that v(x) is a vector of densities corresponding to the case of normal development and we cut a piece of the body column in x1 and x2. The initial conditions for the cutting experiment are u0 ðxÞ ¼ vðx1 þ xðx2 x1 ÞÞ
for x 2 ½0; 1; 0 6 x1 < x2 6 1.
ð25Þ
Since we assume that x = 0 corresponds to the upper end of the body column (head) and x = 1 to the lower end (foot), the case x1 = 0, x2 < 1 simulates a surgical removal of the foot, while the case x1 > 0, x2 = 1 – removal of the head. 2 Let us recall that for the dimensionless models the domain size equals 1 and c ¼ Ld l is a scaling coefficient depending on the domain size of the original problem. The cutting experiment changes the domain size of the problem. The ‘new’ domain is smaller. Thus, to obtain a model correspond2 ing to the cutting experiments, one has to also rescale the coefficient c. Since c ¼ Ld l , we obtain 2 cnew = c(x2 x1) . Since regeneration in hydra is morphallactic (i.e. does not involve the domain growth), we model the cutting experiments with a constant domain size (constant c). The function v(x) is determined by the values obtained from simulations of normal development (gradient-like solutions). We have numerically solved the model equations for initial conditions corresponding to various cutting experiments and found that the model with hysteresis cannot explain all cutting experiments. If we simulate cutting of the head end, we can observe reorganization. But this happens only in the case when the density of ligands is high enough for the size of production (pl) to stay in the basin of attraction of the higher steady state. Similarly if we simulate cutting of the foot, we observe the formation of the gradient-like solution if the lower end values of the initial state are close enough to the lower steady state. However, it seems to be meaningful that there is experimental evidence that after cutting a hydra in half, both parts simultaneously regenerate [18]. This means that to model this experiment we should obtain a gradient-like pattern formation for the initial conditions corresponding to both parts of the initial domain (as is the case for the models with diffusion-driven instability described in [12]). It is an open problem whether such behavior can be obtained using a model of the form (8) for any set of parameters and production kinetics. A similar observation holds for the receptor-based model proposed by Sherratt et al. [27]. Results of the simulations described by these authors only concern cutting and regeneration of the head or foot. Simulations of the cutting experiments described in [27] concern cutting of the upper or lower 1/4 of the body column. This means that new simulations are performed for the initial conditions corresponding to 3/4 of the gradient-like solution. For such initial conditions model (8) also reproduces the results of cutting experiments.
112
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
Fig. 10. The solutions of the model for initial condition corresponding to different transplantations of a head.
5.2.2. Grafting experiments Now we formulate the initial conditions describing grafting experiments. Assuming again that v(x) is a solution corresponding to the case of normal development, we put vðxÞ for x 2 ½0; x1 [ ½x2 ; 1; ð26Þ u0 ðxÞ ¼ vðx x1 þ xgraft Þ for x 2 ½x1 ; x2 . This means that, instead of the portion of tissue with coordinates [x1, x2], we transplant the portion of tissue of the same length from the position starting in xgraft. Corresponding simulations can be interpreted as predicting the outcome of grafting experiments. Fig. 10(a) and (b) show the spatio-temporal evolution of the density of bound receptors from initial conditions corresponding to the transplantation of tissue from the upper part of the body column to the lower one. We observe the formation of a second peak corresponding to a second head. The formation and persistence of the second maximum is the result of the bistability of the reaction terms. For such a model we observe various stable solutions, which interpolate between the two stable steady states.
6. Stationary homogeneous steady states In [12] we showed that a receptor-based model (1)–(4) with diffusion-driven instabilities can fully reproduce the outcome of the cutting experiments. Thus, one of the questions arising is whether it is also possible to obtain, using model (8), diffusion-driven instabilities, and how the corresponding parameter values are related to the parameter space for which we obtained the wave solutions. rb . We consider the case with the production of free receptors given by the function pr ¼ m1 1þr b We also discuss results for the constant production pr = m1.
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
113
Solving the algebraic equations for the spatially homogeneous steady state (for details see the Appendix), we obtain the two equations (33) and (34) given in the Appendix. It is interesting to find the steady states for the values of parameters used in simulations for which we obtained spatial patterns i.e. parameter values from the parameter sets I and II (see Table 1). In the case pe = 0, only one real positive value of rb satisfies Eq. (33). Now we assume that pe 5 0. We solve (33) and (34) for the parameter set II. Eq. (33) has a solution which is illustrated in Fig. 11 above using the light-color line. The solution of Eq. (34) is shown in Fig. 11 above using the dark-color line. In summary, for this set of parameters, there is only one positive real solution u1 ¼ ð0:486; 8:456; 17:383; 0:876; 13:208; 0:262Þ. There is also one real solution with u2 ¼ ð2:538; 6:041; 2:379; 0; 0:627; 0ÞÞ. For the complete system we have also the trivial pe = 0 ð solution ð u3 ¼ ð0; 0; 0; 0; 0; 0ÞÞ.
20 Eq. 28 Eq. 29
rb
15
10
5
0 0
2
4
6
8
6
8
pe 20 Eq. 28 Eq. 29
rb
15
10
5
0 0
2
4
pe
Fig. 11. Solutions of (34) crossing the solution of (33) for the parameter set II and a Michaelis–Menten production of m1 r b free receptors pr ¼ 1þr – on the picture above, and constant pr = m1 – on the picture below. We observe that there is b only one common point.
114
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119 Eq. 28 16
Eq. 29
rb
12
8
4
0 0
2
4
6
8
pe
Fig. 12. Solutions of (34) crossing the solution of (33) for the parameter set I and a production of free receptors pr = m1. We observe that there are three common points.
We perform a linear stability analysis of system (27) and find that the trivial solution, u3 , is stable, and that the solution u1 is also stable. The solution u2 is unstable. We investigate whether the steady 2state u1 can be destabilized by diffusion. Thus, we study the eigenvalues of the matrix A D lcm , where A is a linearization of the kinetic system at the state u1 , D is a matrix of diffusion coefficients and l2n;k are wavenumbers obtained from the eigenproblem for the Laplacian. We find that there exist values of the diffusion coefficients and of c for which we can select an unstable mode lm. The simulations performed for such values show that the steady state is indeed destabilized. However, the solution tends to the trivial equilibrium and no spatially inhomogeneous spatial pattern is obtained. We stress here that for the values of diffusion coefficients (and c) for which we obtained a spatially inhomogeneous pattern, the steady state u1 is always stable. We obtained qualitatively similar results for the constant production function pr = m1. For the parameter set I, there exist three positive steady states of (27) (see Fig. 12). The picture looks qualitatively the same as for the constant production of bound receptors. In both cases, two of these steady states are unstable and one is stable and remains stable also for the system with diffusion. We did not find numerically any parameter set such that there is a diffusion-driven instability in system (33), (34) leading to the formation of spatially inhomogeneous patterns.
7. Discussion – A Turing-type model versus a model with hysteresis In the framework of reaction–diffusion systems, there are essentially two ways in which a system of initially identical cells can start to differentiate: • There is a critical number of cells (related to the size of the domain), above which the spatially homogeneous attractor loses stability which leads to ‘spontaneous’ spatial patterning. This is the case for the models with a Turing-type instability.
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
115
• There is an external inducing signal or forcing term which drives the system into a new, spatially inhomogeneous state. Such a signal originates from a transplant formed by a group of already differentiated cells. In such a case the signal must be strong enough to trigger differentiation. This corresponds to a sufficiently large initial perturbation of the homogeneous steady state. This type of initialization of the pattern-forming mechanism is involved in the model with hysteresis. In [12] and in this work we approach the question of mechanisms responsible for the formation of new structures and regeneration of missing ones using the specific example of a fresh-water polyp, the hydra. We consider both mechanisms of pattern formation i.e. diffusion-driven instability and a hysteresis-driven mechanism and demonstrate their advantages and constraints in explaining the different aspects of pattern formation and regulation. The model proposed in this work and the previous Turing-type model of pattern formation in hydra [12] are based on the idea that patterns are controlled by specific cell-surface receptors, which transmit to the cells signals responsible for their differentiation. Our studies show that by taking into account the receptor reactions on the cell surface, activator–inhibitor models of pattern formation can be considerably improved in the sense that they cease to require the artificial hypothesis of the existence of two diffusible morphogens and the autocatalytic production of the molecules. In the previous Turing-type model of pattern formation in hydra [12], it was shown that an arbitrary, ‘random’ initial condition leads to a gradient-like solution. In our model organism, the hydra, such a model can explain the results of cutting experiments. For a certain range of domain size values, we observe the formation of a gradient-like pattern from any initial conditions corresponding to the distributions of molecules on the different parts of a body column of a complete organism. This mimics well the regeneration of the complete animal from an arbitrary piece of the body column which is sufficiently large. However, models with diffusion-driven instabilities (and also the activator–inhibitor models) fail to reproduce the grafting experiments. The present work shows that our model with multiple steady states and hysteretic behavior can explain the results of grafting experiments, because in such a model multiple spatially inhomogeneous solutions may coexist. The solution depends on the initial condition in a manner that is similar to how the results of the grafting experiments depend on the graft position and not on the size of the animal. We have shown that stationary fronts already arise in a simplified model which does not include any equations describing the dynamics of the enzyme or its production. For the complete model we have found numerically similar stationary fronts. However, the dynamics of the complete system is more complicated and we have also observed spatio-temporal patterns oscillating in time. Numerical simulations show that the model with hysteresis cannot reproduce the results of all types of cutting experiments. Only regeneration of a limited interval of a cut head or foot can be simulated (as in [27]). Experiments in which the animal is cut in half and both parts regenerate (compare Fig. 1) are not reproduced by the model with hysteresis (only one of the two halves of the animal is regenerated). However, such experiments are reproduced by the receptor-based system which exhibits a Turing-type instability [12]. The question arises as to whether these two kinds of experiments can be explained using the same mechanism or whether it could be a combination of previously considered mechanisms.
116
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
We checked that in model (8) diffusion-driven instabilities do not lead to the formation of a gradient-like pattern. However, the behavior of the system for different parameter sets and different kinetics describing the production of ligands and enzyme is still an open problem. As in many other cases the complexity of the biological process makes it difficult to give a rigorous derivation of the model equations. In this work we have proposed heuristic kinetics for the production of new ligands and enzyme molecules. These kinetics aim to summarize the interactions of a complicated cascade of biochemical intracellular and extracellular reactions. The main aim is to check whether a hysteretic mechanism arising from the existence of multiple steady states and switches in the dynamics is capable of explaining the observed experimental results which so far have not been able to be explained by any existing models. It would be interesting to confirm experimentally the existence of a hysteretic dependence of ligand production and their concentration in the regulatory system of hydra, in a way analogous to the work of Tyson and Novak [32] in the context of the cell cycle. This would manifest itself in rapid changes of the production rate under a continuous change of the concentration of the molecules. Further biological results are needed to gain insight into the molecular nature of cell communication and of positional information. To understand morphogenesis in hydra it is necessary to bridge the gap between experimental observations at the cellular level and at the genetic and biochemical levels. So far, about 20 ligand-activated signaling pathways (among others Wnt, TGF-b, Hedgehog and Notch) have been identified in hydra, in which at least one component is likely to play a developmental role [29]. In some cases, there is definitive experimental evidence for a developmental role as a head activator for example. In other cases, such a role is assumed based on the results obtained for other organisms. Recently, the Wnt signaling cascade has become the focus of much research attention and fragments of candidate receptor-genes were found in the cDNA library [5]. The question of how the Wnt regulatory network may be responsible for spatial patterning is the subject of our further work. We aim to build more detailed models based on the experimental results on the Wnt signaling cascades in hydra cells performed in the lab of Prof. Holstein (Zoology Institute, University of Heidelberg). Studies with heuristic kinetics suggest possible a pattern formation scenario and show that the existence of even one diffusing substance in the complicated signaling network can lead to the formation of different spatial structures. On the other hand, the constraints of the models suggest that communication over a distance may involve cell-to-cell signaling mechanisms different from simple diffusion. This is why there arises the need to explain how the intracellular and extracellular interactions are translated into macroscopic form and pattern. The complicated structures of the cell cytoplasm, of the cell nucleus and of the extracellular matrix make it evident that modeling of transport processes must be handled carefully, in many cases reflecting the porous structure of the medium [8]. This current work shows that it is justified to look for complicated nonlinear dynamics in the signaling pathways leading to synthesis of new ligands. The models proposed in this paper have shown how the nonlinearities of the signaling pathways give rise to the final spatial patterning. They may be interesting for biologists who are looking for a possible scenario of pattern formation in hydra, and also for mathematicians looking for the mechanisms involved in reaction–diffusion systems which would lead to the formation of spatially inhomogeneous solutions consistent with experiments.
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
117
Acknowledgment I would like to thank Prof. Willi Ja¨ger and Prof. Werner Mu¨ller for supporting me and Prof. Marek Kimmel for valuable comments on the manuscript. This work was supported by the grant Graduiertenkolleg 484: ‘Signalling systems and gene expression in model systems of developmental biology’. Appendix Stationary homogeneous steady states of the model (8) with a Michaelis–Menten production of rb satisfy the system, the free receptors, pr ¼ m1 1þr b rb 0 ¼ lf rf brf l þ drb þ m1 ; 1 þ rb 0 ¼ lb rb þ brf l drb ; 0¼
1 o2 l ll l brf l þ drb þ pl be le; c ox2
d 2 o2 e le e þ p e ; c ox2 pl m2 lrb þ 0 ¼ dl ; 1 þ p2l ð1 þ rl p2l bl pl Þð1 þ al rb Þ pe m3 le þ . 0 ¼ de 2 2 1 þ pe ð1 þ re pe be pe Þð1 þ ae rb Þ
ð27Þ
0¼
We solve the 1st, 2nd and 4th equation with respect to rf, l and e. Assuming that rb 5 0, we obtain p e¼ e; ð28Þ le rb ðlb þ lb rb m1 Þ ; ð29Þ rf ¼ lf ð1 þ rb Þ l¼
ðlb þ dÞlf ð1 þ rb Þ . bðlb þ lb rb m1 Þ
ð30Þ l d ðb2 þr l2 l2 b l l b Þ
For rb = 0 we can calculate that the steady states are (0, 0, 0, 0, 0, 0) and ð0; 0; e e em3ðbe 2eþll 2 l2eÞ e l e , e l e lbel ; 0; lbe le l Þ, the second of which is non-positive. Substituting (28)–(30) in (27) we reduce system (27) to the following system of three equations: ll ðlb þ dÞlf ð1 þ rb Þ be ðlb þ dÞlf ð1 þ rb Þpe ; þ drb þ bðlb þ lb rb m1 Þ bðlb þ lb rb m1 Þle m2 rb ðlb þ dÞlf ð1 þ rb Þ dl p l 0¼ ; 2 ð1 þ pl Þ bðlb þ lb rb m1 Þð1 þ al rb Þð1 þ rl p2l bl pl Þ m3 ðlb þ dÞlf ð1 þ rb Þpe de p e . 0¼ 2 ð1 þ pe Þ bðlb þ lb rb m1 Þle ð1 þ ae rb Þð1 þ re p2e be pe Þ
0 ¼ pl rb ðlb þ dÞ þ
ð31Þ
118
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
We solve the first of the above equations with respect to pl and obtain pl ¼ ðbe lf pe lb be lf pe d be lf pe lb rb be lf pe drb þ rb ble l2b þ r2b ble l2b rb ble lb m1 ll lf le lb ll lf le d ll lf le lb rb ll lf le drb Þ=bðlb þ lb rb m1 Þle .
ð32Þ
Finally we obtain two equations 0 ¼ dl ðbe lf pe lb belf pe d be lf pe lb rb belf pe drb þ rb ble l2b þ r2b ble l2b rb ble lb m1 ll lf le lb ll lf le d ll lf le lb rb ll lf le drb Þ=b=ðlb þ lb rb m1 Þ=le =ð1 þ ðbe lf pe lb belf pe d belf pe lb rb belf pe drb þ rb ble l2b 2
þ r2b ble l2b rb ble lb m1 ll lf le lb ll lf le d ll lf le lb rb ll lf le drb Þ =b2 =ðlb þ lb rb m1 Þ2 =l2e Þ m2 rb ðlb þ dÞlf ð1 þ rb Þ=b=ðlb þ lb rb m1 Þ =ð1 þ al rb Þ=ð1 þ rl ðbe lf pe lb be lf pe d be lf pe lb rb be lf pe drb þ rb ble l2b þ r2b ble l2b rb ble lb m1 ll lf le lb ll lf le d ll lf le lb rb 2
2
ll lf le drb Þ =b2 =ðlb þ lb rb m1 Þ =l2e bl ðbelf pe lb be lf pe d be lf pe lb rb be lf pe drb þ rb ble l2b þ r2b ble l2b rb ble lb m1 ll lf le lb ll lf le d ll lf le lb rb ll lf le drb Þ=b=ðlb þ lb rb m1 Þ=le Þ; 0¼
m3 ðlb þ dÞlf ð1 þ rb Þpe de pe þ . 2 ð1 þ pe Þ bðlb lb rb þ m1 Þle ð1 þ ae rb Þð1 þ re p2e be pe Þ
ð33Þ ð34Þ
References [1] A. Babloyantz, M. Sanglier, Chemical instabilities of ‘‘All-or-None’’ type in b-galactosidase induction and active transport, FEBS Lett. 23 (1972) 346. [2] F.R. Cross, V. Archambault, M. Miller, M. Klovstad, Testing a mathematical model of the yeast cell cycle, Mol. Biol. Cell 13 (2002) 52. [3] A. Gierer, H. Meinhardt, A theory of biological pattern formation, Kybernetik 12 (1972) 30. [4] S. Heinze, B. Schweizer, H. Schwetlick, Existence of front solutions in degenerate reaction diffusion systems, SFB preprint, University of Heidelberg, 2004. [5] B. Hobmayer, F. Rentzsch, K. Kuhn, C.M. Happel, C.C. Laue, P. Snyder, U. Rothba¨cher, T.W. Holstein, Wnt signaling and axis formation in the diploblastic metazoan Hydra, Nature 407 (2000) 186. [6] F. Hoppensteadt, W. Ja¨ger, Pattern formation by bacteria, in: S. Levin (Ed.), Lecture Notes in Biomathematics: Biological Growth and Spread, Springer, Heidelberg, 1980. [7] F. Hoppensteadt, W. Ja¨ger, C. Po¨ppe, A hysteresis model for bacterial growth patterns, in: S. Levin (Ed.), Lecture Notes in Biomathematics: Modelling of Patterns in Space and Time, Springer, Heidelberg, 1983. [8] U. Hornung, W. Ja¨ger, A. Mikelic, Reactive transport through an array of cells with semi-permeable membranes, Mathematical Modeling and Numerical Analysis 28 (1994) 59. [9] J. Keener, J. Sneyd, Mathematical Physiology, Springer, New York, 1998. [10] Ch.T. Klein, Hysteresis-driven structure formation in biochemical networks, J. Theor. Biol. 194 (1998) 263. [11] D.A. Lauffenburger, J.J. Linderman, Receptors. Models for Binding, Traficking, and Signaling, Oxford University, New York, 1993.
A. Marciniak-Czochra / Mathematical Biosciences 199 (2006) 97–119
119
[12] A. Marciniak-Czochra, Receptor-based models with diffusion-driven instability for pattern formation in hydra, J. Biol. Syst. 11 (2003) 293. [13] J.-L. Martiel, A. Goldbeter, A model based on receptor desensitization for cyclic AMP signaling in Disctyostelium cells, Biophys. J. 57 (1987) 807. [14] H. Meinhardt, A model for pattern formation of hypostome, tentacles and foot in hydra: how to form structures close to each other, how to form them at a distance, Dev. Biol. 157 (1993) 321. [15] P.B. Monk, H.G. Othmer, Cyclic AMP oscillations in suspensions of Disctyostelium Discoideum, Philos. Trans. R. Soc. Lond. B323 (1989) 185. [16] J.D. Murray, Mathematical Biology, Springer, Berlin, 1993. [17] W.A. Mu¨ller, Ectopic head and foot formation in hydra. Diacylglycerol induced increase in potential value and assistance of the head in foot formation, Differentiation 42 (1990) 131. [18] W.A. Mu¨ller, Pattern control in hydra: basic experiments and concepts, in: H.G. Othmer, P.K. Maini, J.D. Murray (Eds.), Experimental and Theoretical Advances in Biological Pattern Formation, Plenum, New York, 1993. [19] W.A. Mu¨ller, Competition of factors and cellular resources as a principle of pattern formation in hydra, Dev. Biol. 167 (1995) 175. [20] W.A. Mu¨ller, Developmental Biology, Springer, New York, 1997. [21] S.C. Mu¨ller, G. Venzl, Pattern formation in precipitation processes, in: S. Levin (Ed.), Lecture Notes in Biomathematics: Modelling of Patterns in Space and Time, Springer, Heidelberg, 1983. [22] B. Novak, J. Tyson, Numerical analysis of a comprehensive model of M-phase control in Xenopus oocyte extracts and intact embryos, J. Cell. Sci. 106 (1993) 1153. [23] B. Novak, Z. Pataki, Mathematical model of cell division cycle of fission yeast, Chaos 11 (2001) 277. [24] F. Rothe, Global Solutions of Reaction–Diffusion Systems, Springer, Berlin, 1984. [25] F.F. Seelig, B. Denzle, Hysteresis without autocatalysis in enzyme systems with substrate inhibition, FEBS Lett. 24 (1972) 283. [26] W. Sha, J. Moore, K. Chen, A. Lassaletta, Y. Chung-Seon, J. Tyson, Hysteresis drives cell-cycle transitions in Xenopus Laevis egg extracts, PNAS 100 (2003) 975. [27] J.A. Sherratt, P.K. Maini, W. Ja¨ger, W. Mu¨ller, A receptor-based model for pattern formation in hydra, Forma 10 (1995) 77. [28] J. Smoller, Shock-Waves and Reaction–Diffusion Equations, Springer, New York, 1994. [29] R. Steele, Developmental signaling in hydra: what does it take to build a ‘‘simple’’ animal? Developmental Biology 248 (2002) 199. [30] U. Technau, B. Hobmayer, F. Rentzsch, T. Holstein, Molecular and cellular analysis of de novo pattern formation in Hydra, in: A. Deutsch (Ed.), Function and Regulation of Cellular Systems: Experiments and Models, Birkha¨user-Verlag, Basel, 2002. [31] C.D. Thron, Bistable biochemical switching and the control of the events of the cell cycle, Oncogene 15 (1997) 317. [32] J.J. Tyson, B. Novak, Regulation of the eukaryotic cell-cycle: molecular antagonism, hysteresis, and irreversible transitions, J. Theor. Biol. 210 (2001) 249. [33] L. Wolpert, Positional information and the spatial pattern of cellular differentiation, J. Theor. Biol. 25 (1969) 1.