Mathematical Biology Lecture Notes, Lecture notes of Biology

Download Mathematical Biology Lecture Notes and more Lecture notes Biology in PDF only on Docsity!

Further Mathematical Biology

Lecture Notes

Prof Helen M Byrne

Michaelmas Term 2017

Contents

• 1 Introduction
  • 1.1 References
• 2 Enzyme kinetics
  • 2.1 Introduction
  • 2.2 The Law of Mass Action
  • 2.3 Michaelis-Menten kinetics
    • 2.3.1 Non-dimensionalisation
    • 2.3.2 Singular perturbation investigation
  • 2.4 More complex systems
    • 2.4.1 Several enzyme reactions and the pseudo-steady state hypothesis
    • 2.4.2 Allosteric enzymes
    • 2.4.3 Autocatalysis and activator-inhibitor systems
• 3 Ions and Excitable Systems
  • 3.1 Introduction
    • 3.1.1 Background
    • 3.1.2 Basic Concepts
  • 3.2 Channel gating
    • 3.2.1 Simple Gates
    • 3.2.2 Multiple gates
    • 3.2.3 Non-identical gates
  • 3.3 The Fitzhugh Nagumo equations
    • 3.3.1 Deducing the FitzHugh Nagumo Equations
    • 3.3.2 A brief analysis of the Fitzhugh Nagumo equations
  • 3.4 Modelling nerve signal propagation - NON-EXAMINABLE
    • 3.4.1 The cable model
• 4 Introduction to spatial variation
  • 4.1 Derivation of the reaction-diffusion equations
  • 4.2 Chemotaxis
  • 4.3 Positional Information and Pattern Formation
• CONTENTS
• 5 Travelling waves
  • 5.1 Fisher’s equation: an investigation
    • 5.1.1 Key points
    • 5.1.2 Existence and the phase plane
    • 5.1.3 Relation between the travelling wave speed and initial conditions
  • 5.2 Models of epidemics
    • 5.2.1 The SIR model (revision from Short Course)
    • 5.2.2 An SIR model with spatial heterogeneity
• 6 Pattern formation
  • 6.1 Minimum domains for spatial structure
    • 6.1.1 Domain size
  • 6.2 Diffusion-driven instability
    • 6.2.1 Linear analysis
  • 6.3 Detailed study of the conditions for a Turing instability
    • 6.3.1 Stability without diffusion
    • 6.3.2 Instability with diffusion
    • 6.3.3 Summary
    • 6.3.4 The threshold of a Turing instability.
  • 6.4 Extended example
    • 6.4.1 The influence of domain size
  • 6.5 Extended example
• 7 Domain Growth in Biology
  • 7.1 Simplest model: 1D solid tumour growth
  • 7.2 Model Reduction: Nondimensionalise
  • 7.3 Cell death at low nutrient concentration
  • 7.4 Revised model: proliferation and necrosis
  • 7.5 Summary
• 8 From Discrete to Continuum Models
  • 8.1 Introduction
  • 8.2 Biased Random Walks and Advection-Diffusion Equations
    • 8.2.1 Boundary Conditions
    • 8.2.2 Model Extensions
  • 8.3 Age-Structured Populations
    • 8.3.1 Structured Models for Proliferating Cells
    • 8.3.2 Age-dependent epidemic models (not examinable)
• A The phase plane
  • A.1 Properties of the phase plane portrait
  • A.2 Equilibrium points
    • A.2.1 Equilibrium points: further properties
  • A.3 Summary
  • A.4 Investigating solutions of the linearised equations
