Reaction-Diffusion Models of Biological Pattern Formation: Effects of Domain Growth and Time Delays

1. Reaction Diffusion Models of Biological Pattern Formation:The Effects of Domain Growth and Time Delays EA Gaffney Collaborators: NAM Monk, E Crampin, PK Maini

2. Summary of Background • In 1952 Turing proposed that • Pattern Formation during morphogenesis might arise through • an instability in systems of reacting chemicals (morphogens), • which is driven by diffusion. • Heterogeneous concentrations of these chemicals form a “pre-pattern”. • Subsequent differentiation of tissue/cell type is in response to whether • or not one of these morphogens exceeds some threshold locally. • The equations describing this for two reacting constituents on a • stationary domain are of the form

3. In the simplest setting there is a unique homogeneous steady state at (a0, b0), given by solutions of The Jacobean at the stationary point is of the form The kinetics are always chosen such that, in the absence of diffusion, the homogeneous steady state is stable (and thus the instability is diffusion driven) In the presence of diffusion, for a sufficiently large domain, and quite reasonably assuming the components of the Jacobean are O(1) compared to the scales ε, 1/ε the homogeneous steady state is unstable if

4. For • a domain larger than the critical size, again assuming the components of the Jacobean are O(1) compared to the scales ε , 1/ε • ε <<1sufficiently small • the rate of growth of the fastest growing mode, μ, is given by • where

5. A key feature of all Turing-Pair kinetics is • Morphogen induced production of morphogen • Short range activation, long range inhibition. A specific example: Schnakenberg Kinetics p = 0.9, q = 0.1, ε << 1. 1/[T0 ] = λ is the decayrate of the activator, b In particular the activator production is equivalent to a law of mass action rule

6. Potential Examples • Avian feather bud formation • HS Jung, …, L Wolpert, ... et al, Developmental Biology 1998 • Vertebrate limb formation • CM Leonard et al, Developmental Biology, 1991. • TGF-β2 possible activator; inhibitor undetermined • Zebrafish mesendoderm induction • L Solnica-Kreznel, Current Biology, 2003 • Nodal (Squint) gene product is an activator. • Lefty gene product is an inhibitor • Evidence that range of Lefty's influence exceeds Nodal's. • Molecular details remain to be uncovered of their interaction (though progress is rapid on this point) • Molecular details remain to be determined on the differential in • their range of influence.

7. Difficulties • Reaction Diffusion Patterns can be observed to be very sensitive to • Noise (in initial conditions and generally). • Perturbations of the domain shape. • Turing Morphogens are hard to find. • However, more and more molecular data is being produced in developmental studies. These indicate that a possible Turing Pair are Nodal and Lefty in Zebrafish mesendodermal induction. • The molecular data also indicates that Nodal and Lefty and other putative Turing pairs induce each others' production by signal transduction and gene expression. • For example, in situ hybridisation reveals mRNA transcripts of the proteins speculated to be Turing pairs. • The extracellular domain is complicated and tortuous. • The precise details of the kinetic functions are only ever speculated.

8. Formulating a model … Complicated extracellular domain

9. Robustness and Reaction diffusion on growing domains The pattern produced by an RD system can be sensitive to the details of the initial conditions 0.3 frequency 0.2 0.1 Mode Number 5 10 15 20 25 • Numerical Simulations of Kondo, Asai, Goodwin and others • indicate that • domain growth can lead to robust pattern formation, ie. an • insensitivity to noise and randomness in the initial conditions. • This has previously motivated a detailed investigation of • the stabilising influence of domain growth • the mechanisms by which it produces robustness • conditions for which one may expect robust pattern formation

10. Model Formulation: Incorporating uniform domain growth Uniform growth Rescaling gives

11. Exponential Domain Growth Self similarity arguments indicate this behaviour will continue indefinitely in time. The pattern is insensitive to details of the initial conditions These observations hold over five - six orders of magnitude of domain growth. The robustness is insensitive to the details of the kinetics (providing pattern initially forms). Linear Domain growth Frequency doubling behaviour breaks down more readily. No self - similarity arguments. Activator Exponential Growth Schnakenberg Kinetics 50 Space Space 30 10 500 4000 Time Activator Linear Growth Schnakenberg Kinetics 30 Space 10 Time 2

