Non-Autonomous Reaction Diffusion Equations Biological Pattern

Document Sample
Non-Autonomous Reaction Diffusion Equations Biological Pattern Powered By Docstoc
					Reaction Diffusion Models of Biological Pattern Formation:
     The Effects of Domain Growth and Time Delays


                       EA Gaffney

     Collaborators: NAM Monk, E Crampin, PK Maini
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
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
For
      • a domain larger than the critical size, again assuming the components
        of the Jacobean are O(1) compared to the scales ε , 1/ε
      • ε <<1 sufficiently small
the rate of growth of the fastest growing mode, μ, is given by




where
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 decay rate of the activator, b

In particular the activator
production is equivalent to
a law of mass action rule
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.
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.
Formulating a model …




 Complicated extracellular domain
Robustness and Reaction diffusion on growing domains
The pattern produced
                                           0.3




                           frequency
by an RD system can
be sensitive to the                        0.2
details of the initial
conditions                                 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
Model Formulation: Incorporating uniform domain growth




Uniform growth



Rescaling



gives
Exponential Domain Growth                                 Activator
Self similarity arguments indicate     50            Exponential Growth
this behaviour will continue                        Schnakenberg Kinetics
indefinitely in time.




                                            Space
                                            Space
                                       30
The pattern is insensitive to
details of the initial conditions
These observations hold over           10
five - six orders of magnitude of
domain growth.                               500             Time            4000

The robustness is insensitive to
the details of the kinetics                                Activator
(providing pattern initially forms).   30               Linear Growth
                                                     Schnakenberg Kinetics
                                            Space
Linear Domain growth
Frequency doubling behaviour           10
breaks down more readily. No
self - similarity arguments.
                                                            Time               2
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.
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.
Time Delays
Question: How are the effects of signal transduction and
gene expression time delays incorporated into a model?

                         Signal
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
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
                 Model Investigations
Linear Analysis.                               (a*, b*) is the homogeneous
                                               steady state solution

Substitute


into the model equations to obtain (on neglecting O(η2))




      –
In the absence of a time delay, with or without domain growth, substitute




into


 • The critical value of the domain length is given via             Schnakenberg
                                                                    Kinetics

                                          =

 • 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
In the absence of domain growth, substitute                             into


                                               +

 to obtain




 For ε <<1 sufficiently 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
In the absence of domain growth, substitute                           into


                                              +

 to obtain




 For ε <<1 sufficiently 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
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.

We focus on non-oscillatory                                                     1
modes (ν = 0), as                                                               0.8
oscillations are not observed
                                                                                0.6
in Schnakenberg kinetics
                                                                                0.4
One can solve the above in                                                      0.2
terms of the LambertW
functions on neglected the                                                      0
O(ε2) terms.
On the left is a plot of the
ratio
μ(τ, ν = 0) / μ(τ = 0, ν = 0)
                                                                                 -1
as a function of
      τq22, τ(p22-ε2π2/γ)                           τ(p22-ε2π2/γ)
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.

We focus on non-oscillatory                                                     1
modes (ν = 0), as                                                               0.8
oscillations are not observed
                                                                                0.6
in Schnakenberg kinetics
                                                                                0.4
One can solve the above in                                                      0.2
terms of the LambertW
functions on neglected the                                                      0
O(ε2) terms.
On the left is a plot of the
ratio
μ(τ, ν = 0) / μ(τ = 0, ν = 0)
                                                                                 -1
as a function of
      τq22, τ(p22-ε2π2/γ)                           τ(p22-ε2π2/γ)
Explicitly Solving the Linear Equations
  to obtain insight for when there are
both Time delays and Domain Growth
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.
                                                                   τ/τ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 τδ.
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
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]
      Numerical Simulations
    of the Nonlinear equations
    with Schankenberg kinetics
and domain growth and time delays
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 τ, ε, δ.
                              τ = 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   τ = τ0, IC = 1    τ = τ0, IC = 2
the initial conditons.
τ0 corresponds to a delay
of 12 minutes in the
dimensional model
A delay of τ0 induces a       τ = 4τ0, IC = 1   τ = 4τ0, IC = 2
patterning lag of about
60τ0
A delay of 4τ0 induces a
patterning lag of about
240τ0
                                    x                  x
                                               103       τ = 0, IC = 1 103    τ = 0, IC = 2
Growing Domain




                             γ/γ(t = 0)
Gray Scale plots of the
activator (Schnakenberg)
Domain doubling time of
2 days                                                τ = τ0, IC = 1           τ = τ0, IC = 2
                                               103                    103
There are no oscillations.


                             γ/γ(t = 0)
τ0 corresponds to a delay
of 12 minutes in the
dimensional model
Time delays delay the                          103 τ = 4τ0, IC = 1    103 τ = 4τ0, IC = 2
onset of peak doubling
                                  γ/γ(t = 0)




                                                      x                        x
                                               103       τ = 0, IC = 1 103    τ = 0, IC = 2
Growing Domain




                             γ/γ(t = 0)
Gray Scale plots of the
activator (Schnakenberg)
Domain doubling time of
2 days                                                τ = τ0, IC = 1           τ = τ0, IC = 2
                                               103                    103