• CONTENTS - A.4.1 Case I - A.4.2 Case II - A.4.3 Case III
  • A.5 Linear stability
    • A.5.1 Technical point
  • A.6 Summary

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 focus initially on systems where spatial variation is either absent or, at least, not important. In such cases only temporal evolution needs to be described, typically via ordinary differential equations.
• We are inevitably confronted with systems of ordinary differential equations, and thus we will study analytical techniques for extracting information from such equa- tions.
• We will then consider systems where there is explicit spatial variation. The resulting models must also incorporate spatial effects.
• In ecological and biological applications the main physical phenomenon governing spatial variation is typically, but not exclusively, diffusion. Thus we invariably con- sider parabolic partial differential equations. Mathematical techniques will be de- veloped 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: 1. enzyme-substrate dynamics and other biochemical reactions; 2. Trans-membrane ion channels and nerve pulses; 3. epidemics; 4. the propagation of an advantageous gene through a population; 5. biological pattern formation mechanisms; 6. chemotaxis; 7. tumour growth.

Acknowledgements: these lecture notes build on material originally developed by Pro- fessors Philip Maini, Ruth Baker, Eamonn Gaffney, Jon Chapman and Andrew Fowler. I am extremely grateful to them for allowing me to re-use and extend their notes.

5

Chapter 2

Enzyme kinetics

2.1 Introduction

Figure 2.1: How do enzymes work? An enzyme has an active site where the substrates and enzyme fit together so that the substrates react. After the reaction, the products are released and the enzyme assumes its original shape.

Biochemical reactions are extremely important for correct biological function. For example, biochemical reactions are involved in:

• Metabolism and its control;
• Immunological responses;
• Cell-signalling processes.

Biochemical processes are often controlled by enzymes. Enzymes are proteins that catalyse biochemical reactions by lowering activation energy. Even when present in very small amounts, enzymes can have a dramatic effect on a system (see Figure 2.3). Table 2. illustrates how effective enzymes can be at accelerating reactions in biochemical systems.

Chapter 2. Enzyme kinetics 8

Figure 2.2: The enzyme catecholase catalyses a reaction between the molecule catechol and oxygen. The product of this reaction is polyphenol, the brown substance that accumulate when apples are exposed to air.

Enzyme Substrate Product Rate without Rate with Accel. enzyme enzyme due to enzyme Hexokinase Glucose Glucose < 0. 0000001 1300 > 13 billion 6-Phosphate Phosphorylase – – < 0. 000000005 1600 > 320 billion Alcohol Ethanol Acetaldehyde < 0. 000006 2700 > 450 million Dehydrogenase Creatine Creatine Creatine < 0. 003 40 > 13 , 000 Kinase Phosphate

Table 2.1: Examples illustrating the impact that enzymes can have on reaction rates.

In this chapter we will focus on enzyme kinetics. These can be thought of as a special case of an interacting species model (for details, see lecture notes on short course ’Mathematical Modelling in Biology’). In all cases we will neglect spatial variation.

References.

• J. D. Murray, Mathematical Biology, 3rd edition, Volume I, Chapter 6 [?].
• J. P. Keener and J. Sneyd, Mathematical Physiology, Chapter 1 [?].

2.2 The Law of Mass Action

Throughout this chapter, we will consider reactions involving m chemical species C 1 ,... , Cm.

• The concentration of Ci, denoted ci, is defined to be the number of molecules of Ci per unit volume.

Chapter 2. Enzyme kinetics 10

Example: Stoichiometry. Suppose m molecules of A react reversibly with n molecules of B to create C:

mA + nB

k 1

k− 1

C

Then the Law of Mass Action takes the form da dt =^ −mk^1 a

mbn (^) + mk− 1 c db dt dc =^ −nk^1 ambn^ +^ nk−^1 c dt =^ k^1 ambn^ −^ k−^1 c

as m molecules of A and n molecules of B must collide to produce one molecule of C.

Note: Mass conservation supplies

a + mc = constant, b + nc = constant.

2.3 Michaelis-Menten kinetics

Michaelis-Menten kinetics approximately describe the dynamics of a number of enzyme systems. The reactions are

S + E

k 1 F^ GGGGGGGGGGBGG k− 1

C, (2.4)

C

k 2 GGGGGGA P + E, (2.5)