12. Pattern formation for logistic growth (Schnakenberg Kinetics). • For the exponential phase of logistic growth, robustness is observed. • However, there is the possibility of a loss of robustness as the domain growth saturates, as above for the centre plot. • However, the saturation domain size increases by a factor of 1.0015 on moving from left to right. • This indicates an extreme fine tuning of parameters is required to loose robustness for logistic growth.

13. Conclusions: Robustness and Domain Growth • Including GROWTH in 1D reaction diffusion models leads to • Robustness to noise in the initial conditions over 4-5 orders • of magnitude of the growth rate for exponential growth • Semi-scale invariance. No need for parameter fine tuning or • feedback between the domain size and the kinetics. • Persistence of robustness for logistic saturating growth. • Robustness independent of the exact details of the model • providing pattern initially forms. Dismissal of Reaction Diffusion as a Pattern Formation mechanism on the grounds of robustness not necessarily founded if slow domain growth is relevant.

14. Time Delays

15. Question: How are the effects of signal transduction and gene expression time delays incorporated into a model? Signal

16. In more detail ... Transcription and Translation delay dependent on size on protein. It is at least 10 minutes and can be several hours. Signal transduction serves to only increase this delay

17. For a suitable non-dimensionalisation, the Schnakenberg reaction diffusion equations, in the presence of domain growth and time delays, can be written in the form A, 2B lost from reaction at time t 3B gained from reaction at time t-τ u is the velocity field of the domain growth y takes values in [0, L(t)], where L(t) is the (non-dimensionalised) domain length τ is the gene expression time delay where

19. Model Investigations Linear Analysis. (a*, b*) is the homogeneous steady state solution Substitute into the model equations to obtain (on neglecting O(η2)) –

20. In the absence of a time delay, with or without domain growth, substitute • The critical value of the domain length is given via • There are no oscillations at this critical point, nor for the fastest growing mode (at least providing ε2 is sufficiently small). • The rate of growth of the fastest growing mode, μ, is again given by into Schnakenberg Kinetics =

21. In the absence of domain growth, substitute into + to obtain • For ε <<1sufficiently small and q11 = q12 = 0, q22 > 0, p22+q22 > 0, as with Schnakenberg kinetics, one has • There are no oscillations at the critical domain length • The critical length for the onset of instability is unchanged by the time delay

22. In the absence of domain growth, substitute into + to obtain • For ε <<1sufficiently small and q11 = q12 = 0, q22 > 0, p22+q22 > 0, as with Schnakenberg kinetics, one has • Writing λ = μ + iν, the rate of growth of the fastest growing mode, μ, is given implicitly by

23. We have the fastest growing mode has a growth rate given by while in the absence of domain growth the fastest growth rate is given by the above expression with τ = 0. 1 We focus on non-oscillatory modes (ν = 0), as oscillations are not observed in Schnakenberg kinetics One can solve the above in terms of the LambertW functions on neglected the O(ε2) terms. On the left is a plot of the ratio μ(τ, ν = 0) / μ(τ = 0, ν = 0) as a function of τq22, τ(p22-ε2π2/γ) 0.8 0.6 0.4 0.2 0 -1 τ(p22-ε2π2/γ)

24. We have the fastest growing mode has a growth rate given by while in the absence of domain growth the fastest growth rate is given by the above expression with τ = 0. 1 We focus on non-oscillatory modes (ν = 0), as oscillations are not observed in Schnakenberg kinetics One can solve the above in terms of the LambertW functions on neglected the O(ε2) terms. On the left is a plot of the ratio μ(τ, ν = 0) / μ(τ = 0, ν = 0) as a function of τq22, τ(p22-ε2π2/γ) 0.8 0.6 0.4 0.2 0 -1 τ(p22-ε2π2/γ)