There are no oscillations.


                             γ/γ(t = 0)
τ0 corresponds to a delay
of 12 minutes in the
dimensional model
Time delays delay the                          103 τ = 8τ0, IC = 1    103 τ = 8τ0, IC = 2
onset of peak doubling
                                  γ/γ(t = 0)




Larger time delays result
in the absence of the
Turing instability
                                                      x                        x
                      τ = 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 minutes induces a lag of a domain doubling
       time, i.e. 2 days
                               103      τ = 4τ0, IC = 1 103       τ = 4τ0, IC = 1




                  γ/γ(t = 0)
Domain Doubling
Time: 2 days


                                              x                          x
                               103   τ = τ0, IC = 1     = 4τ τ = 0
                                                       τ10330, IC= τ1, IC = 2
                                                        10
                               5
                  γ/γ(t = 0)


Domain Doubling
Time: 12 hours


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




                  γ/γ(t = 0)
Domain Doubling
Time: 2 days



                               103   τ = τ0, IC = 1     = 4τ τ = 0
                                                       τ10330, IC= τ1, IC = 2
                                                        10
                               5
                  γ/γ(t = 0)


Domain Doubling
Time: 12 hours


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




                  γ/γ(t = 0)
Domain Doubling
Time: 2 days



                               103         τ = 16τ0, IC = 1       = 4τ τ = 1
                                                                 τ103 0, IC= 16τ0, IC = 2




                                      γ/γ(t = 0)
                  γ/γ(t = 0)


Domain Doubling
Time: 8 days


                                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]
x   x   x   x
Irregular behaviour also possible, along with oscillations, before the
   failure of the Turing instability as one increases the time delay.




              x                  x                 x                     x
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]
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 …
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.
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!
Complicated extracellular domain. An exercise in homogeneisation theory:
The End
                           Caveats

2 Reaction Diffusion Patterns are observed to be very sensitive to

   • Perturbations of the domain shape.
   • Noise (in initial conditions and generally).
                       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.


   – Heterogenous 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.
        Typical Patterns Formed by This Mechanism


  Such a mechanism can often
  produce plausible looking results

      – Mammal coat patterns
                           (Murray)



And interesting, non-trivial, predictions

    – You shouldn’t be able to find a
      mammal with a striped body and
      a spotted tail.
                   Domain Growth

Numerical Simulations of Kondo, Asai, Goodwin and others
indicate that

   • domain growth can lead to robust pattern formation, ie.
     insensitivity to noise and randomness in the initial conditions.

& motivates 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
      Arcuri & Murray Conclusions

By observations of numerical simulations

   • domain growth yields a tendency towards robustness, but
     this is far from universal.

This had led some theoretical biologists to dismiss Reaction Diffusion
Even as a possible/plausible pattern formation mechanism (eg Saunders
& Ho, Bull Math Bio, 1995).

BUT
   • One should derive the Reaction Diffusion Equations with domain
     growth from conservation laws, not crude scaling arguments.
If at  = * we find
                c( x,  * )  c(T ( x),  * / 4)            “Lock in”

then, given PDE uniqueness,

          c( x,  * )  c(T ( x),  * / 4),    * .

Thus we have frequency doubling cascades of self-similar patterns.

Of course, one cannot expect


                   c( x,  * )  c(T ( x),  * / 4)
to hold exactly in practice so a form of “start-near, stay near”
stability is also required for a frequency doubling cascade.
   Domain Growth: Initial Conclusions


Incorporation of exponential domain growth initially appears to give

   • Robust pattern Formation, ie. No dependence on details
     of initial conditions

   • No need to fine tune the parameters

   • No need to fine tune the model.
       Breakdown of Mode Doubling


Mode Doubling breaks for high and low growth rates.

   • Exponential Growth, Schnakenburg Kinetics




                              Growth Rate
   d  0.01                                0.05
                              Reaction Rate
          Breakdown of Mode Doubling

Breakdown at High Growth rates appears to be due to the fact

   • Further peak reorganisation occurs before splitting peak
     reaches quasi-steady state.

   • This in turn prevents Lock in, ie a point where

                c( x,  * )  c(T ( x),  * / 4)
   • As     1 we get a total breakdown of pattern.
For Low growth rates we also get breakdown …
   Low Growth Rate Mode Doubling Breakdown


Low Growth rate frequency doubling breakdown is

   • noise dependent
   • occurs after lock in




   • d  0.01                 10-6
          General Conclusion


Robust Pattern Formation occurs for exponential
growth with Schnakenburg kinetics over 4-5 orders
of magnitude of the growth rate.
   • Do not expect self similar cascade to proceed indefinitely
     for linear growth

   • Given Stability assumptions, if there is a   *   such that

              c(T ( x),  * / 4,  / 2) c( x,  * ,  )
   or equivalently

               c(T ( x),  * / 4,  ) c( x,  * , 2  )
   • then these coincide for     *.
Prediction:
   • If a sequence with a linear growth rate  undergoes M
     frequency doublings before breakdown then

       – A sequence with a linear growth rate of       2  will undergo
         M+1 frequency doublings.