# Full text of "Mathematical Biology and Ecology Lecture Notes"

## See other formats

Mathematical Biology and Ecology Lecture Notes Dr Ruth E. Baker Michaelmas Term 2011 Contents 1 Introduction 5 1.1 References 6 2 Spatially independent models for a single species 7 2.1 Continuous population models for single species 7 2.1.1 Investigating the dynamics 8 2.1.2 Linearising about a stationary point 11 2.1.3 Insect outbreak model 12 2.1.4 Harvesting a single natural population 15 2.2 Discrete population models for a single species 18 2.2.1 Linear stability 20 2.2.2 Further investigation 20 2.2.3 The wider context 25 3 Continuous population models: interacting species 27 3.1 Predator-prey models 27 3.1.1 Finite predation 29 3.2 A look at global behaviour 30 3.2.1 Nullclines 31 3.2.2 The Poincare-Bendixson Theorem 31 3.3 Competitive exclusion 32 3.4 Mutualism (symbiosis) 35 3.5 Interacting discrete models 35 4 Enzyme kinetics 36 4.1 The Law of Mass Action 36 4.2 Michaelis-Menten kinetics 37 4.2.1 Non-dimensionalisation 38 4.2.2 Singular perturbation investigation 38 4.3 More complex systems 40 4.3.1 Several enzyme reactions and the pseudo-steady state hypothesis . . 40 4.3.2 Allosteric enzymes 41 4.3.3 Autocatalysis and activator-inhibitor systems 41 2 CONTENTS 3 5 Introduction to spatial variation 43 5.1 Derivation of the reaction-diffusion equations 44 5.2 Chemotaxis 46 6 Travelling waves 48 6.1 Fisher's equation: an investigation 48 6.1.1 Key points 48 6.1.2 Existence and the phase plane 50 6.1.3 Relation between the travelling wave speed and initial conditions . . 53 6.2 Models of epidemics 54 6.2.1 The SIR model 55 6.2.2 An SIR model with spatial heterogeneity 56 7 Pattern formation 59 7.1 Minimum domains for spatial structure 59 7.1.1 Domain size 60 7.2 Diffusion-driven instability 61 7.2.1 Linear analysis 62 7.3 Detailed study of the conditions for a Turing instability 65 7.3.1 Stability without diffusion 65 7.3.2 Instability with diffusion 66 7.3.3 Summary 67 7.3.4 The threshold of a Turing instability. 68 7.4 Extended example 1 68 7.4.1 The influence of domain size 69 7.5 Extended example 2 69 8 Excitable systems: nerve pulses 72 8.1 Background 72 8.1.1 Resistance 73 8.1.2 Capacitance 73 8.2 Deducing the Fitzhugh Nagumo equations 74 8.2.1 Space-clamped axon 74 8.3 A brief look at the Fitzhugh Nagumo equations 76 8.3.1 The (n, phase plane 76 8.4 Modelling the propagation of nerve signals 78 8.4.1 The cable model 78 A The phase plane 81 A.l Properties of the phase plane portrait 82 A. 2 Equilibrium points 82 A. 2.1 Equilibrium points: further properties 83 A. 3 Summary 84 A. 4 Investigating solutions of the linearised equations 84 A.4.1 Case I 85 CONTENTS 4 A.4.2 Case II 87 A.4.3 Case III 87 A. 5 Linear stability 90 A. 5.1 Technical point 90 A. 6 Summary 91 Chapter 1 Introduction An outline for this course. • We will observe that many phenomena in ecology, biology and biochemistry can be modelled mathematically. • We will initially focus on systems where the spatial variation is not present or, at least, not important. Therefore only the temporal evolution needs to be captured in equations and this typically (but not exclusively) leads to difference equations and/or ordinary differential equations. • We are inevitably confronted with systems of non-linear difference or ordinary dif- ferential equations, and thus we will study analytical techniques for extracting in- formation from such equations. • We will proceed to consider systems where there is explicit spatial variation. Then models of the system must additionally incorporate spatial effects. • In ecological and biological contexts the main physical phenomenon governing the spatial variation is typically, but again not exclusively, diffusion. Thus we are in- variably required to consider parabolic partial differential equations. Mathematical techniques will be developed to study such systems. • These studies will be in the context of ecological, biological and biochemical appli- cations. In particular we will draw examples from: — enzyme dynamics and other biochemical reactions; — epidemics; — interaction ecological populations, such as predator-prey models; — biological pattern formation mechanisms; — chemotaxis; — the propagation of an advantageous gene through a population; — nerve pulses and their propagation. 5 Chapter 1. Introduction 6 1.1 References The main references for this lecture course will be: • J. D. Murray, Mathematical Biology, 3rd edition, Volume I [8]. • J. D. Murray, Mathematical Biology, 3rd edition, Volume II [9]. Other useful references include (but are no means compulsory): • J. P. Keener and J. Sneyd, Mathematical Physiology [7]. • L. Edelstein-Keshet, Mathematical Models in Biology [2]. • N. F. Britton, Essential Mathematical Biology [1]. Chapter 2 Spatially independent models for a single species In this chapter we consider modelhng a single species in cases where spatial variation is not present or is not important. In this case we can simply examine the temporal evolution of the system. References. • J. D. Murray, Mathematical Biology, 3rd edition. Volume I, Chapter 1 and Chapter 2 [8]. • L. Edelstein-Keshet, Mathematical Models in Biology, Chapter 1, Chapter 2 and Chapter 6 [2]. • N. F. Britton, Essential Mathematical Biology, Chapter 1 [1]. 2.1 Continuous population models for single species A core feature of population dynamics models is the conservation of population number, i.e. rate of increase of population = birth rate — death rate (2-1) + rate of immigration — rate of emigration. We will make the assumption the system is closed and thus there is no immigration or emigration. Let N{t) denote the population at time t. Equation (2.1) becomes ^ = fiN) = Ng{N), (2.2) where g{N) is defined to be the intrinsic growth rate. Examples include: 7 Chapter 2. Spatially independent models for a single species 8 The Malthus model. This model can be written as: g{N) = b-d''^^r, (2.3) where b and d are constant birth and death rates. Thus ^ = rN, (2.4) and hence N{t) = Noe^'. (2.5) The Verhulst model. This model is also known as the logistic growth model: giN) =r(l-^Y (2.6) Definition. In the logistic growth equation, r is defined to be the linear birth rate and K is defined to be the carrying capacity. For N <^ K, we have ~ riV ^ AT ~ ATqc^*. (2.7) However, as N tends towards K, dN , , -^0, (2.8) the growth rate tends to zero. We have dN ( N\ , , (2.9) and hence NnKe''^ N{t) = — -f— -^K as t^oo. (2.10) K + Noie^-^-l) ^ ' Sketching N{t) against time yields solution as plotted in Figure 2.1: we see that solutions always monotonically relaxes to as t — t- oo. Aside. The logistic growth model has been observed to give very good fits to popula- tion data in numerous, disparate, scenarios ranging from bacteria and yeast to rats and sheep [8]. 2.1.1 Investigating the dynamics There are two techniques we can use to investigate the model ^ = f{N) = Ng{N). (2.11) Chapter 2. Spatially independent models for a single species 9 120 K 40 300 600 time (t) 900 300 600 time (t) 900 Figure 2.1: Logistic growtli for Nq < K (left-liand) and Nq > K (riglit-liand). Parameters are as follows: r = 0.015 and K ~ 100. Method (i): analytical solution For the initial conditions N{t = 0) = A^O) with A'o fixed, we can we formally integrate equation (2.2) to give N(t) = N*{t), where N*{-) is the inverse of the function F{-) defined by F{x) 1 ■ds. (2.12) However, unless integrating and finding the inverse function is straightforward, there is an easier way to determine the dynamics of the system. Method (ii): plot the graph Plot dN/dt = f{N) = Ng{N) as a function of N. For example, with f{N) = Ng{N) = N{6N'^ - - UN + 6) = N{N - 1){N - 2)(3 - N), (2.13) we have the plot shown in Figure 2.2. t 1 2 population (N) Figure 2.2: Growth according to the dynamics f{N) = N{N - 1){N - 2) (3 - N). Chapter 2. Spatially independent models for a single species 10 Note 1. For a given initial condition, Nq, the system will tend to the nearest root of f{N) = Ng{N) in the direction of /{Nq). The value of |iV(t)| will tend to infinity with large time if no such root exists. For f{N) = Ng{N) = N{N - 1){N - 2)(3 - N), we have: • when A'^o ^ (0,2] the large time asymptote is N{oo) = 1; • for A'^o > 2 the large time asymptote is N(oo) = 3; • N{t) = Vt if N{0) = 0. Note 2. On more than one occasion we will have a choice between using a graphical method and an analytical method, as seen above. The most appropriate method to use is highly dependent on context. The graphical method, Method (ii), quickly and simply gives the large time behaviour of the system and stability information (see below). The analytical method, Method (i), is often significantly more cumbersome, but yields all information, at a detailed quantitative level, about the system. Definition. A stationary point, also known as an equilibrium point, is a point where the dynamics does not change in time. Thus in our specific context of dN/dt = f{N) = Ng{N), the stationary points are the roots of f{N) = 0. Example. For dN/dt = f{N) = Ng{N) = N{N - 1){N - 2)(3 - A^), the stationary points are A^ = 0,1,2,3. (2.14) Definition. A stationary point is stable if a solution starting sufficiently close to the stationary point remains close to the stationary point. Non-examinable. A rigorous definition is as follows. Let A^Ary(t) denote the solution to dN/dt = f{N) = Ng{N) with initial condition N{t = 0) = A^q. A stationary point, A'^s, is stable if, and only if, for all e > there exists a 6 such that if [A'^s — A'ol < 5 then \NN,{t)-Ns\ <e. Exercise. Use Figure 2.2 to deduce which of stationary points of the system dN f{N) = Ng{N) = N{N - 1){N - 2)(3 - A^), (2.15) are stable. Solution. Figure 2.2 shows that both A^^ = 1 and A^^ = 3 are stable. Chapter 2. Spatially independent models for a single species 11 2.1.2 Linearising about a stationary point Suppose Ns is a stationary point of dN/dt = f{N) and make a small perturbation about N{t) = Ns + 7i{t), n{t)<^Ns. (2.16) We have, by using a Taylor expansion of f{N) and denoting ' = d/dN, that f{N{t)) = f{N, + n{t)) = f{N,) + n{t)f'{N,) + ^n{tff"{Ns) + ..., (2.17) and hence ^ = ^ = f{N{t)) = f{Ns) + n{t)f'{N,) + ln{t)'f"{Ns) + ... (2.18) The linearisation of dN/dt = f{N) about the stationary point Ns is given by neglecting higher order (and thus smaller) terms to give ^=m)n(,). The solution to this linear system is simply n{t) = n(t = 0) exp t^iNs) dN^ (2.19) Definition. Let Ns denote a stationary point of dN/dt = f{N), and let n{t) = n{t = 0) exp (2.20) be the solution of the linearisation about Ns- Then Ns is linearly stable if n{t) — ?• as t — 7- oo. In other words, Ns is linearly stable if ^(A^.)<0. (2.21) Exercise. By algebraic means, deduce which stationary points of the system dN — = f{N) = Ng{N) = N{N - l)(iV - 2)(3 - N), (2.22) are linearly stable. Can your answer be deduced graphically? Solution. Differentiating f{N) with respect to gives f{N) = 2- 22N + ISiV^ - 4iV^ (2.23) and hence /'(O) = 6 (unstable), /'(I) = —8 (stable) etc. Consider the graph of f{N) to deduce stability graphically — steady states with negative gradient are linearly stable c.f. Figure 2.2. Chapter 2. Spatially independent models for a single species 12 population (N) Figure 2.3: Growth according to the dynamics f{N) = (1 — 7V)^. Exercise. Find a function f{N) such that dN/dt = f{N) has a stationary point which is stable and not linearly stable. Solution. The function f{N) = {l-Nf, gives /'(I) = and is therefore not linearly stable (see Figure 2.3). (2.24) 2.1.3 Insect outbreak model First introduced by Ludwig in 1978, the model supposes budworm population dynamics to be modelled by the following equation: (2.25) The function p{N) is taken to represent the effect upon the population of predation by birds. Plotting p{N) as a function of N gives the graph shown in Figure 2.4. 250 population (N) 500 Figure 2.4: Predation, p{N), in the insect outbreak model. Parameters are as follows: A = 150, B = 0.5. Chapter 2. Spatially independent models for a single species 13 Non-dimensionsionalisation Let N = N*u, t = Tt, (2.26) where A^*, N have units of biomass, and t, T have units of time, with N* , T constant. Then N*du ( N*u\ B{N*fu^ -Yd^ = "^^H ^^J"^4^T(ivW ^ ^ ^ du ( N*u\ BTN*u'^ ^ -r = rsTu 1 - -— - TTT-^. 2.28 dr V Kb J + {N*Yv? ^ ' Hence with AT* A rj. ^ rp rsA Kb Kb , . N =A T=-, r = rBT = ^, ^ = _ = (2.29) we have du f 1 ^e/ Thus we have reduced the number of parameters in our model from four to two, which substantially simplifies our subsequent study. Steady states The steady states are given by the solutions of \ 2 U\ U ru\\--\- -— ^ = 0. (2.31) V q) 1 + Clearly w = is a steady state. We proceed graphically to consider the other steady states which are given by the intersection of the graphs /i(^z)=r(l--') and f^{u) = (2.32) The top left plot of Figure 2.5 shows plots of fi{u) and /2('u) for different values of r and q. We see that, depending on the values of r and q, we have either one or three non-zero steady states. Noting that d/(u;r, q) = r>0, (2.33) u=0 du typical plots of du/dr vs. u are shown in Figure 2.5 for a range of values of r and q Definition. A system displaying hysteresis exhibits a response to the increase of a driving variable which is not precisely reversed as the driving variable is decreased. Remark. Hysteresis is remarkably common. Examples include ferromagnetism and elasticity, amongst others. See http://en.wikipedia.org/wiki/Hysteresis for more details. Chapter 2. Spatially independent models for a single species 14 0.0 scaled population (u) scaled population (u) scaled population (u) scaled population (u) Figure 2.5: Dynamics of the non-dimensional insect outbreak model. Top left: plots of the functions fi{u) (dashed line) and f2{u) (solid line) with parameters r = 0.2,0.4,0.6, q = 10, 15, 20, respectively. Top right: plot of f{u; r, q) with parameters r = 0.6, q = 0.6. Bottom left: plot of /(w; r, q) with parameters r = 0.6, q = Q. Bottom right: plot of /(u; r, q) with parameters r = 0.6, q = 10. Extended Exercise • Fix r = 0.6. Explain how the large time asymptote of u, and hence N, changes as one slowly increases q from q <^ 1 to q ^ 1 and then one decreases q from g ^ 1 to g <C 1. In particular, show that hysteresis is present. Note for this value of r, there are three non-zero stationary points for q £ {qi, q2) with 1 < qi < q2 < 10. Solution. For small values of q there is only one non-zero steady state, S\. As q is increased past qi, three non-zero steady states exist. Si, S2, S3, but the system stays at Si. As q is increased further, past q2, the upper steady state ^3 is all that remains and hence the system moves to 53. If g is now decreased past 52, three non-zero steady states (5*1, ^2, ^3) exist but the system remains at ^3 until q is decreased past qi. Figure 2.6 shows f{u; r, q) for different values of q. The dashed line shows a plot for q = qi whilst the dash-dotted line shows a plot for q = q2- • What is the biological interpretation of the presence of hysteresis in this model? Chapter 2. Spatially independent models for a single species 15 scaled population (u) Figure 2.6: Left-hand plot: dw/dr = /(u; r, q) in the non-dimensional insect outbreak model as q is varied. For small q there is one, small, steady state, for q G (91,92) there are three non-zero steady states and for large q there is one, large, steady state. Right-hand plot: the steady states plotted as a function of the parameter q reveals the hysteresis loop. Solution. If the carrying capacity, q, is accidentally manipulated such that an out- break occurs (Si —7- 5*3) then reversing this change is not sufficient to reverse the outbreak. 2.1.4 Harvesting a single natural population We wish to consider a simple model for the maximum sustainable yield. Suppose, in the absence of harvesting, we have We consider a perturbation from the non-zero steady state, N = K. Thus we write = K + and find, on linearising, — - = -rn =^ n = noe"""*. (2.35) Hence the system returns to equilibrium on a timescale of 7^(0) = 0{l/r). We consider two cases for harvesting: • constant yield, Y; • constant effort, E. Constant yield For a constant yield, Y = Yq, our equations are !!^ = riv(l-^)-r„^l//(A';r„). (2,36) Chapter 2. Spatially independent models for a single species 16 25 50 75 100 population (N) Figure 2.7 : Dynamics of the constant yield model for Yq = 0.00, 0.15, 0.30. As Yq is increased beyond a critical value the steady states disappear and — > in finite time. Parameters are as follows: K = 100 and r = 0.01. Plotting dN/dt as a function of reveals (see Figure 2.7) that the steady states disappear as Iq is increased beyond a critical value, and then A^ — )• in finite time. The steady states are given by the solutions of rN*-^-Yo = N* = ^ ^^^^ °^ . (2.37) Therefore extinction will occur once rK Yo > — . (2.38) Constant effort For harvesting at constant effort our equations are dN ( N\ def rN"^ ^ = rN[l--yEN f{N; E) = N{r - E) - (2.39) where the yield is Y{E) = EN. The question is: how do we maximise Y{E) such that the stationary state still recovers? The steady states, N* , are such that f{N*;E) = (see Figure 2.8). Thus N*{E)=^^^:^^=(i-^)k, (2.40) and hence Y*{E) = EN*{E) = (^1- KE. (2.41) Thus the maximum yield, and corresponding value of A^*, are given by the value of E such that dY* r rK K -^ = ^ E = -, Y:,,, = —, N*^,, = -. (2.42) Chapter 2. Spatially independent models for a single species 17 Linearising about the stationary state N*{E) we have = N*{E) + n with n + . . . = — (r — E)n + . . . , dn , df{N;E) — ~ fp N ) H — — ■ — - dt ~ ' ^ dN and hence the recovery time is given by N=N* Tr{E) ~ o r-E (2.43) (2.44) Defining the recovery time to be the time for a perturbation to decrease by a factor of e according to the hnearised equations about the non-zero steady state, then rfl(o) = Hence, at the maximum yield state, 1 Tr{E) 1 r-E 2 r Tr(E) = — since E = — at maximum yield. r 2 (2.45) (2.46) As we measure Y it is useful to rewrite E in terms of Y to give the ratio of recovery times in terms of the yield Y{E) and the maximum yield Ym- Tb{Y) Tr{0) Derivation. At steady state, we have K (2.47) E^ - KE + Y* = as Y* = EN* =KE{1 r \ r E (2.48) This gives E ±r^l-AY*/K7 r-E 2 2 Substituting into equation (2.45) gives the required result itWi (2.49) Plotting Tji(Y)/Tr{0) as a function of Y/Ym yields some interesting observations, as shown in Figure 2.8. Note. As Tr increases the population recovers less quickly, and therefore spends more time away from the steady state, N* . The biological implication is that, in order to maintain a constant yield, E must be increased. This, in turn, implies Tr increases, resulting in a positive feedback loop that can have disastrous consequences upon the population. Chapter 2. Spatially independent models for a single species 18 0.3 0.2 0.1 0.0 /' / y // \ / / / ' A // ^ / / // / 15 10 ■ \ ■ \ \ \ \ \ \ — . 25 50 75 population (N) 100 0.0 0.2 0.4 0.6 0.8 Y/ Y. 1.0 Figure 2.8: Dynamics of the constant effort model. The Icft-liand plot shows the logistic growth curve (solid line) and the yield, Y = EN (dashed lines), for two values of E. The right-hand plot shows the ratio of recovery times, Tii{Y)/Tii{0), with the negative root plotted as a dashed line and the positive root as a solid line. Parameters are as follows: K = 100 and r = 0.01. 2.2 Discrete population models for a single species When there is no overlap in population numbers between each generation, we have a discrete model: Nt+i = Ntf{Nt) = H{Nt). (2.50) A simple example is Nt+i = rNt, (2.51) which implies oo r > 1 Nt = r*iVo ^ "! No r = l . (2.52) r < 1 Definition. An equilibrium point, N*, for a discrete population model satisfies N* = N*f{N*) = H{N*). (2.53) Such a point is often known fixed point. An extension of the simple model, equation (2.51), called the Ricker model includes a reduction of the growth rate for large Nt: Nt' ' Nt+i = Nt exp or, in non-dimensionalised form, r 1 K , r > K>0, def ut+i=utex-p[r{l-ut)] = H{ut). (2.54) (2.55) We can start developing an idea of how this system evolves in time via cobwebbing, a graphical technique, as shown in Figure 2.9. Chapter 2. Spatially independent models for a single species 19 50 100 150 200 2 4 6 8 10 ^\ generation (t) Figure 2.9: Dynamics of tlie Ricker model. The left-hand plot shows a plot of Nt+i = Nt exp [r (1 — Nt/K)] alongside iVt+i = Nt with the cobwebbing technique shown. The right- hand plot shows Nt for successive generation times t = 1, 2, . . . , 10. Parameters are as follows: iVo = 5, r = 1.5 and K = 100. In particular, it is clear that the behaviour sufficiently close to a fixed point, u* , depends on the value of H'{u*). For example: • -1< H'{u*) < t • H'{u*) = -1 t Chapter 2. Spatially independent models for a single species 20 • H'{u*) < -1 t 2.2.1 Linear stability More generally, to consider the stability of an equilibrium point algebraically, rather than graphically, we write ut = u* + vt, (2.56) where u* is an equilibrium value. Note that u* is time-independent and satisfies u* = H{u*). Hence ut+i = u* + vt+i = H{u* + Vt) = H{u*) + vtH\u*) + o{v^). (2.57) Consequently, we have vt+i = H'{u*)vt where H'{u*) is a constant, independent of t, (2.58) and thus vt=[H'{u*)\\Q. (2.59) This in turn enforces stability if \H'{u*)\ < 1 and instability if \H'{u*)\ > 1. Definition. A discrete population model is linearly stable if \H'{u*)\ < 1. 2.2.2 Further investigation The equations are not as simple as they seem. For example, from what we have seen thus far, the discrete time logistic model seems innocuous enough. Nt+i =rNt(^l-^^, r > K > 0. (2.60) If we put in enough effort, one could be forgiven for thinking that the use of cobwebbing will give a simple representation of solutions of this equation. However, the effects of increasing r are stunning. Figure 2.10 shows examples of cobwebbing when r = 1.5 and r = 4.0. It should now be clear that even this simple equation does not always yield a simple solution! How do we investigate such a complicated system in more detail? Chapter 2. Spatially independent models for a single species 21 Definition. A bifurcation point is, in the current context, a point in parameter space where the number of equilibrium points, or their stability properties, or both, change. We proceed to take a closer look at the non-dimensional discrete logistic growth model: ut+i = rut{l - ut) = H{ut), (2.61) for different values of the parameter r, and, in particular, we seek the values where the number or stability nature of the equilibrium points change. Note that we have equilibrium points at M* = and u* = (r — l)/r, and that H'{u) = r — 2ru. For < r < 1, we have: • = is a stable steady state since |-ff'(0)| = \r\ < 1; • the equilibrium point at u* = (r — l)/r is unstable. It is also unreachable, and thus irrelevant, for physical initial conditions with no > 0. For 1 < r < 3 we have: • u* = is an unstable steady state since |ff'(0)| = |r| > 1; • u* = {r — l)/r is an stable steady state since |if'((r — l)/r)| = |2 — r| < 1. In Figure 2.11 we plot this on a diagram of steady states, as a function of r, with stable steady states indicated by solid lines and unstable steady states by dashed lines. When r = 1 we have (r — l)/r = 0, so both equilibrium points are at u* = 0, with H'[u* = 0) = 1. Clearly we have a switch in the stability properties of the equilibrium points, and thus r = 1 is a bifurcation point. Chapter 2. Spatially independent models for a single species 22 0.8 0.8 0.4 0.2 0.0 ' 12 3 r Figure 2.11: Bifurcation diagram for the non-dimensional discrete logistic model. The non- zero steady state is given, for r > 1, by N* = {r — l)/r. What happens for r > 3? We have equihbrium points at n* = 0, u = (r — l)/r and H'{u* = (r — l)/r) < —1; both equilibrium points are unstable. Hence if the system approaches one of these equilibrium points the approach is only transient; it quickly moves away. We have a switch in the stability properties of the equilibrium points, and thus r = 3 is a bifurcation point. To consider the dynamics of this system once r > 3 we consider Ut+2 = H{ut+i) = H [H{ut)] =^ H2{ut) = r [rut{l - Uf)] [1 - rut{l - ut)] . (2.62) Figure 2.12 shows H2{ut) for r = 2.5 and r = 3.5 and demonstrates the additional steady states that arise as r is increased past r = 3. Figure 2.12: Dynamics of the non-dimensional discrete logistic model in terms of every second iteration. The left-hand plot shows results for r = 2.5 whilst the right-hand plot shows results for r = 3.5. Note. The fixed points of H2 satisfy U2 = H2{u2), which is a quartic equation in However, we know two solutions, the fixed points H{-), i.e. and (r — l)/r. Using standard Chapter 2. Spatially independent models for a single species 23 techniques we can reduce the quartic to a quadratic, which can be solved to reveal the further fixed points of H2, namely ul = '^±^[{r-lf-A]"\ (2.63) Zr Zr ^ J These roots exist if (r — 1)^ > 4, i.e. r > 3. Definition. The m}^ composition of the function H is given by Hm{u)'^= [H ■ H ...H ■ H]{u). (2.64) ^ V ' m times Definition. A point u is periodic of period Hm{u) = u, Hi{u) / tt, m for the function H if i € {1,2, ... m - 1}. (2.65) Thus the points u^ = !:±l±l[(r-l)^-4]^/^ (2.66) are points of period 2 for the function H. Problem. Show that the U2 are stable with respect to the function H2 for r > 3, (r-3) < 1. Let ^^dgr + l^l 2_ 1,2^ ni = F(uo), U2 = H2{u,), (2.67) Zr Zr ^ J and let d_ du Then \ = — [H2{u)]\u=u,. (2.68) \ = ^[H- H{u)] \u=uo = H'{uo)H'{ui). (2.69) Thus for stability we require \H'{uq)H'{ui)\ < 1. Exercise. Finish the problem: show that the steady states are stable for the dynamical system ut+i = H2{ut), with r > 3, (r — 3) <C 1. Chapter 2. Spatially independent models for a single species 24 Exercise. Suppose uq is an equilibrium point of period m for the function H. Show that uq is stable if ni^o^ [H'{u,)] < 1, (2.71) where Ui = Hi{uo) for i G {1, 2, . . . , m — 1}. Solution. Defining A in a similar manner as before, we have d A du d_ du u=uo [H{Q{u)] where Q{u) = Hm-i{u), U=Uq d = H (itm-l)^-f^m-l Hence, by iteration, we have the result. (2.72) (2.73) (2.74) (2.75) u=uo We plot the fixed points of H2, which we now know to be stable, in addition to the fixed points of Hi, in Figure 2.13. The upper branch, U2u, is given by the positive root of equation (2.70) whilst the lower branch, U21, is given by the negative root. We have = H{u2u), u^u = H{u2i). Thus a stable, period 2, oscillation is present, at least for (r — 3) ^ 1. Any solution which gets sufficiently close to either U2u or ut^^^ stays close. Figure 2.13: Bifurcation diagram for the non-dimensional discrete logistic model with inclu- sion of the period 2 solutions. For higher values of r, there is a bifurcation point for H2; we can then find a stable fixed point for H/i[u) : : H2[H2{u)] in a similar manner. Increasing r further there is a bifurcation point for H4^(u). Again, we are encountering a level of complexity which is too much to deal with our current method. To bring further understanding to this complex system, we note the following definition. Chapter 2. Spatially independent models for a single species 25 Definition. An orbit generated by the point uq are the points {uo,ui,U2, ,U3, . . .} where Ui = Hi{u) = H{ui-i). We are primarily interested in the large time behaviour of these systems in the context of biological applications. Thus, for a fixed value of r, we start with a reasonable initial seed, say u* = 0.5, and plot the large time asymptote of the orbit of m*, ie. the points Hi{u*) once i is sufficiently large for there to be no transients. This gives an intriguing plot; see Figure 2.14. 1.0 0.8 0.6 0.4 0.2 0.0 3.0 3.2 3.4 3.6 3.8 4.0 r Figure 2.14: The orbit diagram of the logistic map. For each value of r € [3,4] along horizontal axis, points on the large time orbits of the logistic map arc plotted. In particular, we have regions where, for r fixed, there are three points along the ver- tical corresponding to period 3 oscillations. This means any period of oscillation exists and we have a chaotic system. This can be proved using Sharkovskii's theorem. See P. Glendinning, Stability, Instability and Chaos [4] for more details on chaos and chaotic systems. Note. A common discrete population model in mathematical biology is Nt+i = * , . (2.76) Models of this form for the Colorado beetle are within the periodic regimes, while Nichol- son's blowfly model is in the chaotic regime [8]. 2.2.3 The wider context In investigating the system ut+i = rut (1 - ut) H{ut), (2.77) we have explored a very simple equation which, in general, exhibits greatly different be- haviours with only a small change in initial conditions or parameters {i.e. linear growth rate, r). Such sensitivity is a hallmark feature of chaotic dynamics. In particular, it makes Chapter 2. Spatially independent models for a single species 26 prediction very difficult. There will always be errors in a model's formalism, initial condi- tions and parameters and, in general, there is no readily discernible pattern in the way the model behaves. Thus, assuming the real system behaviour is also chaotic, using statistical techniques to extract a pattern of behaviour to thus enable an extrapolation to predict future behaviour is also fraught with difficulty. Attempting to make accurate predictions with models containing chaos is an active area of research, as is developing techniques to analyse seemingly random data to see if such data can be explained by a simple chaotic dynamical system. Chapter 3 Continuous population models: interacting species In this chapter we consider interacting populations, but again in the case where spatial variation is not important. Appendix A contains relevant information for phase plane analysis that may be useful. References. • J. D. Murray, Mathematical Biology, 3rd edition, Volume I, Chapter 3 [8]. • L. Edelstein-Keshet, Mathematical Models in Biology, Chapter 6 [2]. • N. F. Britton, Essential Mathematical Biology, Chapter 2 [1]. There are three main forms of interaction: Predator-prey An upsurge in population I (prey) induces a growth in population II (predator). An upsurge in population II (predator) induces a decline in population I (prey). Competition An upsurge in either population induces a decline in the other. Symbiosis An upsurge in either population induces an increase in the other. Of course, there are other possible interactions, such as cannibalism, especially with the "adult" of a species preying on the young, and parasitism. 3.1 Predator-prey models The most common predator-prey model is the Lotka-Volterra model. With the number of prey and P the number of predators, this model can be written dA^ dP dF cNP - dP, aN - bNP, (3.1) (3.2) 27 Chapter 3. Continuous population models: interacting species 28 with a, b, c, d positive parameters and c < b. Non-dimensionalising with u = {c/d)N, V = {b/a)P, T = at and a = d/a, we have 1 ddu ad bda du ^ rr x /o on — ; — = U —UV =^ —— = U — UV = U(l — V) = t(U,V), (o.o) l/acdT c c b dr ^ j j \ ^ \ j 1 adv da a dv — ■j- — — = c——uv — d—v, =^ — = a\uv — V) = aviu — 1) = q{u,v), 3.4 l/abdT cb b dr ^ ' ^ i ^ \ j There are stationary points at {u,v) = (0,0) and {u,v) = (1, 1). Exercise. Find the stabihty of the stationary points {u,v) = (0,0) and {u,v) = (1,1). The Jacobian, J, is given by fu fv \ ^ I 'i^-v -u Qu Qv I \ av a{u - 1) (3.5) At (0, 0) we have -a with eigenvalues 1, —a. Therefore the steady state (0,0) is an unstable saddle. At (1, 1) we have -1 a (3.7) with eigenvalues ziz^^/a. Therefore the steady state (1, 1) is a centre (not linearly stable). These equations are special; we can integrate them once, as follows, to find a conserved constant: du u(l — v) /"u— 1, f 1 — V , . A- = -r — rf ^ / / — • dv a{u — l)v J u J av Hence H = const = au + V — alnu — Inv. (3.9) This can be rewritten as / pV\ / pU\a^ ' ^ ^ e^, (3.10) V / \ u from which we can rapidly deduce that the trajectories in the (n, v) plane take the form shown in Figure 3.1. Thus u and v oscillate in time, though not in phase, and hence we have a prediction; predators and prey population numbers oscillate out of phase. There are often observations of this e.g. hare-lynx interactions. Chapter 3. Continuous population models: interacting species 29 0.0 1.0 2.0 3.0 4.0 5.0 4 8 12 16 u Time Figure 3.1: Dynamics of tlie non-dimensional Lotka-Volterra system for a = 1.095 and H — 2.1,2.4,3.0,4.0. The left-hand plot shows the dynamics in the {u,v) phase plane whilst the right-hand plot shows the temporal evolution of u and v. 3.1.1 Finite predation The common predator-prey model assumes that as — ?• oo the rate of predation per predator becomes unbounded, as does the rate of increase of the predator's population. However, with an abundance of food, these quantities will saturate rather than become unbounded. Thus, a more realistic incorporation of an abundance of prey requires a refine- ment of the Lotka-Volterra model. A suitable, simple, model for predator-prey interactions under such circumstances would be (after a non-dimensionalisation) du „, s. , . auv — = f{u,v) = u{l-u)--—, 3.11 dr a + u g{u,v)=bv(l--), (3.12) dr \ uy where a, b, d are positive constants. Note the effect of predation per predator saturates at high levels of u whereas the predator levels are finite at large levels of prey and drop exceedingly rapidly in the absence of prey. There is one non-trivial equilibrium point, (u*,?;*), satisfying OjU* V* = V* where il — u*) = , (3.13) and hence * u = - 2 The Jacobian at [u* ,v* ) is -(a + d - 1) + V(a + d - 1)2 + 4d . (3.14) ( fu fu \ \ gu gv J (3.15) {u* ,v*) where fu{u ,v ) = l-2u - — — - + =-u + — — -2 . (3.16 d + u* [d + u*)'^ [d + u*)'' Chapter 3. Continuous population models: interacting species 30 b(v*)^ 9u{u*y) = ^^ = b, (3.18) g,{u\v*) = b[l-2^—j =-h. (3.19) The eigenvalues satisfy (A-/n)(A-5.)- A'5« = X^-{fu + 9v)X + {fu9v-fvgu) = 0, (3.20) and hence A^-aA + /3 = =, A= "^^f (3.21) where Note that 0=1- «(^*)' = 1 _ = (n*+d)-n* + (u*)^ ^ d+(n*)^ > .3 (n* + d)2 {u*+d) u*+d d + u* ■ > Thus, if a < the eigenvalues A are such that we have either: • a stable node (a^ - 4/3 > 0); • stable focus (q^ - 4/3 < 0); at the equilibrium point {u*,v*). If a > we have an unstable equilibrium point at (n* , f * ) . 3.2 A look at global behaviour This previous section illustrated local dynamics: we have conditions for when the dynamics will stably remain close to the non-trivial equilibrium point. One is also often interested in the global dynamics. However, determining the global dynamics of a system, away from its equilibrium points, is a much more difficult problem compared to ascertaining the local dynamics, sufficiently close to the equilibrium points. For specific parameter values, one can readily solve the ordinary differential equations to consider the behaviour of the system. One is also interested in the general properties of the global behaviour. This is more difficult, and we will consider one possible approach below. There are many potential tools available: nullcline analysis, the Poincare-Bendixson The- orem, the Poincare Index and the Bendixson-Dulac Criterion. The Poincare-Bendixson Theorem is a useful tool for proving that limit cycles must exist, while Poincare indices and the Bendixson-Dulac Criterion are useful tools for proving a limit cycle cannot exist. Chapter 3. Continuous population models: interacting species 31 We will briefly consider nullclines and the Poincare-Bendixson Theorem in detail. Please refer to P. Glendinning, Stability, Instability and Chaos: An Introduction to the Theory of Nonlinear Differential Equations [4], or D. W. Jordan and P. Smith, Mathematical Techniques: An Introduction for Engineering, Mathematical and Physical Sciences [6], for further details than considered here. 3.2.1 Nullclines Definition. Consider the equations ^ = f{u,v), f^=9{n,v). (3.24) The nullclines are the curves in the phase plane where f{u, v) = and g{u, v) = 0. Reconsider du X / N auu — = f{u,u) = u{l-u)- -—, 3.25 dr a + u dv d7 g{u,v) =bv(l-^y (3.26) The u nullclines are given by f(u,v) = u = and v = -(1 - u)(u + d). (3.27) a The V nullclines are given by g{u,v) = v = and v = u. (3.28) A sketch of the nullclines and the behaviour of the phase plane trajectories is shown in Figure 3.2. 3.2.2 The Poincare-Bendixson Theorem For a system of two first order ordinary differential equations, consider a closed bounded region D. Suppose a positive half path, H, lies entirely within D. Then one of the following is true: 1. H is a closed trajectory, e.g. a limit cycle; 2. H asymptotically tends to a closed trajectory, e.g. a limit cycle; 3. H terminates on a stationary point. Therefore, if D does not have a stationary point then there must be a limit cycle. For a proof see P. Glendinning, Stability, Instability and Chaos: An Introduction to the Theory of Nonlinear Differential Equations [4] . Chapter 3. Continuous population models: interacting species 32 1.0 0.8 0.6 > 0.4 0.2 0.0 0.0 0.2 0.4 0.6 0.8 1.0 u Figure 3.2: The (u, v) phase-plane for the finite predation model when the steady state is stable. The u nullclines are plotted in red and the v nuUclines in green. Trajectories for a number of different initial conditions are shown as dashed lines. Parameters are as follows: a = 2.0, b = 0.1, d 2.0. Exercise. Explain why a > in the previous example (see equation (3.22)) implies we have limit cycle dynamics. What does this mean in terms of the population levels of predator and prey? Solution. For a > the steady state is an unstable node or spiral. Further, we can find a simple, closed boundary curve, C, in the positive quadrant of the {u,v) plane, such that on C phase trajectories always point into the domain, T>, enclosed by C. Applying the Poincare-Benedixon Theorem to the domain gives the existence of a limit cycle. See J. D. Murray, Mathematical Biology Volume I [8] (Chapter 3.4) for more details. 3.3 Competitive exclusion We consider an ordinary differential equation model of two competitors. An example might be populations of red squirrels and grey squirrels [8]. Here, both populations compete for the same resources and a typical model for their dynamics is and N2 with grey squirrels in our example. In particular, given a range of parameter values and some initial values for A'^i and at the time t = 0, we would typically like to know if the final outcome is one of the following possibilities: (3.29) (3.30) where Ki, K2, ri, r2, 612, ^21 are positive constants. Let us associate A''i with red squirrels • the reds become extinct, leaving the greys; Chapter 3. Continuous population models: interacting species 33 • the greys become extinct, leaving the reds; • both reds and greys become extinct; • the reds and greys co-exist. If this system is perturbed in any way will the reds and greys continue to coexist? After a non-dimensionalisation (exercise) we have u'l = ui{l-ui-ai2U2)'^= fi{ui,U2), (3.31) U2 = PU2{1-U2- a2lUi)'^= f2{ui,U2), (3.32) where p = r2/ri. The stationary states are (nt,n^) = (0,0), «,^) = (1,0), (nt,^) = (0,1), (3.33) and ^ {ul,u*2) = - (1 - ai2, 1 - 021), (3.34) 1 — ai2Q;2i if ai2 < 1 and 021 < 1 or Q12 > 1 and a2i > 1- The Jacobian is 1 - 2ui — Qi2ti2 — Qi2ttl -pa2iU2 p{l - 2u2 - a2iui) (3.35) It is a straightforward application of phase plane techniques to investigate the nature of these equilibrium points: Steady state (nj,n2) = (0,0). (' p- X J - AI = ( ^ ^ ^ \ X = l,p. (3.36) Therefore (0, 0) is an unstable node. Steady state (nj,n2) = (1,0). J-AI=("\"^ , ,) A = -l,Kl-a2i). (3.37) y p{i- a2i) -Ay Therefore (1,0) is a stable node if a2i > 1 and a saddle point if a2i < 1. Steady state (n^,n2) = (0, 1). J - AI = I ^ " " ^ ° I A = -p, 1 - ai2. (3.38) \^ -yoa2i -p - X J Therefore (0, 1) is a stable node if ai2 > 1 and a saddle point if ai2 < 1. Chapter 3. Continuous population models: interacting species 34 Steady state K,n^) = Tr^7^(l - ai2,l - 021] J-AI 1 / "21 — 1 — A ai2(ai2 — 1) 1 - 012^21 y pa2i(a2i - 1) p{a2i - 1) - A Stability depends on Q12 and Q21. (3.39) There are several different possible behaviours. The totality of all behaviours of the above model is reflected in how one can arrange the nullclines within the positive quadrant. However, for competing populations these straight line nullclines have negative gradients. Figure 3.3 shows the model behaviour for different sets of parameter values. 0.0 0.4 1.2 0.8 0.4 0.0 / / / / / / ^ / / ' / y / / ^ " " \Nv — 0.0 0.4 1.2 Figure 3.3: Dynamics of the non-dimensional competitive exclusion system. Top left: Q12 = 0.8 < 1, a2\ — 1.2 > 1 and is excluded. Top right: a\2 — 1.2 > 1, 0,21 — 0.8 < 1 and u\ is excluded. Bottom left: ayi = 1.2 > 1, ol2\ = 1-2 > 1 and exclusion is dependent on the initial conditions. Bottom right: oi\2 = 0.8 < 1, ol2\ = 0.8 < 1 and we have coexistence. The stable steady states are marked with *'s and p = 1.0 in all cases. The red lines indicate /i = whilst the green lines indicate /2 = 0. Note. In ecology the concept of competitive exclusion is that two species competing for exactly the same resources cannot stably coexist. One of the two competitors will always have an ever so slight advantage over the other that leads to extinction of the second competitor in the long run (or evolution into distinct ecological niches). Chapter 3. Continuous population models: interacting species 35 3.4 Mutualism (symbiosis) We consider the same ordinary differential equation model for two competitors, i.e. dNi diV2 / N2 ^, iVi' (3.40) (3.41) where Ki, K2, ri, r2, 612, &21 are positive constants or, after non-dimensionalisation, def Ui{l - Ui + ai2U2) = fliui,U2), def PU2{1 - U2 + a2lUi] f2{ui,U2). (3.42) (3.43) In symbiosis, the straight line nullclines will have positive gradients leading to the following two possible behaviours shown in Figure 3.4. Figure 3.4: Dynamics of the non-dimensional symbiotic system. The left-hand figure shows population explosion (ai2 = 0.6 = 0:21) whilst the right-hand figure shows population coex- istence (ai2 = 0.1 — a2i). The stable steady states are marked with *'s and p = 1.0 in all cases. The red lines indicate /i = whilst the green lines indicate /2 = 0. 3.5 Interacting discrete models It is also possible, and sometimes useful, to consider interacting discrete models which take the form ut+i = f{ut,vt), (3.44) vt+i = g{ut,vt), (3.45) and possess steady states at the solutions of u* = f{u*,v*), v* = g{u*,v*). (3.46) It is interesting and relevant to study the linear stability of these equilibrium points, and the global dynamics, but we do not have time to pursue this here. Chapter 4 Enzyme kinetics In this chapter we consider enzyme kinetics, which can be thought of as a particular case of an interacting species model. In all cases here we will neglect spatial variation. Throughout, we will consider the m chemical species Ci, . . . , Cm- • The concentration of Cj, denoted q, is defined to be the number of molecules of Ci per unit volume. • A standard unit of concentration is moles m~^, often abbreviated to mol m~^. Recall that 1 mole = 6.023 x 10^'^ molecules. References. • J. D. Murray, Mathematical Biology, 3rd edition, Volume I, Chapter 6 [8]. • J. P. Keener and J. Sneyd, Mathematical Physiology, Chapter 1 [7]. 4.1 The Law of Mass Action Suppose Ci, . . . , Cm undergo the reaction XlCi + X2C2 + . . . + XmCm V UiCi + U2C2 + . . . + l^mCm- (4.1) h The Law of Mass Action states that the forward reaction proceeds at rate kj4^4\..ct, (4.2) while the back reaction proceeds at the rate kf,C^^ ■ ■ ■ Cm ^ (4-3) where kj and kh are dimensional constants that must be determined empirically. 36 Chapter 4. Enzyme kinetics 37 Note 1. Strictly, to treat kj, kb above as constant, we have to assume that the tem- perature is constant. This is a very good approximation for most biochemical reactions occurring in, for example, physiological systems. However, if one wanted to model re- actions that produce extensive heat for example, burning petrol, one must include the temperature dependence in /cj and kf, and subsequently keep track of how hot the sys- tem gets as the reaction proceeds. This generally makes the modelling significantly more difficult. Below we assume that we are dealing with systems where the temperature is approximately constant as the reaction proceeds. Note 2. The Law of Mass Action for chemical reactions can be derived from statistical mechanics under quite general conditions (see for example L. E. Riechl, A Modern Course in Statistical Physics [11]). Note 3. As we will see later, the Law of Mass Action is also used in biological scenarios to write down equations describing, for example, the interactions of people infected with, and people susceptible to, a pathogen during an epidemic. However, in such circumstances its validity must be taken as an assumption of the modelling; in such scenarios one cannot rely on thermodynamic/statistical mechanical arguments to justify the Law of Mass Action. 4.2 Michaelis-Menten kinetics Michaelis-Menten kinetics approximately describe the dynamics of a number of enzyme systems. The reactions are ki S + E . SE, (4.4) k-i k2 SE ^ P + E. (4.5) Letting c denoting the concentration of the complex SE, and s, e, p denoting the con- centrations of S, E, P, respectively, we have, from the Law of Mass Action, the following ordinary differential equations: ds — = —kise + k^ic, (4.6) dc — = kise — k^ic — k2C, (4.7) de — = —kise + k-ic + k2C, (4-8) I = k,c. (4.9) Note that the equation for p decouples and hence we can neglect it initially. The initial conditions are: s{0) = so, e(0) = eo < sq, c(0) = 0, p(0) = 0. (4.10) Chapter 4. Enzyme kinetics 38 Key Point. In systems described by the Law of Mass Action, linear combinations of the variables are often conserved. In this example we have A(e + c) = ^ e = eo-c, (4.11) at and hence the equations simplify to: ds — = -ki{eo - c)s + k_ic, (4.12) dc — = ki{eo - c)s - {k_i + k2)c, (4.13) with the determination of p readily achievable once we have the dynamics of s and c. 4.2.1 Non-dimensionalisation We non-dimensionalise as follows: , , S C k2 def eo def k^^rk2 , . T = kieot, u = — , V = — , A = , e = — < 1, K = — ; , (4.14) So eo fciSQ So kiso which yields u' = -u + {u + K - X)v, (4.15) ev' = u-{u + K)v, (4.16) where u(0) = 1, v{0) = and e ^ 1. Normally e ~ 10 ^. Setting e = yields u (4.17) u + K' which is inconsistent with the initial conditions. Thus we have a singular perturbation problem; there must be a (boundary) region with respect to the time variable around t = where v' oo 0{\). Indeed for the initial conditions given we find v'(0) ~ 0(l/e), with m(0), f (0) < 0(1). This gives us the scaling we need for a singular perturbation investigation. 4.2.2 Singular perturbation investigation We consider a = ^, (4.18) with u(t, e) = u{(T,e) = uo{(7) + eui{(T) + . . . , (4.19) v{t,€) = v{a,e)=VQ{a) + evi{a) + .... (4.20) Proceeding in the usual way, we find that uq, vq satisfy — = =^ uq = constant = 1, (4.21) dcr and ^=uo-{l + K)vo = l-{l + K)vo vo = -— , (4.22) dcr 1 + K Chapter 4. Enzyme kinetics 39 which gives us the inner solution. To find the outer solution we expand u(r,e) = v{T,e) = within the equations to find that and This gives no(T) + eni(T) + . . . , Voir) + evi{T) + . . . , ev dug dr = -u+{u + K-X)v, = u — {u + K)v, = -no + {uq + K - X)vo, = no - (no + i^)no. no dno Ano — and — — = uq + K dr no + i^ In order to match the solutions as cr — )• oo and r — )• we require lim no = hm no = 1 and lim no = 1™ vq = — . (T— ^OO T— >0 IT— >0O T— >-0 1 -\- K (4.23) (4.24) (4.25) (4.26) (4.27) (4.28) (4.29) (4.30) Thus the solution looks like that shown in Figure 4.1. 0.1 0.2 time (t) \ \ \ \ \ \- \ 0.0 0.5 1.0 time (t) 1.5 Figure 4.1: Numerical solution of the non-dimensional Michaelis-Menten equations clearly illustrating the two different time scales. The u dynamics arc indicated by the solid line and the V dynamics by the dashed line. Parameters are e = 0.01, K = 0.1 and A = 1.0. Chapter 4. Enzyme kinetics 40 Often the initial, fast, transient is not seen or modelled and one considers just the outer equations with a suitably adjusted initial condition (ultimately determined from consis- tency/matching with the inner solution). Thus one often uses Michaelis-Menten kinetics where the equations are simply: du Xu , . N , u , . — = with u(0) = 1 and v = -. (4.31) dt u + K ^ ^ u + K ^ ^ Definition. We have, approximately, that dv/dr ~ using Michaelis-Menten kinet- ics. Taking the temporal dynamics to be trivial, ^ 0, (4.32) when the time derivative is fast, i.e. of the form e^=giu,v), (4.33) where e ^ 1, g{u,v) ~ 0{1), is known as the pseudo-steady state hypothesis and is a common assumption in the literature. We have seen its validity in the case of enzyme kinetics about at least away from the inner region. Note. One must remember that the Michaelis-Menten kinetics derived above are a very useful approximation, but that they hinge on the validity of the Law of Mass Action. Even in simple biological systems the Law of Mass Action may breakdown. One (of many) reasons, and one that is potentially relevant at the sub-cellular level, is that the system in question has too few reactant molecules to justify the statistical mechanical assumptions underlying the Lass of Mass Action. Another reason is that the reactants are not well-mixed, but vary spatially as well as temporally. We will see what happens in this case later in the course. 4.3 More complex systems Here we consider a number of other simple systems involving enzymatic reactions. In each case the Law of Mass Action is used to write down a system of ordinary differential equations describing the dynamics of the various reactants. See J. Keener and J. Sneyd, Mathematical Physiology [7], for more details. 4.3.1 Several enzyme reactions and the pseudo-steady state hypothesis We can have multiple enzymes. In general the system of equations reduces to = f{u,Vi,...,Vn), = gi{u,Vi,...,Vn), (4.34) (4.35) Chapter 4. Enzyme kinetics 41 for i G {1, . . . , n}, while the pseudo-steady state hypothesis gives a single ordinary differ- ential equation u' = f{u,vi{u), . . . ,Vn{u)), (4.36) where vi{u), . . . , Vn{u) are the appropriate roots of the equations gi{u,vi, . . . ,Vn) = 0, i£{l,...,n}. (4.37) 4.3.2 Allosteric enzymes Here the binding of one substrate molecule at one site affects the binding of another substrate molecules at other sites. A typical reaction scheme is: S + E . Ci > P + E (4.38) S + Ci ^— ^ C2 > Ci + E. (4.39) k-3 Further details on the investigation of such systems can be found in J. D. Murray, Mathe- matical Biology Volume I [8], and J. P. Keener and J. Sneyd, Mathematical Physiology [7]. 4.3.3 Autocatalysis and activator-inhibitor systems Here a molecule catalyses its own production. The simplest example is the reaction scheme A + b\2B, (4.40) though of course the positive feedback in autocatalysis is usually ameliorated by inhibition from another molecule. This leads to an example of an activator-inhibitor system which have a very rich behaviour. Other examples of these systems are given below. Example 1 This model qualitatively incorporates activation and inhibition: — = cu, (4.41 dt b + v ' ^ ' dv — = du-ev. (4.42) Example 2 This model is commonly referred to as the Gierer-Meinhardt model [3]: du , , , — = a-bu^ , (4.43) at V — = u^-v. (4.44) at Chapter 4. Enzyme kinetics 42 Example 3 This model is commonly referred to as the Thomas model [8]. Proposed in 1975, it is an empirical model based on a specific reaction involving uric acid and oxygen: — = a - u - pR{u,v), (4.45) dv — = a{b-v) - pR{u,v), (4.46) where represents the interactive uptake. Chapter 5 Introduction to spatial variation We have initially considered biological, biochemical and ecological phenomena with neg- ligible spatial variation. This is, however, often not the case. Consider a biochemical reaction as an example. Suppose this reaction is occurring among solutes in a relatively large, unstirred solution. Then the dynamics of the system is not only governed by the dy- namics of the rate at which the biochemical react, but also by the fact there can be spatial variation in solute concentrations, which entails that diffusion of the reactants can occur. Thus modelling such a system requires taking into account both reaction and diffusion. We have a similar problem for population and ecological models when we wish to incor- porate the tendency of a species to spread into a region it has not previously populated. Key examples include modelling ecological invasions, where one species invades another's territory (as with grey and red squirrels in the UK [10]), or modelling the spread of dis- ease. In some, though by no means all, of these ecological and disease-spread models the appropriate transport mechanism is again diffusion, once more requiring that we model both reaction and diffusion in a spatially varying system. In addition, motile cells can move in response to external influences, such as chemical concentrations, light, mechanical stress and electric fields, among others. Of particular interest is modelling when motile cells respond to gradients in chemical concentrations, a process known as chemotaxis, and we will also consider this scenario. Thus, in the following chapters, we will study how to model such phenomena and how (when possible) to solve the resulting equations in detail, for various models motivated from biology, biochemistry and ecology. References. • J. D. Murray, Mathematical Biology Volume I, Chapter 11 [8]. • N. F. Britton, Essential Mathematical Biology, Chapter 5 [1]. 43 Chapter 5. Introduction to spatial variation 44 5.1 Derivation of the reaction-diffusion equations Let i S {1, . . . ,m}. Suppose the chemical species Cj, of concentration q, is undergoing a reaction such that, in the absence of diffusion, one has dcj = Ri{ci,C2,...,Cm)- (5.1) Recah that Ri{ci,C2, ■ ■ ■ ,Cm) is the total rate of production/destruction of Cj per unit volume, i.e. it is the rate of change of the concentration q. Let t denote time, and x denote the position vector of a point in space. We define • c{x,t) to be the concentration of (say) a chemical (typically measured in mol m~'^). • q{x,t) to be the flux of the same chemical (typically measured in mol m~^ s~^)- Recall that the flux of a chemical is defined to be such that, for a given infinitesimal surface element of area dS and unit normal n, the amount of chemical flowing through the surface element in an infinitesimal time interval, of duration dt, is given by n • q dSdt. (5.2) Definition. Pick's Law of Diffusion relates the fiux q to the gradient of c via q = -DVc, (5.3) where D, the diffusion coefficient, is independent of c and Vc. Bringing this together, we have, for any closed volume V (fixed in time and space), with bounding surface dV, dt Hence ^[cidV = -[ q-ndS+ [ Ri{ci,C2, ■ ■ ■ ,0^) dV, iG {!,..., m}. (5.4) Jv JdV JV ^^Qdy = J^V ■qdV + J^Riici,C2,...,Cm)dV (5.5) / {V ■{DVci)+Ri{ci,C2,...,Cm)} dV, (5.6) Jv Chapter 5. Introduction to spatial variation 45 and thus for any closed volume, V, with surface dV, one has ^1^- V-pVQ)-i?i|dy = 0, {!,..., m}. (5.7) Hence (JC' = V ■{DVci) + Ri, xeV, (5.8) ot which constitutes a system of reaction-diffusion equations for the m chemical species in the finite domain D. Such equations must be supplemented with initial and boundary conditions for each of the m chemicals. Warning. Given, for example, that r-2TT / cos6'd6l = 7^ cos 61 = 0, 6'G[0,27r], (5.9) Jo are you sure one can deduce equation (5.8)? Suppose (JC ' -^-V.{DVc,)-R,^0, (5.10) at some x = x*. Without loss of generality, we can assume the above expression is positive i.e. the left-hand side of equation (5.10) is positive. Then 3 e > such that _2 _ V . (DVq) - i?, > 0, (5.11) for aU X G Be{x*). In this case Br ^ - V • (DVq) - R^ 'B,{X*) contradicting our original assumption, equation (5.7). dV > 0, (5.12) Hence our initial supposition is false and equation (5.8) holds for x £ T>. Remark. With one species, with a constant diffusion coefficient, in the absence of reac- tions, we have the diffusion equation which in one dimension reduces to dc „ d^c For a given length scale, L, and diffusion coefficient, D, the timescale of the system is T = L'^/D. For a ceh, L ~ 10~^m = 10"^ cm, and for a typical protein D ~ 10 cm^s would not be unreasonable. Thus the timescale for diffusion to homogenise spatial gradients of this protein within a cell is L2 lO"*' cm2 r~ — 5 ^~10s, (5.14 D 10-7 cm2 s~i ^ ' therefore we can often neglect diffusion in a cell. However, as the scale doubles the time scale squares e.g. L x 10 ^ T x 100 and L x 100 ^ T x 10^. Chapter 5. Introduction to spatial variation 46 Note. The above derivation generalises to situations more general than modelling chem- ical or biochemical diffusion. For example, let I{x, y, t) denote the number of infected people per unit area. Assume the infectives, on average, spread out via a random walk mechanism and interact with susceptibles, as described in Section (6.2.1). One has that the flux of infectives, Qj, is given by qj = -DiVI, (5.15) where Dj is a constant, with dimensions of (length)^ (time)"^. Thus, one has, via precisely the same ideas and arguments as above, that dl — = \/.(DjVI) + rIS-aI, (5.16) where S{x, y, t) is the number of susceptibles per unit area, and r and a have the same interpretation as in Section 6.2.1. Fisher's Equation. A very common example is the combination of logistic growth and diffusion which, in one spatial dimension, gives rise to Fisher's Equation: Dirj, + ru(l--], (5.17) dt dx^ V K which was first proposed to model the spread of an advantageous gene through a popula- tion. See Section 6.1 for more details. 5.2 Chemotaxis As briefly mentioned earlier, motile cells can move in response to gradients in chemical concentrations, a process known as chemotaxis. This leads to slightly more complicated transport equations, as we shall see. The diffusive flux for the population density of the cells, n, is as previously: J d = —DVn. The flux due to chemotaxis (assuming it is an attractant rather than a repellent) is taken to be of the form: Jc = nx{a)Va = nV^{a), (5.18) where a is the chemical concentration and <&(a) increases monotonically with a. Clearly x(o) = *^'(a); the cells move in response to a gradient of the chemical in the direction in which the function ^{a) is increasing at the fastest rate. Thus the total flux is Jr) + Jc = -DVn + nx{a)Va. (5.19) Combining the transport of the motile cells, together with a term describing their repro- duction and/or death, plus an equation for the chemical which also diffuses and, typically, Chapter 5. Introduction to spatial variation 47 is secreted and degrades leads to the following equations — = V-{DVn)-V-{nx{a)Va) + f{7i,a), (5.20) da — = V • (DaVa) + An - fxa. (5.21) In the above the above /(n, a) is often taken to be a logistic growth term while the function x(o) describing the chemotaxis has many forms, including X(a) = ^2, (5,22) X(a) = (5.23) where the latter represents a receptor law, with $(a) taking a Michaelis-Menten form [5]. Chapter 6 Travelling waves Certain types of models can be seen to display wave-type behaviour. Here we will be interested in travelling waves, those that travel without change in shape and at constant speed. References. • J. D. Murray, Mathematical Biology Volume I, Chapter 13 [8]. • J. D. Murray, Mathematical Biology Volume II, Chapter 1 [9]. • N. F. Britton, Essential Mathematical Biology, Chapter 5 [1]. 6.1 Fisher's equation: an investigation Fisher's equation, after suitable non-dimensionalisation, is ^ = ^ + (6.1) where (3, z, t are all non-dimensionalised variables. Clearly the solution of these equations will depend on the initial and boundary conditions we impose. We state these conditions for the time being as /3(2;, t) — )• /3±oo as z — )• ±00 and f3{z,T = 0) = /3q{z), (6.2) where /3±oo, Po, are constants. 6.1.1 Key points • We will investigate whether such a wave solution exists for the above equations which propagates without a change in shape and at a constant (but as yet unknown) speed V. Such wave solutions are defined to be travelling wave solutions. 48 Chapter 6. Travelling waves 49 The investigation of the potential existence of a travelhng wave solution will be substantially easier to investigate on performing the transformation to the moving coordinate frame y = z — vt as, by the definition of a travelling wave, the wave profile will be independent of time in a frame moving at speed v. Using the chain rule and noting that we seek a solution that is time independent with respect to the y variable, we have d£_d£dy d£dT dp_ _ d^dy_ d£dT_ dt dy dt dr dt dz dy dz dr dz' i.e. dt dy dr dz dy' Assuming /3 = I3{y) so that d/S/dr = the partial differential equation, (6.1), reduces to /3" + v(3' + /3(1 - /?) = where ' = ^. (6.5) dy One must choose appropriate boundary conditions at ±00 for the travelling wave equations. These are the same as the boundary conditions for the full partial differ- ential equation (but rewritten in terms of y), i.e. P{y) P±oo as y ^ ±00, (6.6) where /3±oo) are the same constants as specified in equation (6.2). One must have that /3+cxd, /3-oo only take the values zero or unity: [/3" + vP' + /3(1 - /?)] dy = 0, (6.7) J —00 gives /oo /3(l-/3)dy = 0. (6.8) -00 If we want /3 — )■ constant as y — )■ ±00 and /3, /3' finite for Vy we must have either /3— T-Oor/?— T-lasy— )-oo and similarly for y — )■ —00. • With the boundary conditions (/?(— 00), /3(oo)) = (1,0), we physically anticipate u > 0. Indeed there are no solutions of the Fisher travelling wave equations for these bound- ary conditions and f < 0. Chapter 6. Travelling waves 50 • Solutions to equations (6.1) and (6.2) are unique. The proof would be an exercise in the theory of partial differential equations. • The solutions of the travelling wave equations are not unique. One may have solu- tions for different values of the unknown v. Also, if f3(y) solves (6.5) for any fixed value of V then, for the same value of v, so does /?(?/ + A), where A is any constant. For both v and A fixed the solution of the travelling wave equations are normally unique. • Note that the solutions of the travelling wave equations, (6.5), can only possibly be solutions of the full partial differential equation, when considered on an infi- nite domain. [Realistically one requires that the length scale of variation of the system in question is much less than the length scale of the physical domain for a travelling wave to (have the potential to) be an excellent approximation to the reaction-diffusion wave solutions on the physical, i.e. finite, domain]. • One "loses" the partial differential equation initial conditions associated with (6.1) and (6.2). The solution of the travelling wave equations given above for P are only a solution of the full partial differential equation, (6.1), for all time if the travelling wave solution is consistent with the initial conditions specified in (6.2). • However, often (or rather usually!), one finds that for a particular choice of v the solutions of the full partial differential equation system, (6.1) and (6.2), tend, as t — )• CO, to a solution of the travelling wave equations (6.5), with fixed v and A, for a very large class of initial conditions. • The Russian mathematician Kolmogorov proved that solutions of the full partial differential equation system, (6.1) and (6.2), do indeed tend, as t — )■ oo, to a solution of the travelling wave equations for v = 2 for a large class of initial conditions. 6.1.2 Existence and the phase plane We will investigate the existence of solutions of Fisher's equation, equation (6.5), with the boundary conditions (/3(— oo), /3(oo)) = (1, 0) and u > 0, by means of an extended exercise involving the phase plane (/?', /3). Consider the travelling wave equation d^/3 d/3 ^+.^ + /3(l-/3) = 0, (6.9) with V > and the boundary conditions (/?(— oo), /3(oo)) = (1,0). Exercise 1. Show that the stationary point at (/?,/?') = (1,0) is always a saddle point and the stationary point at (/3, /?') = (0, 0) is a stable node for v > 2 and a stable focus for V < 2. Chapter 6. Travelling waves 51 Solution. Writing /3' = 7 gives dy I 7 / dy\ (3' [ -vj- /3(1 - /3) (6.10) The Jacobian, J, is given by df dl f dj dg dg_ dp d-f At (0, 0) we have det(J - AI) = det I ^ , | \^ + vX + l = 0, (6.12) \ -1 -v-X and hence A = 1 . (6.13) Therefore: • if f < 2 we have A = —v/2 it ifi and hence a stable spiral; • if u > 2 we have A = —v/2 it /_f and hence a stable node; • if f = 2 we have A = — 1 and hence a stable node. At (1,0) we have det( J - AI) = det ( ^ \ A^ + uA - 1 = 0, (6.14) and hence A = (6.15) Therefore (1,0) is a saddle point. Exercise 2. Explain why solutions of Fisher's travelling wave equations must tend to phase plane stationary points as y — )• ±00. Hence explain why solutions of (6.9) with V < 2 are unphysical. Solution. {f3,-f) = (/3, f3') will change as y increases, unless at a stationary point. Therefore they will keep moving along a phase space trajectory as y — )• 00 unless the y — )■ 00 limit evolves to a stationary point. To satisfy lim^^o /3(y) = 0, we need to be on a phase space trajectory which "stops" at /3 = 0. Therefore we must be on a phase space trajectory which tends to a stationary point with /3 = as y — )• 00. Hence we must tend to (0, 0) as y — )■ 00 to satisfy liniy^oo /3(y) = as y — )• 00. An analogous argument holds as y — )• — 00. If v < 2 then /9 < at some point on the trajectory which is unphysical: Chapter 6. Travelling waves 52 7 13 Exercise 3. Show that the gradient of the unstable manifold at {f3,l3') = (1,0) is given by ^ ^ ^ (6.16) Sketch the qualitative form of the phase plane trajectories near to the stationary points for V > 2. Solution. We require the eigenvectors of the Jacobian at (1,0): ' 1 \ 1 1 -V 1 Q± Ah q± q± = X± and 1 — vq-i- = X±q±. (6.17) Hence v± 1 (6.18) Exercise 4. Explain why any physically relevant phase plane trajectory must leave {/3,/3') = (1,0) on the unstable manifold pointing in the direction of decreasing /3. Solution. Recall that, close to the stationary point, (3 7 /3* 7* a-e^-yv_+a+e^+yv+. (6.19) The solution moves away from the saddle along the unstable manifold, which corresponds to a_. Chapter 6. Travelling waves 53 def Exercise 5. Consider v > 2. With 7 = /5' show that for /3 € (0, 1] the phase plane trajectories, with gradient d/3 7 ^ satisfy the constraint dj/dfS < —1 whenever 7 = — /?. Solution. .20) d7 d^ -v + il-/3)< {-V + 2) - (1 + /3) < -1. .21) /9=-7 Exercise 6. Hence show that with v >2 the unstable manifold leaving (/3,/3') = (1,0) and entering the region /?' < 0, /3 < 1 enters, and can never leave, the region n=^{{P,j) I 7<0, /3g [0,1], 7>/3}. (6.22) Solution. Along £1 = {(/3,7) | 7 = 0, /3 G (0, 1)} the trajectories point vertically into TZ as 00 as we approach £1 and 7' = -/3(1 - /3) < 0. (6.23) d7 d/3 Along £2 = {(/3,7) I ,5 = 1, 7 E (-1,0)} we have d7 dp C2 /3(l-/3) 7 -V < 0. (6.24) Hence trajectories that enter TZ cannot leave. There any trajectory must end at a station- ary point and trajectories are forced to the point (/3,7) = (0,0). Exercise 7. Thus prove that that there exists a monotonic solution, with /? > 0, to equation (6.9) for every value of > 2 and, with v > 2 fixed, the phase space trajectory is unique. Solution. The above analysis is valid for v > 2. For v fixed a trajectory enters the region TZ along the unstable manifold (only one unstable manifold enters TZ). The solution is monotonic as 7 < throughout TZ. Figure 6.1 shows the results of numerical simulation of the Fisher equation (6.1) with initial and boundary conditions given by (6.2) at a series of time points. 6.1.3 Relation between the travelling wave speed and initial conditions We have seen that, for v fixed, the phase space trajectory of Fisher's travelling wave equation is unique. The non-uniqueness associated with the fact that if /3(y) solves Fisher's travelling wave equation then so does /3{y + A) for A constant simply corresponds to a shift along the phase space trajectory. This, in turn, corresponds simply to translation of the travelling wave. Chapter 6. Travelling waves 54 1.0 0.8 0.6 0.4 0.2 0.0 50 100 150 z Figure 6.1: Solution of the Fisher equation (6.f) witli initial and boundary conditions given by (6.2) at times t = 10, 20, 30, 40, 50. Key question. Do solutions of the full system, equations (6.1) and (6.2), actually evolve to a travelling wave solution, and if so, what is its speed? Non-Examinable: initial conditions of compact support Kolmogorov considered the equation = g + (6.25) with the boundary conditions ■0(z,t)— 7-1 as z — 7- — CX3 and ^/'(z,r)— t-O as z — t- oo, (6.26) and non-negative initial conditions satisfying the following: there is a K, with < K < oo, such that il;(z,T = 0) = for z>K and ilj{z,T = 0) = 1 for z < -K. (6.27) He proved that ^l^{z, r) tends to a Fisher travelling wave solution with f = 2 as t — t- oo. This can be applied to equations (6.1) and (6.2) providing the initial conditions are non- negative and the initial condition for /3 satisfies the above constraint, i.e. there is a K, with < i^T < oo, such that P{z, r = 0) = for z> K and f3{z, r = 0) = 1 for z < -K. (6.28) Under such constraints (3 also tends to a Fisher travelling wave solution with v = 2. 6.2 Models of epidemics The study of infectious diseases has a long history and there are numerous detailed models of a variety of epidemics and epizootics (i.e. animal epidemics). We can only possibly scratch the surface. In the following, we consider a simple, framework model but even this is capable of highlighting general comments about epidemics and, in fact, approximately describes some specific epidemics. Chapter 6. Travelling waves 55 6.2.1 The SIR model Consider a disease for which the population can be placed into three compartments: • the susceptible compartment, S, who can catch the disease; • the infective compartment, I, who have and transmit the disease; • the removed compartment, R, who have been isolated, or who have recovered and are immune to the disease, or have died due to the disease during the course of the Assumptions • The epidemic is of short duration course so that the population is constant (counting those who have died due to the disease during the course of the epidemic) . • The disease has a negligible incubation period. • If a person contracts the disease and recovers, they are immune (and hence remain in the removed compartment). • The numbers involved are sufficiently large to justify a continuum approximation. • The 'dynamics' of the disease can be described by applying the Law of Mass Action The model Then the equations describing the time evolution of numbers in the susceptible, infective and removed compartments are given by epidemic. to: r a (6.29) d5 d^ d/ di dR rIS — al al. rIS, (6.30) (6.32) subject to S{t = 0) = 5o, I{t 0) = lo, R{t = 0) = 0. Note that Key questions in an epidemic situation are, given r, a, 5*0 and Iq Chapter 6. Travelling waves 56 1. Will the disease spread, i.e. will the number of infectives increase, at least in the short-term? Solution. — = —rIS =^ is decreasing and therefore S < Sq. (6.35) ^ = I{rS -a) < I{rSo - a). (6.36) Therefore, if Sq < a/r the infectives never increase, at least initially. 2. If the disease spreads, what will be the maximum number of infectives at any given time? Solution. d/ {rS-a) p def a - = ^ = -1 + ^ where p = -. (6.37) Integrating gives I + S — plnS = Iq + So — plnSo, (6.38) and so, noting that dl/dS = for S = p, the maximum number of infectives is given by Ima. = I ^° ^° - ^ . (6.39) \^ Io + So- plnSo- plnp- p Sq > p 3. How many people in total catch the disease? Solution. From 2, / — t- as t — t- oo. Therefore the total number who catch the disease is R{oo) = No- S{oo) - /(oo) = No- 5(oo), (6.40) where S'(oo) < So is the root of ^oo - yolnS'oo = iVo - yolnS'o, (6-41) obtained by setting S = Soo and No = lo + So in equation (6.38). 6.2.2 An SIR model with spatial heterogeneity We consider an application to fox rabies. We will make the same assumptions as for the standard SIR model, plus: • healthy, i.e. susceptible, foxes are territorial and, on average, do not move from their territories; • rabid, i.e. infective, foxes undergo behavioural changes and migrate randomly, with an effective, constant, diffusion coefficient D; • rabies is fatal, so that infected foxes do not return to the susceptible compartment but die, and hence the removed compartment does not migrate. Chapter 6. Travelling waves 57 t X. X 20 40 60 80 100 S Figure 6.2: Numerical solution of the SIR model, equations (6.30)-(6.32), where the solid lines indicate the phase trajectories and the dashed line S' + / = 5o + /o- Parameters are as follows: r = 0.01 and a = 0.25. Taking into account rabid foxes' random motion, the SIR equations become — = -rIS, (6.42) _ = DV^I + rIS -al, (6.43) f = («-^^> The / and 5 equations decouple, and we consider these in more detail. We assume a one-dimensional spatial domain x E (—00,00) and apply the following scalings/non- dimensionalisations , h = ^, 5"* = -^, x^ = J^x, U=rSot, A = -^, (6.45) 60 60 V rSo rSo where Sq is the population density in the absence of rabies, to obtain ^ = V^I + I{S-X), (6.47) where asterisks have been dropped for convenience in the final expression. Travelling waves We seek travelling wave solutions with S{x,t) = S{y), I{x,t) = I{y), y = x - ct, c> 0, (6.48) which results in the system = cS'-IS, (6.49) = I" + cl' + I{S - X), (6.50) Chapter 6. Travelling waves 58 where ' = d/dy. We assume A = a/(r5o) < 1 below. This is equivalent to the condition for disease spread in the earlier SIR model. Boundary conditions We assume a healthy population as ?/ — )■ oo: 5^1 and / ^ 0, (6.51) and as y — )■ — oo we require 1^0. (6.52) Bound on travelling wave speed We write S = 1 — P and linearise about the wavefront: -cP'-I = and I" + cl' + X). (6.53) The / equation decouples and analysis of this equation gives a stable focus at (/, /') = (0, 0) if the eigenvalues ,= Z^±VIIWER, (0.54) are complex. This requires c > 2Vl - A. (6.55) Severity of epidemic S{oo) is a measure of the severity of the epidemic. We have / = cS' / S and therefore ^(/' + cl) + cS' (^^^) = 0. (6.56) Therefore (/' + cl) + c(5 - A In S) = constant = c, (6.57) by evaluating the equation as — )• oo. In this case S'(-oo) - AlnS'(-oo) = 1, where S'(-oo) < 1, (6.58) gives the severity of the epidemic. Further comments on travelling wave speed Typically, the wave evolves to have minimum wave speed: C ^ Cmin- (6.59) Chapter 7 Pattern formation Examples of the importance of spatial pattern and structure can be seen just about every- where in the natural world. Here we will be concerned with building and analysing models which can generate patterns; understanding how self-organising principles may lead to the generation of shape and form. References. • J. D. Murray, Mathematical Biology Volume II, Chapter 2 and Chapter 3 [9]. • N. F. Britton, Essential Mathematical Biology, Chapter 7 [1]. 7.1 Minimum domains for spatial structure Consider the non-dimensionalised, one dimensional, budworm model, but with a diffusive spatial structure: We also suppose that exterior to our domain conditions are very hostile to budworm so that we have the boundary conditions where L > is the size of the domain. Note that /'(O) = r > 0. Question. Clearly u = is a solution. However, if we start with a small initial distri- bution of budworm, will we end up with no budworm, or an outbreak of budworm? In particular, how does this depend on the domain size? Solution. For initial conditions with < u{x,t = 0) ^ 1, sufficiently small, we can approximate f{u) by f'{0)u at least while u{x,t) remains small. Thus our equations are, approximately. (7.3) 59 Chapter 7. Pattern formation 60 We look for a solution of the form (invoking completeness of Fourier series) : oo ■r-^ . /nvra u[x, t) = an{t) sm 1^— (7.4) n=l This gives that the time-dependent coefficients satisfy dan Dn^vr^ dt L2 -«n + /'(0)an = O-nfln, and hence ^{x, t) = ^ «n ^ exp n=l sm nvrx (7.5) (7.6) For the solution to decay to zero, we require that all Fourier modes decay to zero as t — )■ oo, and hence we require that a„ < Vn /'(0)-^^^<0 Vn, and thus that /'(O) < L < L2 Dtt^ 7W) def L crit- (7.7) (7.8) Hence there is a critical lengthscale, Lcrit-, beyond which an outburst of budworm is possible in a spatially distributed system. 7.1.1 Domain size On first inspection one probably should be surprised to see that Lent increases linearly with the diffusion coefficient, i.e. diffusion is destabilising the zero steady state. We can further investigate how the nature of a steady state pattern depends on the diffusion coefficient. Suppose L > Lcrit and that the steady state pattern is of the form: u , i ^max -L/2 X - = L/2 X We therefore have = Du^^r. + f{u)- Multiplying by Ux and integrating with respect to x, we have = y DuxUxx dx + J Uxf{u) dx. (7.9) (7.10) Chapter 7. Pattern formation 61 Thus we have -Dul + F{u) constant = F{umax) where F'{u) = f{u) (7.11) We can therefore find a relation between L, D, integrals of F{u)=^ r f{y)dy, and max(u), the size of the outbreak, as follows: 2 \ 2 D (7.12) \/ F{umax) — F{u) since x > and therefore Ux < 0. (7-13) Integrating, gives rL/2 ^ n 2 dx = -(2L»)2 / Jo Ju L = {2D) and hence 1 -tax '\J F (Ufyiax) 1 \/F{Umax) - F{u) F{u) du. : dtt, Therefore Umax is a function of L/^2D and the root of equation (7.15). (7.14) (7.15) Figure 7.1: Numerical simulation of the u,n-L space, equation (7.15) with r = 0.6, q ~ 6.2 and D = 0.1. 7.2 Diffusion- driven instability Consider a two component system ut = DuV^u + f{u,v), (7.16) vt = D^V^v + g{u,v), (7.17) Chapter 7. Pattern formation 62 for a; e 0, t G [0, oo) and O bounded. The initial conditions are u{x,0) = uq{x), v{x,0) = vo{x), (7.18) and the boundary conditions are either Diriclilet, i.e. u = ub, V = vb, X G dil, (7.19) or homogeneous Neumann, i.e. n • Vn = 0, n • = 0, for x e dU, (7.20) where n is the outward pointing normal on dU. Definition. Patterns are stable, time-independent, spatially heterogeneous solutions of equations (7.16)-(7.17). Definition. A diffusion- driven instability, also referred to as a Turing instability, occurs when a steady state, stable in the absence of diffusion, goes unstable when diffusion is present. Remark. Diffusion-driven instabilities, in particular, can drive pattern formation in chemical systems and there is significant, but not necessarily conclusive, evidence that it can drive pattern formation in a variety of biological systems. A key point is that this mechanism can drive the system from close to a homogeneous steady state to a state with spatial pattern and structure. The fact that diffusion is responsible for this is initially quite surprisingly. Diffusion, in isolation, disperses a pattern; yet diffusion, when in combination with the kinetic terms, often can drive a system towards a state with spatial structure. 7.2.1 Linear analysis We wish to understand when a diffusion-driven instability occurs. Using vector and matrix notation we define uJ-]. F(u)J '^"•"l), oJ'i' » V (7.21) V ^ 9{u,v) y ' y ^ and write the problem with homogeneous Neumann boundary conditions as follows: ut = DV'^u + F{u), (7.22) l(:)^("o"«>^(:)^(:i"::;). <-> with n • Vu = 0, xedn, (7.24) Chapter 7. Pattern formation 63 I.e. n.Vu = = n.Vv x G dn. (7.25) Let u* be such that F{u*) = 0. Implicit in this definition is the assumption that u* is a constant vector. Let w = u — u* with \w\ <^ 1. Then we have duo — — = DV^w + F{u*) + Jw + higher order terms, (7.26) ot where / df df \ du dg dv dg \ du dv J (7.27) u=u* is the Jacobian of F evaluated at u = u* . Note that J is a constant matrix. Neglecting higher order terms in \w\, we have the equation wt = DV^w + Jw, n ■ Vw = 0, x e 80.. (7.28) This is a linear equation and so we look for a solution in the form of a linear sum of separable solutions. To do this, we first need to consider a general separable solution given by w{x,t) = A{t)p{x), (7.29) where A{t) is a scalar function of time. Substituting this into equation (7.28) yields 1 dA - — p = DV^p + Jp. (7.30) Clearly to proceed, with p dependent on x only, we require A/A to be time independent. It must also be independent of a; as ^ is a function of time only. Thus A/A is constant. We take A = XA, where A is as yet an undetermined constant. Thus A = Aoexp{\t), (7.31) for Aq ^ constant. Hence we require that our separable solution is such that [Xp - Jp- UVV] = 0. (7.32) Suppose p satisfies the equation V^p + k^p = 0, n • Vp = 0, a; G (7.33) where A; G M. This is motivated by the fact in one-dimensional on a bounded domain, we have p" + k'^p = 0; the solutions are trigonometric functions which means one immediately has a Fourier series when writing the sum of separable solutions. Chapter 7. Pattern formation 64 Then we have [Xp -Jp + Dk^p] = 0, (7.34) and thus [\I - J + Dk'^]p = Q, (7.35) with IpI not identically zero. Hence det [\I-J + k'^D] = 0. (7.36) This can be rewritten as y -Qu X- gv + Di,k^ J which gives the following quadratic in A: + [{Du + D,)k^ - {fu + g,)] A + h{e) = 0, (7.38) where hik") = DuD^k^ - {D„fu + D^g,)^ + {fuQv - 9ufv)- (7.39) Note 1. Fixing model parameters and functions (i.e. fixing Z)„, D^, /, (7), we have an equation which gives A as a function of k"^. Note 2. Thus, for any A:^ such that equation (7.33) possesses a solution, denoted Pk{x) below, we can find a A = A(/i;^) and hence a general separable solution of the form Aoe^^^'^'Pkix). (7.40) The most general solution formed by the sum of separable solutions is therefore J]Ao(fc2)e^('=')*Pfc(^), (7.41) if there are countable k"^ for which equation (7.33) possesses a solution. Otherwise the general solution formed by the sum of separable solutions is of the form j Ao{k^)e^^''"^'p,.{x)de, (7.42) where k'^ is the integration variable. Unstable points If, for any /c^ such that equation (7.33) possesses a solution, we find Re{X{k^)) > then: • u* is (linearly) unstable and perturbations from the stationary state will grow; • while the perturbations are small, the linear analysis remains valid; thus the per- turbations keep growing until the linear analysis is invalid and the full non-linear dynamics comes into play; Chapter 7. Pattern formation 65 • a small perturbation from the steady state develops into a growing spatially hetero- geneous solution which subsequently seeds spatially heterogeneous behaviour of the full non-linear model; • a spatially heterogeneous pattern can emerge from the system from a starting point which is homogeneous to a very good approximation. Stable points If, for all fc^ such that equation (7.33) possesses a solution, we find Re{\{k'^)) < then: • u* is (linearly) stable and perturbations from the stationary state do not grow; • patterning will not emerge from perturbing the homogeneous steady state solution u*- • the solution will decay back to the homogeneous solution^. 7.3 Detailed study of the conditions for a Turing instability For a Turing instability we require the homogeneous steady state to be stable without diffusion and unstable with diffusion present. Here we analyse the requirements for each of these conditions to be satisfied. 7.3.1 Stability without diffusion We firstly require that in the absence of diffusion the system is stable. This is equivalent to i?e(A(0)) < 0, (7.43) for all solutions of A(0), as setting /c^ = removes the diffusion-driven term in equation (7.36) and the preceding equations. We have that A(0) satisfies A(0)2 - [/„ + g,] A(0) + [UQv - fvQu] = 0. (7.44) Insisting that i?e(A(0) < 0) gives us the conditions fu + gv < (7.45) fu9v-fv9u > 0. (7.46) '^Technical point: Strictly, this conclusion requires completeness of the separable solutions. This can be readily shown in ID on bounded domains. (Solutions of p" + k^p = on bounded domains with Neumann conditions are trigonometric functions and completeness is inherited from the completeness of Fourier series). Even if completeness of the separable solutions is not clear, numerical simulations of the full equations are highly indicative and do not, for the models typically encountered, contradict the linear analysis results. With enough effort and neglecting any biological constraints on model parameters and functions, one may well be able to find Du, D^, f, g where there was such a discrepancy but that is not the point of biological modelling. Chapter 7. Pattern formation 66 The simplest way of deducing (7.45) and (7.46) is by brute force. The roots of the quadratic are given by w„^ ifu + 9v) ± V {fu + QvY - ^{fu9v - fvOu) , A(Uj± = . (7.47) 7.3.2 Instability with diffusion Now consider the effects of diffusion. In addition to i?e(A(0)) < 0, we are required to show, for diffusion-driven instabihty, that there exists k"^ such that Re{\{k^)) > 0, (7.48) so that diffusion does indeed drive an instabihty. We have that \{k'^) satisfies a2 + [{Du + D,)k' - ifu + g,)] X + hiP) = 0, (7.49) where h{k^) = DuD.k" - {D,fu + D^g,)k^ + (/„5„ - gufv), (7.50) and a = {fu + gv)-{Du + D,)k^ <0. (7.51) Thus Re{X{k^)) > requires that Re (^a ± ^a^ - 4h{k^)j > ^ h{k'^) < 0. (7.52) Hence we must find k'^ such that h{k^) = DuD.k'' - {D,fu + D^g,)k^ + if^g, - g^h) < 0, (7.53) so that we have /c^ G [A;^,A;^] where hi^kj.) = 0. Figure 7.2 shows a plot of a caricature /i(A;2). This gives us that we have an instability whenever 'a - _ B A + _ B k^G [kl,kl], (7.54) where A = Dyfu + Dugv and B = ADuD^{fugv - gufv) > 0, (7.55) and there exists a solution of the following + k"^? = 0, n • Vp = 0, x£dn, (7.56) for A;2 in the above range. Chapter 7. Pattern formation 67 Figure 7.2: A plot of a caricature h{k'^). Insisting that k is real and non-zero (we have considered the k = case above) we have A>0 and A"^ - B > 0, (7.57) which gives us that when Re{X{k'^)) > 0, the following conditions hold: ^ > : D^fu + DuQv > 0, (7.58) A^-B>0: {DJu + A.5t.) > 2^DMfu9v - fv9u)- (7.59) 7.3.3 Summary We have found that a diffusion-driven instability can occur when conditions (7.45), (7.46), (7.58), (7.59) hold whereupon the separable solutions, with k'^ within the range (7.54) and for which there is a solution to equation (7.33), will drive the instability. Key point 1. Note that constraints (7.45) and (7.58) immediately gives us that ^ D2. Thus one cannot have a diffusion-driven instability with identical diffusion coefficients. Key point 2. Prom constraints (7.45), (7.46), (7.58) the signs of /„, must be such that J takes the form (7.60) Key point 3. A Turing instability typically occurs via long-range inhibition, short-range activation. In more detail, suppose (7.61) Chapter 7. Pattern formation 68 Then we have fu>0 and ^r^, < by the signs of J. In this case Dyf^ + D^g^ > D2 > Du- Hence the activator has a lower diffusion coefficient and spreads less quickly than the inhibitor. 7.3.4 The threshold of a Turing instabiHty. The threshold is defined such that equation (7.39), i.e. D^Dykt - {D,U + Dugv)kl + UnQv - 9ufv) = 0, (7.62) has a single root, k^. Thus we additionally require A^ = B i.e. {Dyfu + Dugv? = ADMfu9v-gufv)>Q, (7.63) whereupon ,2 ^ ^ ^ ^-"f^ + ^^9v ^ .X Strictly one also requires that a solution exists for V^p + /c^p = 0, n • Vp = 0, X e dn, (7.65) when k"^ = k^. However, the above value of k1 is typically an excellent approximation. 7.4 Extended example 1 Consider the one-dimensional case: ut = DuUxx + f{u,v), (7.66) vt = DyVxx + g{u,v), (7.67) for a; € [0, L], t G [0, 00) and zero flux boundary conditions at x = and x = L. The analogue of V^p + k'^p = 0, n • Vp = 0, xe dn, (7.68) is Pxx + k^P = 0, p'{0) = p'{L) = 0, (7.69) which gives us that TITT Pk{x) = Ak cos (kx) , k = —, n G {1,2, . . .}, (7.70) where A^ is fc-dependent in general but independent of t and x. Thus the separable solution is of the form cos (kx) , (7.71) k where the sum is over the allowed values of k i.e. fe = ^, nE{l,2,...}. (7.72) Chapter 7. Pattern formation 69 Activator Inhibitor 400 400 I I 20 40 60 80 distance (x) Figure 7.3: Numerical simulation of the Giercr-Meinhardt model for pattern formation. 7.4.1 The influence of domain size If the smallest allowed value of k"^ = vr^/L^ is such that vr = — > then we cannot have a Turing instability. (7.73) Thus for very small domains there is no pattern formation via a Turing mechanism. However, if one slowly increases the size of the domain, then L increases and the above constraint eventually breaks down and the homogeneous steady state destabilises leading to spatial heterogeneity. This pattern formation mechanism has been observed in chemical systems. It is regularly hypothesised to be present in biological systems {e.g. animal coat markings, fish markings, the interaction of gene products at a cellular level, the formation of ecological patchiness) though the evidence is not conclusive at the moment. 7.5 Extended example 2 Consider the two-dimensional case with spatial coordinates x = (x,y)^, x G [0, a], y € [0, 5], and zero flux boundary conditions. We find that the allowed values of /c^ are k^ 62 with /mvrxx /"mry Pm,n{^) = Am,n COS \ COS (^-^ excluding the case where n, m are both zero. n,m £ {0, 1,2,...}, (7.74) (7.75) Chapter 7. Pattern formation 70 Suppose the domain is long and thin, b <ti a. We may have a Turing instabihty if — ^ + 62 G [k'i,kl] where = 0. (7.76) For b sufficiently small, this requires n = and therefore no spatial variation in the y direction. This means we have that the seed for pattern formation predicted by the linear analysis is a separable solution which is "stripes" ; this typically invokes a striped pattern once the non-linear dynamics sets in. For a large rectangular domain, 5 ~ a sufficiently large, it is clear that a Turing instability can be initiated with n, m > 0. This means we have that the seed for pattern formation predicted by the linear analysis is a separable solution which is "spots". This typically invokes a spotted pattern once the non-linear dynamics sets in. 10 10 Figure 7.4: Changes in patterning as the domain shape changes. Figure 7.4 shows how domain size may affect the patterns formed. On the left-hand side the domain is long and thin and only a striped pattern results, whilst the on the right-hand side the domain is large enough to admit patterning in both directions. Suppose we have a domain which changes its aspect ratio from rectangular to long and thin. Then we have the following possibilities: Chapter 7. Pattern formation 71 This leads to an interesting prediction, in the context of animal coat markings, that if it is indeed driven by a Turing instability, then one should not expect to see an animal with a striped body and a spotted tail. Figure 7.5: Animal coat markings which are consistent with the predictions of pattern formation by a Turing instability. Common observation is consistent with such a prediction (see Figure 7.5) but one should not expect universal laws in the realms of biology as one does in physics (see Figure 7.6). More generally, this analysis has applications in modelling numerous chemical and biochemical reactions, in vibrating plate theory, and studies of patchiness in ecology and modelling gene interactions. Figure 7.6: Animal coat markings which are inconsistent with the predictions of pattern formation by a Turing instability. Chapter 8 Excitable systems: nerve pulses Many cells communicate with one another via nerve impulses, also known as action po- tentials. Action potentials are brief changes in the membrane potential of a cell produced by the flow of ionic current across the cell membrane. Such type of communication is not limited to neurons by make occur in other cells, for example cardiac and muscle cells. See http://en.wikipedia.org/wiki/Action_potential and related links for more de- tails. References. • J. P. Keener and J. Sneyd, Mathematical Physiology, Chapter 4 and Chapter 8 [7]. 8.1 Background Here we outline the background physics required to write down a model to describe a nerve impulse. Firstly, we note that: • numerous fundamental particles, ions and molecules have an electric charge, e.g. the electron, e~, and the sodium ion, Na~^; • it is an empirical fact that total charge is conserved; • electric charges exert electrical forces on one another such that like charges repel and unlike charges attract. The electric potential, denoted V, is the potential energy of a unit of charge due to such forces and is measured in volts; • a concentration of positive particles has a large positive potential, while a concen- tration of negative particles has a large, but negative potential; • electric current is defined to be the rate of flow of electric charge, measured in Amps. 72 Chapter 8. Excitable systems: nerve pulses 73 8.1.1 Resistance Ohms Law, ISV = IR, holds in most situations, where AV is the change in potential, / is the current flowing and R, which may depend on material properties and geometries but not on I nor V is the resistance. > V + AV R I V Key point. Suppose one uses a wire of low resistance to connect a region with a con- centration of positive charges to a region with a concentration of negative charges. The charges will, very quickly, flow onto/off the wire until the potential is constant and there is no further flow of charge. 8.1.2 Capacitance A simple example of capacitor is two conducting plates, separated by an insulator, for example, an air gap. + V Connecting a battery to the plates, as illustrated, using wires of low resistance leads to charge flowing onto/off the plates. It will equilibrate (very quickly!); let Qeqm denote the difference in the total value of the charge stored on the two plates. The capacitance of the plates, C, is defined to be C = ^>0, (8.1) where C is a constant, independent of V. Thus the higher the capacitance, the better the plates are at storing charge, for a given potential. Chapter 8. Excitable systems: nerve pulses 74 V{t) Suppose now the potential is a function of time, V = V{t). Charge will flow on and off the plates in response to time dependent changes in V{t). If the time for a significant change in V{t) is far longer than the time it takes for the difference in the total value of the charge stored on the two plates, Q, to reach its equilibrium value, Qegm, as is essentially always the case, one has Q = Qeqm. = CV{t). (8.2) Hence, the current, J, i.e. the rate of flow of charge on/off the plates is given by J = Q = CV. (8.3) 8.2 Deducing the Fitzhugh Nagumo equations An axon is a part of nerve cell: The nerve signal along the axon is, in essence, a propagating pulse in the potential dif- ference across the plasma (i.e. outer) membrane of the axon. This potential difference, V, arises due to the preferential permeability of the axon plasma membrane which allows potassium and sodium ions, K+ and iVa+, to pass through the membrane at rates which differ between the two ions and vary with V. In the rest state, V = Vrest — — 70mV (millivolts); in a nerve signal pulse in V rises to a peak of ~ 15mV. It is this pulse we are interested in modelling. The geometry of the axon can be treated as a cylindrical tube. An axon is axisymmetric, so we have no 9 dependence in any of our models of the axon. 8.2.1 Space-clamped axon A common, simplifying, experimental scenario is to space-clamp the axon, i.e. to place a conducting wire along the axon's axis of symmetry. Chapter 8. Excitable systems: nerve pulses 75 • Thus, by conservation of charge, the total current flowing across the axon membrane must be zero. • Note that any changes in the interior due to, for example the axon allowing and Na^ to pass through its membrane, will occur on a much slower timescale and hence one has that the interior of the space-clamped axon has no spatial variation in its potential difference, no current flowing along the inside of the axon, and, most importantly, the total current flowing through the axon membrane is zero. The basic model for the space- clamped axon plasma membrane potential is given by = total transmembrane current per unit area, dV = C— + lNa + Ik + Io + Iapplied{t), (8.4) where • Iappiied{t) is the applied current, i.e. the current injected through the axon plasma membrane in the experiment, which is only function of time in most experimental set-ups. We will take Iappiied{t) = below. Chapter 8. Excitable systems: nerve pulses 76 • cdV/dt is the capacitance current through a unit area of the membrane. Recall: Rate of flow on/off capacitor = C dV/dt. Therefore rate of flow of charge per unit area of membrane = c dV/dt where c is the membrane capacitance per unit area. • iNa, Ik are the voltage dependent Na'^ and currents. Iq is a voltage dependent background current. These currents actually take complicated forms involving numerous other variables which satisfy complex equations, that can be simplified, if somewhat crudely. An excellent account is given in T. F. Weiss, Cellular Biophysics [12]. The resulting equations written in terms of the non-dimensional variables v = {V — yrest)/\yrest\ and T = t/T where T = 6 ms, the time scale of a typical nerve pulse, are dv e— = Av{5 - v){v - 1) - n, (8.5) dr dn — = -^n + v, (8.6) dr where A, 7, e, 6 are positive parameters such that A, 7 ~ 0(1), < e <^ 5 <^ 1. Key point. The spatially independent behaviour of a space-clamped axon is approxi- mated by the above Fitzhugh Nagumo equations, (8.5)-(8.6). 8.3 A brief look at the Fitzhugh Nagumo equations We have = Av(5 - v)(v - 1) - n, (8.7) dr ^ = -7n + v, (8.8) where A, 7, e, 6 are positive parameters such that ^,7 ~ 0(1), < e <^ 5 <^ 1. 8.3.1 The {n,v) phase plane The nullclines of equations (8.5)-(8.6) are the lines where v = Q and n = 0. A plot of the nullclines separates the {v,n) phase plane into four regions, as shown in Figure 8.1. There are several things to note about the dynamics. • There is one stationary point which is a stable focus. • Thus, with initial conditions sufficiently close to the stationary point, the system evolves to the stationary point in a simple manner. Chapter 8. Excitable systems: nerve pulses 77 0.20 0.15 0.10 0.05 0.00 -0.5 0.0 0.5 1.0 1.5 V Figure 8.1: The phase plane for the Fitzhugh-Nagumo equations with the v nullcline shown in red and the n nullcline in green. The trajectories for two different initial perturbations from the steady state are shown as dashed lines. Parameters are as follows: A — 1, j = 0.5, S = 0.1 and e = 0.001. • Consider initial conditions with n ~ 0, but v increased sufficiently. The system does not simply relax back to the equilibrium. However, one can understand how the qualitative behavioural of the system by considering the phase plane. • We anticipate that v = (V — Vrest)/\yrest \ behaves in the manner shown in Figure 8.2 for a sufficiently large perturbation in v. 0.01 0.4 0.6 time (t) 0.00 -0.01 1.5 2.0 time (t) Figure 8.2: Solutions of the Fitzhugh Nagumo equations with v dynamics indicated by the solid line and n dynamics by the dashed line. The right-hand figure shows the oscillations that arise for large t. Parameters are as follows: A = 1, j = 0.5, S = Q.l and e = 0.01. • This is essentially a nerve pulse (although because of the space clamping all the nerve axon is firing at once). Definition. A system which, for a sufficiently large perturbation from a stationary point, undergoes a large change before eventually returning to the same stationary point is referred to as excitable. Chapter 8. Excitable systems: nerve pulses 78 8.4 Modelling the propagation of nerve signals In the following, we generalise the ideas we have seen for modelling the plasma membrane potential of an axon to scenarios where this potential can vary along the axon. 8.4.1 The cable model In the model we are about to develop we make following assumptions. • The cell membrane is a cylindrical membrane separating two conductors of electric current, namely the extracellular and intracellular mediums. These are assumed to be homogeneous and to obey Ohm's law. • The model has no 9 dependence. • A circuit theory description of current and voltages is adequate, i.e. quasi-static terms of Maxwell's equations are adequate; for example, electromagnetic radiation effects are totally negligible. • Currents flow through the membrane in the radial direction only. • Currents flow through the extracellular medium in the axial direction only and the potential in the extracellular medium is a function of z only. Similarly for the potential in the intracellular medium. These assumptions are appropriate for unmyelinated nerve axons. Deriving the model requires considering the following variables: • Ie{z,t) - external current; • Ii{z,t) - internal current; • J{z, t) - total current through the membrance per unit length; • Jioniz,t) - total ion current through the membrance per unit area; • V{z,t) = Vi{z,t) — Ve{z,t) - transmembrane potential; • ri - internal resistance per unit length; • re - external resistance per unit length; • C - membrane capacitance per unit area. Chapter 8. Excitable systems: nerve pulses 79 Veiz) Ie{z)dz Ie{z + dz)dz *^Ve{z + dz) cell exterior -AAAA/ rr,dz 2TraJ^on(z)dz J{z)dz • V^{Z) r.,dz 2-KaC^dz cell membrane h{z)dz V,{z + dz) J{z + dz)dz Ii{z + dz)dz cell interior Consider the axial current in the extracellular medium, which has resistance rg per unit length. We have dV Ve{z + dz)-Ve{z) = -rJ,{z)dz rj,{z) = -^, (8.9) where the minus sign appears because of the convention that positive current is a flow of positive charges in the direction of increasing z. Hence, if Ve{z + dz) > Ve{z) then positive charges flow in the direction of decreasing z giving a negative current. Similarly, TiUz) = (8.10) Using conservation of current, we have Ie{z + (lz,t) - Ie{z,t) = J{z,t)dz = Ii{z,t) - Ii{z + dz,t), (8.11) which gives J(r ^) = - _ dz dz J{z,t) = -^ = ^. (8.12) Hence and so (8 13) |^ = (r.+re)J. (8.14) Putting this all together gives ^_ d{Ii + Ie) _ d fl dVe ^ 1 dVi\ _ / re + rA d'^Ve ^ 1 d'^V ^^^^^ dz dz \re dz Vi dz J \ VeVi J dz"^ ri dz"^ ' and so = ~ ^ = ~ 'JTY + {'re+ri)J{z,t) Ti dz'^ \ Ti J dz Ti \ dz'^ J .16) Chapter 8. Excitable systems: nerve pulses 80 We also have that J{z, t) = 2na (^J^on{V, z, t) + ) , (8.17) and finally, therefore, 1 a^y av This gives an equation relating the cell plasma membrane potential, V ^ to the currents across the cell plasma membrane due to the flow of ions, Jion{V, z, t). Note 1. Note even though, physically, there is no diffusion, we still have a parabolic partial differential equation, so the techniques we have previously studied are readily applicable. Note 2. From the above equation one can model cell plasma membrane potentials given suitable initial and boundary conditions, and a suitable expression for Jjo„(z,t). We use the same expression for Jion{z,t), i.e. the expression for I^a + Ik + Iq as in the Fitzhugh Nagumo model of a space-clamped axon. Thus with V = {V — Vrest)/\yrest\ and X = Kz, where K is a constant, we have = e'^^+Av{6-v){v-l)-n, (8.19) dn _ = -^n + v, (8.20) dr where < ^,7 ~ 0(1), < e < 5 < 1. Note that K has been chosen so that the coefficient infront of the Vxx term is e^. This means, with respect to such variables, the front of the nerve pulse is extremely sharp. Hence, for such a scaling to exist, the extent of the nerve pulse must be less than eL, where L is the length of the axon; this constraint holds true for typical parameter estimates. The reason for the choice of this scaling is simply mathematical convenience in a travelling wave analysis. We are interested in nerve pulses, so we take the boundary conditions to be n, — )■ as X — )• zboo. We thus again have a system of parabolic partial differential equations to solve, and we are particularly interested in travelling pulse solutions. This entails that a travelling wave analysis would be most insightful. With the travelling wave coordinate y = x — ct and v{y) = v{x,t), n{y) = n{x,T), we obtain + ec^ + Av{5 - v){v - I) - n = 0, (8.21) dy dy dzi c- -fn + v = 0. (8.22) dy We have < ^ ~ 0(1), < 7~^ (5,e < 1. One can readily investigate these ordinary differential equations to find that the travelling wave speed is unique, giving a unique prediction for the speed of a nerve pulse in terms of biophysical parameters. Appendix A The phase plane Throughout this appendix we wih be concerned with systems of two coupled, first-order, autonomous, non-Unear ordinary differential equations. Disclaimer. This material should have been covered elsewhere (for example in your course on differential equations) and hence below is intended to review, rather than intro- duce and lecture this topic. We can represent solutions to the equations X{x,y), (A.l) y(x,y), (A.2) as trajectories (or "integral paths") in the phase plane, that is the {x^y) plane. Suppose, for the initial condition x{t = UniUai) = xq, y{t = tinuiai) = Vo we plot, in the {x,y) plane, the solution of (A.l): dx di dy dt We can do exactly the same for all the values of {UniUai -.x initial ■, ViniUai] ■, to build-up a graphical representation of the solutions to the equations (A.l) and (A.2) for many initial conditions. This plot is referred to as the "phase plane portrait". 81 Chapter A. The phase plane 82 A.l Properties of the phase plane portrait The gradient of the integral path through the point (2:0,2/0) is given by dy dy j dx Y{x,y) i^(a;o,yo) Key point 1. Note that if ^(xo, yo) = and X(xo, yo) 7^ then = 0, dy dx (2^0, yo) which corresponds to a horizontal line segments in the phase plane. Key point 2. If y(xo, yo) 7^ and X(xo, yo) = then dy dx 00 as (x,?/) (xo,2/o), which corresponds to a vertical line segment in the phase plane. (A.3) (A.4) (A.5) Key point 3. Assuming that either X(xo,yo) / or y(xo,yo) / 0, then two path integral curves do not cross at the point (xo, yo)- This is because under these circumstances dy/dx takes a unique value, i.e. the following is not possible: A. 2 EquiUbrium points Definition. A point in the phase plane where X(xo,yo) = y(2;o,yo) = is defined to be an equilibrium point, or equivalently, a stationary point. The reason for the above definition is because if (x,y) = (xo,yo) then both dx/dt and dy/dt are zero, and hence (x,y) do not change as t increases; hence x{t), y{t) remain at (xo,yo) for all time. Chapter A. The phase plane 83 Key point 1. Integral curves cannot cross at points whicli are not equilibrium points. Key point 2. If an integral path ends it must end on a stationary point. Key point 3. As we shall see below, equilibrium points are only approached as f — )■ oo or t — >■ —CO. However, what about the gradient of integral paths at (xo,yo)? We informally have , -r- = (A.6) dx ^ ' which is not uniquely defined — the value ultimately depends on the details of how quickly X{x,y) and Y{x,y) approach zero as (x, y) — )• {xo,yo), and this generally depends on the direction upon which {x,y) approaches {xo,yo). A. 2.1 Equilibrium points: further properties Suppose the equations (A.l) and (A. 2) have an equilibrium point at {xo,yo). Thus X{xo,yo) = Y{xQ,yo) = 0. To determine the behaviour of integral paths close to the equilibrium point we write X = Xo+X, y = yQ + y, (A.7) where it is assumed that x, y are sufficiently small to allow the approximations that we will make below. By Taylor expansion, we have dX dX X{x,y) = X{xo + x,yo + y) = X{xo,yo) + x—{xo,yo) + y—{xo,yo) + h..o.t. dX dX = x—{xo,yo)+y-^{xo,yo) + h.o.t., (A.8) using the fact X{xo,yo) = 0. Similarly, we have _dY _dY Y{x,y) = Y{xo + x,yo + y) = Y{xo,yo) + x—{xo,yo) +y-^{xo,yo) + h.o.t., dX dX = x—{xo,yo)+y—{xo,yo)+h..o.t. (A.9) ox oy Note that xq and yo are constant, and hence have zero time derivative. Hence, by use of Taylor expansions and neglecting higher orders (i.e. taking x, y sufficiently small), we can neglect terms of the order O and hence we can write equations (A.l) and (A. 2) in the form ^ = ( 1!'°' 'i 1!'°' 'i ] - Ju where u'il(']. (A.IO) Chapter A. The phase plane 84 Definition. The matrix J = ( ) (A.ll) is defined to be tfie Jacobian matrix at the equihbrium point (a::o,yo)- A. 3 Summary The key points thus far are as follows. 1. We have taken the full non-linear equation system, (A.l) and (A. 2), and expanded about one of its (possibly many) equilibrium points taken to be located at {xo,yo), using Taylor expansions of X{x,y), Y{x,y). 2. We assume that we are sufficiently close to (xo,?/o) to enable us to only consider linear terms of the order of {x — xq), {y — yo)- 3. In this way, we obtain a set of two coupled, linear, autonomous ordinary differential equations, i.e. equation (A. 10) above, which in principle we can solve! 4. This procedure is sometimes referred to as "a linearisation of equations (A.l) and (A. 2) about the point (a;o,yo)"- 5. In virtually all cases the behaviour of the linearised system is the same as the be- haviour of the full non-linear equations sufficiently close to the point (xq, yo)- In this respect one should note that the statement immediately above can be formulated more rigorously and proved for all the types of stationary points except: • centre type equilibrium points, i.e. case [3c] below; • the degenerate cases where Ai = and/or A2 = 0, which are briefly mentioned in item 2 on page (87). These stationary points can be considered non-examinable. The relevant theorem is "Hartmann's theorem", as discussed further in P. Glendin- ning. Stability, Instability and Chaos [4]. 6. However, one should also note that the solution of the linearised equations may behave substantially differently from the solutions of the full non-linear equations, (A.l) and (A. 2), sufficiently far from (xo,yo)- A. 4 Investigating solutions of the linearised equations We now have a set of two coupled, linear, autonomous ordinary differential equations, (A. 10). It is useful to look for a solution of the form U = UqC (A.12) Chapter A. The phase plane 85 for some constant, A. Substituting this into equation (A. 10) we obtain Xuoe^^ = Jtioe^* i.e. (J - XI) uq = 0. (A. 13) For a non-zero solution, we must have uq ^ (0, 0) and hence we require det(J-A/) = 0, (A.14) where / is the 2x2 identity matrix. This quadratic equation has two roots for A, denoted Ai, A2, which are possibly equal and possibly complex; these are, of course, the eigenvalues of J evaluated at the point (a;o,yo)- A.4.1 Case I Ai, A2 real, with Ai 7^ 0, A2 / 0, Ai / A2. Without loss of generality we take A2 > Ai below. We have two distinct, real eigenvalues. Let the corresponding eigenvectors be denoted by ei and 62- We thus have Jei = Aiei, Je2 = A2e2. (A. 15) We seek a solution of the form u = ^161 + ^262. (A. 16) Substituting this into equation (A. 10), we find, by comparing coefficients of ei and 62, that and hence Ai = Ai{t = 0)e^^\ ^2 = ^2(i = 0)e^2i_ ^^.18) Thus we have ^ u = Ai{t = 0)e^i*ei + A2{t = 0)e^^'e2, (A.19) y / which gives us a representation of the solution of (A. 10) for general initial conditions. This information is best displayed graphically, and we do so below according to the values of A2, Ai. Note. The equilibrium point i.e. {x,y) = (0,0) can only be reached either as i — ?• 00 or t — )• —00. 1. Ai < A2 < 0. The phase plot of the linearised equations in the {x, y) plane looks like one of the two possibilities in Figure A.l. Definition. An equilibrium point which results in this case is called a stable node, with the word "stable" referring to the fact that integral paths enter the node, i.e. the equilibrium point at (0,0). Chapter A. The phase plane 86 Figure A.l: Possible phase portraits of a stable node. The equilibrium point in eaeh case is denoted by the large dot. 2. A2 > Ai > 0. We still have u = Ai{t = 0)e-^i*ei + ^2(i = 0)e^^^e2. (A.20) However, the direction of the arrows is reversed as the signs of Ai, A2 are changed. The phase plane portraits are the same as in Figure A.l except the direction of the arrows is reversed. Definition. An equilibrium point which results in this case is called an unstable node, with the word "unstable" referring to the fact that integral paths leave the node, i.e. the equilibrium point at (0,0). 3. A2 > > Ai. Once more, we still have u = Ai{t = 0)e-^i*ei + A2{t = 0)e^^^e2, (A.21) but again the phase plane portrait is slightly different — see Figure A. 2. Definition. An equilibrium point which results in this case is called a saddle point. Definition. The two integral paths originating from the saddle point are some- times referred to as the unstable manifolds of the saddle point. Conversely, the integral paths tending to the saddle point are sometimes referred to as the sta- ble manifolds of the saddle point. This forms part of a nomenclature system commonly used in more advanced dynamical systems theory; see P. Glendinning, Stability, Instability and Chaos [4]. Chapter A. The phase plane 87 Figure A. 2: The phase portrait of a saddle point. The equilibrium point is denoted by the large dot. A.4.2 Case II A2, Ai real. One, or more, of the following also holds: A2 = Ai, Ai = 0, A2 = 0. (A.22) We typically will not encounter these degenerate cases in this course. We briefly note that behaviour of the full equations, (A.l), can be highly nontrivial when the linearisation reduces to these degenerate cases. Further details of such cases can be found in P. Glendin- ning, Stability, Instability and Chaos [4], which is on the reading list for this course. When Ai,A2 = 0, Hartmann's theorem doesn't hold. A.4.3 Case III A2, Ai complex. The complex eigenvalues of a real matrix always occur in complex conju- gate pairs. Thus we take, without loss of generality, X^ = a-ib = X*2, X2 = a + ib = Xl, (A.23) where a, b real, b ^ 0, and * denotes the complex conjugate. We also have two associated complex eigenvectors ei, 62, satisfying Jei = Aiei, Je2 = A2e2, (A. 24) which are complex conjugates of each other, i.e. ei = e^. Using the same idea as in Case I above, we have u = Ai{t = 0)e^i*ei + A2{t = Oje^^tg^, (A.25) though now, in general, Ai{t = 0), Ai, ei, A2(t = 0), A2, 62 are complex, and hence so is u. Chapter A. The phase plane 88 Restricting u to be real gives u = Ai{t = 0)e^i*ei + Al{t = 0)e^^*e2 = Ai{t = 0)e^i*ei + {Ai(t = 0)e^i*ei)*, (A.26) and this is real, as for any complex number z, we have z + z* E M. After some algebra this reduces to u = e''* [Mcos(6t) + Ksm{bt)] = e' where M = (Mi, M2)^ , K = {Ki,K2)^ are real, constant vectors, which can be expressed in terms of Ai{t = 0), A2{t = 0) and the components of the eigenvectors ei and 62- Equivalently, we have X = e"* [cos(6t)Mi + s\Ti{bt)Ki] , y = e"* [cos(6t)M2 + sm{ht)K2\ , (A.28) where Mi, M2, Ki, K2 are real constants. 1. a > 0. We have x, y are, overall, increasing exponentially but are oscillating too. For, example, with Ki = 0, Mi = 1 we have x = e"* cos{bt), which looks like: 2 S h"*') (A.27) Note that the overall growth of x is exponential at rate a. Thus, in general, the phase plane portrait looks like one of the examples shows in Figure A. 3. Note. The sense of the rotation, clockwise or anti-clockwise, is easily determined by calculating dy/dt when y = or dx/di when x = 0. Definition. An equilibrium point which results in the above, is called an un- stable spiral or, equivalently, an unstable focus. The word "unstable" refers to the fact that integral paths leave the equilibrium point. Chapter A. The phase plane 89 Figure A. 3: Possible pliase portraits of a focus. Tlie equilibrium point in each case is denoted by the large dot. 2. a < 0. This is the same as 1. except now the phase plane portrait arrows point towards the equilibrium point as x and y are exponentially decaying as time increases rather than exponentially growing. Definition. An equilibrium point which results in this case is called a stable spiral or, equivalently, a stable focus. The word "stable" refers to the fact that integral paths enter the equilibrium point. 3. a = 0. Thus we have A2 = —ib = — Ai, 6 7^ 0, 6 real, and X = [cos{bt)Mi + sm{bt)Ki] , y = [cos(6t)M2 + sui{bt)K2] , (A.29) where Mi, M2, Ki, K2 are constants. Note that K2X - Kiy = Lcos{bt), - M2X + Miy = L sm{bt) , (A. 30) where L = K2M1 - K1M2. Letting x* = K2X - Kiy and y* = -M2X + Miy, we have (x*)2 + (y*)2 = (A.31) i.e. a circle in the {x*,y*) plane, enclosing the origin, which is equivalent to, in general, a closed ellipse, in the (x,y) plane enclosing the origin. Note. As with 3. above, the sense of the rotation, clockwise or anti-clockwise, is easily determined by calculating dy/dt when y = or dx/dt when x = 0. Definition. An equilibrium point which results in this case, is called a centre. A centre is an example of a limit cycle. Chapter A. The phase plane 90 Definition. A limit cycle is an integral path which is closed (and which does not have any equilibrium points). A. 5 Linear stability Definition. An equilibrium point is linearly stable if the real parts of both eigenvalues Ai, A2 are negative. Prom the expressions for u above, for example u = Ai{t = 0)e^i*ei + A2{t = 0)e^''^e2, (A.32) when Ai, A2 real, we see that any perturbation away from the equilibrium decays back to the equilibrium point. Definition. An equilibrium point is linearly unstable if the real parts of at least one of the eigenvalues Ai, A2 is positive (and the other is non-zero). Other situations are in general governed by the non-linear behaviour of the full equations and we do not need to consider them here. A. 5.1 Technical point The behaviour of the linearised equations and the behaviour of the non-linear equations sufficiently close to the equilibrium point are guaranteed to be the same for any of the equilibrium points [I 1-3], [III 1,2] or [II] with Ai = A2 7^ 0. All these equilibrium points are such that i?e(Ai) Re{\2) 7^ 0. This is the essence of Hartmann's theorem. This guarantee does not hold for centres, or the equilibrium points described in [II] with A1A2 = 0. The underlying reasons for this are as follows. • First, note from the above that integral paths which meet the equilibrium point can either grow/decay at exponential rate i?e(Ai), or exponential rate Re{\2), or consist of the sum of two such terms. Second, note that in the above we took a Taylor expansion. Including higher order terms in this Taylor expansion can lead to a small correction for the rate of exponential decay towards or growth away from the stationary point exhibited by the integral paths. These corrections tend to zero as one approaches the equilibrium point. • Consider the centre equilibrium point, which has Re{\i) = Re{\2) = 0, and an exponential growth/decay of zero. If the corrections arising from the Taylor series are always positive, the exponential growth/decay rate of all integral paths sufficiently near the stationary point is always (slightly) positive. Hence these integral paths grow exponentially away from the stationary point. However, b is non-zero, so x and y are still oscillating. Hence one has the non-linear equations behave like a stable focus. Chapter A. The phase plane 91 • If Re{Xi) , Re{X2) 7^ 0, then for all integral paths reaching the stationary point, the above mentioned corrections, sufficiently close to the equilibrium point, are negligi- ble, e.g. they cannot change exponential growth into exponential decay or vice-versa. This allows one to show that stationary points with i?e(Ai) Re{\2) 7^ are guar- anteed to have the same behaviour for the linearised and the non-linear equations sufficiently close to the equilibrium point. A. 6 Summary We will typically only encounter stationary points [I 1-3], [III 1-3]. Of these stationary points, all but centres exhibit the same behaviour for the linearised and the non-linear equations sufficiently close to the equilibrium point as plotted above and in D. W. Jordan and P. Smith, Mathematical Techniques [6]. Bibliography [1] N. F. Britton. Essential Mathematical Biology. Springer Undergraduate Mathematics Series. Springer, 2005. [2] L. Edelstein-Keshet. Mathematical Models in Biology. SIAM Classics in Applied Mathematics, 2005. [3] A. Gierer and H. Meinhardt. A theory of biological pattern formation. Kybernetik, 12:30-39, 1972. [4] P. Glendinning. Stability, Instability and Chaos: An Introduction to the Theory of Nonlinear Differential Equations. Cambridge Texts in Applied Mathematics, 1999. [5] T. Hillen and K. Painter. A user's guide to PDE models for chemotaxis. J. Math. Biol, 58(1):183-217, 2009. [6] D. W. Jordan and P. Smith. Mathematical Techniques: An Introduction for the En- gineering, Physical and Mathematical Sciences. Oxford University Press, 3rd edition, 2002. [7] J. P. Keener and J. Sneyd. Mathematical Physiology, volume 8 of Interdisciplinary Applied Mathematics. Springer, New York, 1st edition, 1998. [8] J. D. Murray. Mathematical Biology I: An Introduction, volume I. Springer- Verlag, 3rd edition, 2003. [9] J. D. Murray. Mathematical Biology II: Spatial Models and Biochemical Applications, volume II. Springer- Verlag, 3rd edition, 2003. [10] A. Okubo, P. K. Maini, M. H. Wilhamson, and J. D. Murray. On the spread of the grey squirrel in Great Britain. Proc. R. Soc. Lond. B, 238(1291):113-125, 1989. [11] L. E. Reichl. A Modern Course in Statistical Physics. Wiley- VCH, 3rd edition, 2009. [12] T. F. Weiss. Cellular Biophysics, volume 2. MIT Press, 1996. 92