25. Explicitly Solving the Linear Equations to obtain insight for when there are both Time delays and Domain Growth

26. We can see that time delays do greatly increase the time it takes to leave the homogeneous steady state, as indicated by analysis A sufficiently large value of τδ = 2ln(2) [Time Delay/Doubling Time] can result in the large time asymptote of the linear theory decaying to zero. These results appear to be completely general.

27. τ/τ0 τ/τ0 τδ/τ0 τδ/τ0 The minimum of, say, 50 and the large time asymptotic value of the components of An=1 are plotted against δ, τ/τ0,τδ/τ0 for various parameters. Note that the large time asymptote is always small for sufficiently large τδ.

28. Thus one cannot rely on a naive linear analysis predicting an instability via growth away from the homogeneous steady state. The large time asymptote may decay to zero for sufficiently large τδ. Whether the intermediate behaviour triggers pattern formation depends on the non-linear dynamics

29. Conclusions from the linearised equations • Domain Growth, No Time Delay. In the absence of time delays, domain growth does not have much effect on the linear analysis • Time delays, No Domain Growth. There will typically be a substantial patterning lag • The location of the onset of the instability is independent of the time delay • The ratio of the fastest growing modes in the presence of the time and in the absence of the time delay will typically be large. • Domain Growth & Time Delays. • There will typically be a substantial patterning lag • A naive linear analysis is conceptually flawed for the prediction of instability. Whether the intermediate behaviour triggers pattern formation depends on the non-linear dynamics • The behaviour of the large asymptote is governed by the parameter • τδ = 2ln(2) [Time Delay/Doubling Time]

30. Numerical Simulations of the Nonlinear equations with Schankenberg kinetics and domain growth and time delays

31. Initial Conditions for t in the domain [0,τ]: a b x x The initial conditions are typically given by the solid lines above (IC=1) and the dashed lines (IC=2). We also consider multiplying these initial conditions by time dependent factors. One example is [1+ 0.0025 cos (πx) cos (πt/(2τ))] All the behaviour observed below is representative of the numerous initial conditions considered. Similarly for an order of magnitude variation of the parameters τ, ε, δ.

32. τ = 0, IC = 1 τ = 0, IC = 2 Stationary Domain Gray Scale plots of the activator (Schnakenberg) There are no oscillations. The final pattern is sensitive to the details of the initial conditons. τ0 corresponds to a delay of 12 minutes in the dimensional model A delay of τ0 induces a patterning lag of about 60τ0 A delay of 4τ0 induces a patterning lag of about 240τ0 τ = τ0, IC = 1 τ = τ0, IC = 2 τ = 4τ0, IC = 1 τ = 4τ0, IC = 2 x x

33. τ = 0, IC = 1 τ = 0, IC = 2 103 103 Growing Domain Gray Scale plots of the activator (Schnakenberg) Domain doubling time of 2 days There are no oscillations. τ0 corresponds to a delay of 12 minutes in the dimensional model Time delays delay the onset of peak doubling γ/γ(t = 0) τ = τ0, IC = 1 τ = τ0, IC = 2 103 103 γ/γ(t = 0) τ = 4τ0, IC = 1 τ = 4τ0, IC = 2 103 103 γ/γ(t = 0) x x

34. τ = 0, IC = 1 τ = 0, IC = 2 103 103 Growing Domain Gray Scale plots of the activator (Schnakenberg) Domain doubling time of 2 days There are no oscillations. τ0 corresponds to a delay of 12 minutes in the dimensional model Time delays delay the onset of peak doubling Larger time delays result in the absence of the Turing instability γ/γ(t = 0) τ = τ0, IC = 1 τ = τ0, IC = 2 103 103 γ/γ(t = 0) τ = 8τ0, IC = 1 τ = 8τ0, IC = 2 103 103 γ/γ(t = 0) x x

35. τ = 0, IC = 1 τ = τ0, IC = 1 γ/γ(t = 0) τ = 2τ0, IC = 1 τ = 4τ0, IC = 1 γ/γ(t = 0) γ/γ(t = 0) x x • Growing Domain • The time delay induces a delay to patterning • A time delay of τ0, i.e. 12 minutesinduces a lag of a domain doubling time, i.e. 2 days