where C represents the complex SE, and s, e, p and c denote the concentrations of S, E, P and C respectively. From the Law of Mass Action, we can derive the following ordinary differential equations:

ds dt

= −k 1 se + k− 1 c, (2.6) dc dt =^ k^1 se^ −^ k−^1 c^ −^ k^2 c,^ (2.7) de dt

= −k 1 se + k− 1 c + k 2 c, (2.8) dp dt

= k 2 c. (2.9)

Note that the equation for p decouples and, hence, we can neglect it initially.

The initial conditions are:

s(0) = s 0 , e(0) = e 0  s 0 , c(0) = 0, p(0) = 0. (2.10)

Chapter 2. Enzyme kinetics 11

Key Point. In systems described by the Law of Mass Action, linear combinations of the variables are often conserved. In this example we have

d dt

(e + c) = 0 ⇒ e = e 0 − c, (2.11)

and, hence, the equations simplify to:

ds dt

= −k 1 (e 0 − c)s + k− 1 c, (2.12) dc dt =^ k^1 (e^0 −^ c)s^ −^ (k−^1 +^ k^2 )c,^ (2.13)

with the determination of p readily achievable once the dynamics of s and c are known.

2.3.1 Non-dimensionalisation

We non-dimensionalise as follows:

τ = k 1 e 0 t, u =

s s 0 ,^ v^ =^

c e 0 ,^ λ^ =^

k 2 k 1 s 0 ,^ ^

def = e 0 s 0 ^1 ,^ K^

def = k− 1 +^ k 2 k 1 s 0 ,^ (2.14)

which yields

u′^ = −u + (u + K − λ)v, (2.15) v′^ = u − (u + K)v, (2.16)

where u(0) = 1, v(0) = 0 and   1. Normally  ∼ 10 −^6. Setting  = 0 yields

v =

u u + K ,^ (2.17)

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 = 0 where v′^  O(1). Indeed for the initial conditions given we find v′(0) ∼ O(1/), with u(0), v(0) ≤ O(1). This gives us the scaling we need for a singular perturbation investigation.

2.3.2 Singular perturbation investigation

We consider σ = τ , (2.18)

with

u(τ, ) = u˜(σ, ) = ˜u 0 (σ) + u˜ 1 (σ) +... , (2.19) v(τ, ) = ˜v(σ, ) = ˜v 0 (σ) + ˜v 1 (σ) +.... (2.20)

Proceeding in the usual way, we find that ˜u 0 , ˜v 0 satisfy

d˜u 0 dσ

= 0 ⇒ u˜ 0 = constant = 1, (2.21)

and d˜v 0 dσ = ˜u^0 −^ (1 +^ K)˜v^0 = 1^ −^ (1 +^ K) ˜v^0 ⇒^ v˜^0 =

1 − e−(1+K)σ 1 + K ,^ (2.22)

Chapter 2. Enzyme kinetics 13

Often the initial, fast, transient is not seen or modelled: we consider only the outer equations, with a suitably adjusted initial condition (ultimately determined from consis- tency/matching with the inner solution). In particular, we often use Michaelis-Menten kinetics where the equations are simply:

du dt

= − λu u + K

with u(0) = 1 and v = u u + K

Definition. We have, approximately, that dv/dτ ' 0 using Michaelis-Menten kinet- ics. When the time derivative is fast, i.e. of the form

 d dvτ = g(u, v), (2.32)

where   1 and g(u, v) ∼ O(1), taking the temporal dynamics to be trivial, dv dτ

is known as the pseudo-steady state hypothesis. This is a common assumption in the literature. We have seen its validity for enzyme kinetics, at least away from the inner region.

Note. While the Michaelis-Menten kinetics derived above are a very useful approxima- tion, they hinge on the validity of the Law of Mass Action. Even in simple biological systems the Law of Mass Action may break down. 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.

2.4 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 [?], for more details.

2.4.1 Several enzyme reactions and the pseudo-steady state hypothesis