36. τ = 4τ0, IC = 1 τ = 4τ0, IC = 1 103 103 Domain Doubling Time: 2 days γ/γ(t = 0) x x τ = τ0, IC = 2 τ = 4τ0, IC = 1 103 τ = τ0, IC = 1 103 103 5 Domain Doubling Time: 12 hours γ/γ(t = 0) x • Growing Domain • The behaviour of the system appears to be governed by the parameter grouping • τδ = 2ln(2)[Time Delay/Doubling Time]

37. τ = 4τ0, IC = 1 τ = 4τ0, IC = 1 103 103 Domain Doubling Time: 2 days γ/γ(t = 0) τ = τ0, IC = 2 τ = 4τ0, IC = 1 103 τ = τ0, IC = 1 103 103 5 Domain Doubling Time: 12 hours γ/γ(t = 0) x x x • Growing Domain • The behaviour of the system appears to be governed by the parameter grouping • τδ = 2ln(2)[Time Delay/Doubling Time]

38. τ = 4τ0, IC = 1 τ = 4τ0, IC = 1 103 103 Domain Doubling Time: 2 days γ/γ(t = 0) τ = 16τ0, IC = 1 τ = 16τ0, IC = 2 τ = 4τ0, IC = 1 103 103 γ/γ(t = 0) Domain Doubling Time: 8 days γ/γ(t = 0) x x x x • Growing Domain • The behaviour of the system appears to be governed by the parameter grouping • τδ = 2ln(2)[Time Delay/Doubling Time]

39. x x x x

40. Irregular behaviour also possible, along with oscillations, before the failure of the Turing instability as one increases the time delay. x x x x

41. Schakenberg Model. Results. Summary • Stationary Domain • No oscillations generally • Pattern is sensitive to the initial conditions • Time delays can induce a large patterning lag • Growing Domain Results. Summary: • No oscillations generally • Time delays can induce a large patterning lag • Time delays can induce irregular behaviour and a failure of the • Turing instability • The behaviour of the system is governed by the parameter • grouping: τδ = 2ln(2)[Time Delay/Doubling Time]

42. Conclusions. Time Delays and Biological RD Systems • We have • motivated the biophysical need for the inclusion of • signal transduction and gene expression time delays in models • of biological pattern formation • We have demonstrated how these delays can be included in • one of the simplest "long range activation-short range inhbition" • pattern forming reaction diffusion models. • While we have not typically found oscillations, we have found that • Time delays can make a large difference to the patterns emerging • from the models, especially with regard to patterning lags and, • for growing domains, the failure of the Turing instability • Thus …

43. The above observations do not rule out reaction diffusion as a • putative pattern formation mechanism, whether on a stationary • or uniformly growing spatial domain. • However, when considering patterning events especially those for • which rapid establishment of pattern is critical, such as in the tissues • of developing embryos, our results show that • any putative time delays cannot be neglected in general without • careful justification. • our finding that time delays can dramatically increase the • time taken for the reaction diffusion system to initiate patterns • imposes potentially severe constraints on the potential • molecular details of any Turing system that might operate • during developmental patterning.

44. Future Work/Further Questions • Continue investigating the extent to which the results are general, especially for models with "short range activation, long range inhibition" • Kinetics with a negative feedback loop,e.g. Gierer Meinhardt • Kinetics with more than two componenents and multiple time delays. Are the patterning lags cumalative? • Other biological pattern formation mechanisms e.g. the mechanochemical models • If the results are general • We have, in general, potentially severe constraints on the reaction diffusion mechanism and other mechanisms of biological pattern formation • If the results are not general • there is a clear distinction between the patterning forming behaviour of "short range activation, long range inhibition" on the inclusion of time delays. This would have a substantial impact in that the choice of the kinetics really does matter!

45. Complicated extracellular domain. An exercise in homogeneisation theory:

46. The End