We can have multiple enzymes. In general the system of equations reduces to

u′^ = f (u, v 1 ,... , vn), (2.34) iv′ i = gi(u, v 1 ,... , vn), (2.35)

for i ∈ { 1 ,... , n}, while the pseudo-steady state hypothesis gives a single ordinary differ- ential equation u′^ = f (u, v 1 (u),... , vn(u)), (2.36)

Chapter 2. Enzyme kinetics 14

where v 1 (u),... , vn(u) are the appropriate roots of the equations

gi(u, v 1 ,... , vn) = 0, i ∈ { 1 ,... , n}. (2.37)

Exercise. Consider an enzymatic reaction in which an enzyme can be activated or inactivated by a chemical X as follows:

E + X

k 1

k− 1

E 1 , E 1 + X

k 2

k− 2

E 2 , E 1 + S

k 3 → P + Q + E.

Suppose further that X is supplied at a constant rate, and removed at a rate propor- tional to its concentration.

1. Write down differential equations for the evolution of E, E 1 , E 2 , X and S.
2. Show that E + E 1 + E 2 is a conserved quantity, E∗^ say.
3. Nondimensionalise the system, scaling E, E 1 and E 2 with E∗, X and S with X 0 = X(0), and time with 1/(k 1 E∗). Assuming that δ = E∗/X 0  1, use the resulting “quasi-steady” equations for the dimensionless quantities e, e 1 , e 2 to solve for these variables in terms of x and s, and hence obtain the following system of two ODEs for x and s only: dx dτ =^ α^0 −^ ν^4 x^ −^

κ 3 xs μ 1 + κ 3 s + x + κ 2 x^2 /μ 2 ,

ds dτ =^ −^

κ 3 xs μ 1 + κ 3 s + x + κ 2 x^2 /μ 2. Identify all parameters and variables in these equations.

2.4.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

k 1 F^ GGGGGGGGGGBGG k− 1

C 1

k 2 GGGGGGA P + E (2.38)

S + C 1

k 3 F^ GGGGGGGGGGBGG k− 3

C 2

k 4 GGGGGGA C 1 + P. (2.39)

Further details on the investigation of such systems can be found in J. D. Murray, Mathe- matical Biology Volume I [?], and J. P. Keener and J. Sneyd, Mathematical Physiology [?].

2.4.3 Autocatalysis and activator-inhibitor systems

Here a molecule catalyses its own production. The simplest example is the reaction scheme

A + B →k 2 B, (2.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.

Chapter 3

Ions and Excitable Systems

3.1 Introduction

3.1.1 Background

The cell membrane is a phospholipid bilayer separating the cell interior (the cytoplasm) from the extracellular environment. The membrane contains numerous proteins, and is approximately 7.5nm thick. The most important property of the cell membrane is its selective permeability: it allows the passage of some molecules but restricts the passage of others, thereby regulating the passage of materials into and out of the cell. Many substances penetrate the cell membrane at rates reflected by their diffusive behaviour in a pure phospholipid bilayer. However, certain molecules and ions such as glucose, amino acids and Na+^ pass through cell membranes much more rapidly, indicating that the membrane proteins selectively facilitate transport.

Figure 3.1: A schematic of the phospholipid membrane double layer, with a gating protein in one of two configurations, Ce and Ci, spanning the membrane, as part of a passive, carrier- mediated transport system.

The membrane contains water-filled pores with diameters of about 0.8nm, and protein- lined pores, called channels or gates, which allow the passage of specific molecules. Both the intracellular and extracellular environments comprise (among other things) a dilute aqueous solution of dissolved salts, mainly NaCl and KCl, which dissociate into Na+, K+ and Cl−^ ions. The cell membrane acts as a barrier to the free flow of these ions and to the flow of water.

Chapter 3. Ions and Excitable Systems 17

Figure 3.2: Representation of an exchange pump which actively transports across the cell membrane.

The mechanisms that facilitate transport across the cellular membrane can be divided into active and passive processes. Active processes requires energy expenditure, while passive processes result solely from the random motion of molecules, for example, diffusion. Action potentials, or ’nerve impulses’, are brief changes in the membrane potential of a cell produced by the flow of ionic current across the cell membrane. They enable commu- nication by many cell types, including neurons, cardiac and muscle cells. In section 3.3 we will study the Fitzhugh Nagumo equations which describe nerve impulses in axons. First, we outline the basic physical concepts needed to study ion channels and nerve impulses, and introduce basic mathematical models for ionic currents and channel gating.

Note: Much of the material covered in this course is extended and studied in greater detail in the Part C Course Mathematical Physiology.

3.1.2 Basic Concepts

First, we note that:

• numerous fundamental particles, ions and molecules have an electric charge, e.g. the electron, e−, and the sodium ion, N a+;
• 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.

Chapter 3. Ions and Excitable Systems 19

charges stored on the two plates. The capacitance of the plates, C, is defined to be

C = Qeqm V

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.

If I is the ionic current out of the cell (the rate of flow of positive charges outwards),then the stored charge changes according to

I = −

dQ dt.^ (3.2)

Thus, assuming the capacitance is constant

C dV dt

+ I = 0. (3.3)

This equation is the basis for much theoretical electrophysiology. The difference between various models relates to the expression used for the ionic current I. The simplest models assume linear dependent of I on V (as in Ohms law). For a single ion S, with Nernst potential VS , this gives an ionic current

IS = gS (V − VS ),

where the constant gS is the ion-specific membrane conductance, since the current must vanish when V = VS. If more than one ion is present, then the currents from the different ions are added together to produce the total ionic current.

3.2 Channel gating

3.2.1 Simple Gates

In practice, gS is not constant: it depends on V and time t. One explanation for this is that the channels are not always open (they may be open or closed), and the transition rates between open and closed states depends on the potential difference, V. The membrane conductance may then be written as ngS , where gS is the constant conductance that would result if all channels were open, and n is the proportion of open channels. For a generic ion, let n be the proportion of open ion channels. Denoting the open channels by O and the closed channels by C, the reaction scheme is simply

C

α(V )

β(V )

O (3.4)

Chapter 3. Ions and Excitable Systems 20

Figure 3.3: Schematic diagram of channel gating.

where α(V ) and β(V ) represent voltage dependent rates of switching between the closed and open states. Using the law of mass action we obtain

dn dt =^ α(V^ )(1^ −^ n)^ −^ β(V^ )n,^ (3.5) or, equivalently,

τn(V ) dn dt

= n∞(V ) − n, (3.6)

where n∞(V ) = α/(α + β) is the equilibrium value of n and τn(V ) = 1/(α + β) is the timescale for approach to this equilibrium (Note: both n∞ and τn can be determined experimentally).

3.2.2 Multiple gates

The simple model presented above can be generalised to channels with multiple identical subunits, each of which can be in either the open or closed state. We start by assuming that the channel consists of two “gates”, which may both exist in open or closed states. The ion channel is open only if both “gates” are open; the ion channel is closed if any one gate within the ion channel is closed. If we denote by Si (i ∈ { 0 , 1 , 2 }) the proportion of channels with exactly i gates open, then we have

S 0 + S 1 + S 2 = 1 (3.7)

and the reaction scheme

S 0

2 α(V )

β(V )

S 1

α(V )

2 β(V )

S 2. (3.8)

The 2s arise because there are two possible states with one gate open and one gate closed. Since each gate is identical we lump these two states into one variable S 1. Using mass action kinetics gives

dS 0 dt =^ β(V^ )(1^ −^ S^0 −^ S^2 )^ −^2 α(V^ )S^0 ,^ (3.9) dS 2 dt

= α(V )(1 − S 0 − S 2 ) − 2 β(V )S 2 , (3.10)

Figure 3.4: Schematic diagram of two identical gate units.