Mathematical Biology Lecture notes for MATH 4333 (formerly MATH 365)
Jefrey R. Chasnov
The Hong Hong Kong Univer Universit sity y of Science Science and Technolo Technology gy
The Hong Kong University University o Science Science and Technolog Technology y Department o Mathematics Clear Water Bay, Kowloon Hong Kong
c 2009–2012 Copyright ○ 2009–2012 by Jefrey Robert Chasnov Chasnov This work is licensed under the Creative Commons Attribution 3.0 Hong Kong License. To view a copy o this license, license, visit http://crea http://creative tivecommon commons.org/li s.org/licenses/ censes/by/3 by/3.0/hk .0/hk/ / or send a letter to Creative Commons, 171 Second Street, Suite 300, San Francisco, Caliornia, Caliornia, 94105, 94105, USA.
The Hong Kong University University o Science Science and Technolog Technology y Department o Mathematics Clear Water Bay, Kowloon Hong Kong
c 2009–2012 Copyright ○ 2009–2012 by Jefrey Robert Chasnov Chasnov This work is licensed under the Creative Commons Attribution 3.0 Hong Kong License. To view a copy o this license, license, visit http://crea http://creative tivecommon commons.org/li s.org/licenses/ censes/by/3 by/3.0/hk .0/hk/ / or send a letter to Creative Commons, 171 Second Street, Suite 300, San Francisco, Caliornia, Caliornia, 94105, 94105, USA.
Preace What ollows are my lecture notes or Math 4333: Mathematical Mathematical Biology , taught at the Hong Kong Universit University y o Science and Technolog echnology y. This applied mathematics course is primarily or nal year mathematics major and minor students. Other students are also welcome to enroll, but must have the necessary mathematical skills. My main emphasis is on mathematical modeling, with biology the sole application plication area. I assume assume that students students have no knowledge knowledge o biology, biology, but I hope that they will learn a substanti substantial al amount during the course. course. Students Students are required to know diferential equations and linear algebra, and this usually means having having taken two two courses in these subjects. I also touch on topics topics in stochasstochastic modeling, which requires requires some knowledge knowledge o probability probability. A ull course on probability, however, is not a prerequisite though it might be helpul. Biology, as is usually taught, requires memorizing a wide selection o acts and remembering remembering them or exams, exams, sometimes sometimes orgetting orgetting them soon ater. ater. For students exposed to biology in secondary school, my course may seem like a diferent diferent subject. The ability ability to model problems problems using mathematics mathematics requires almost no rote memorization, but it does require a deep understanding o basic principles principles and a wide range o mathematical mathematical techniques. techniques. Biology Biology ofers a rich variety o topics that are amenable to mathematical modeling, and I have chosen specic topics that I have ound to be the most interesting. I, as a UST student, you have not yet decided i you will take my course, please browse these lecture notes to see i you are interested in these topics. Other web surers are welcome to download these notes rom http://www.math.ust.hk/˜machas/mathematical-biology.pd and to use them reely or teaching and learning. I welcome any comments, suggestions, or corrections sent to me by email (
[email protected]). Although most o the material in my notes can be ound elsewhere, I hope that some o it will be considered to be original. Jeffrey R. Chasnov Chasnov
Hong Kong May 2009
iii
Contents 1 Population Dynamics 1.1 The Malthusian growth model . . 1.2 The Logistic equation . . . . . . 1.3 A model o species competition . 1.4 The Lotka-Volterra predator-prey
. . . . . . . . . . . . model
. . . .
2 Age-structured Populations 2.1 Fibonacci’s rabbits . . . . . . . . . . . . . 2.2 The golden ratio Φ . . . . . . . . . . . . . 2.3 The Fibonacci numbers in a sunower . . 2.4 Rabbits are an age-structured population 2.5 Discrete age-structured populations . . . . 2.6 Continuous age-structured populations . . 2.7 The brood size o a hermaphroditic worm
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
3 Stochastic Population Growth 3.1 A stochastic model o population growth . . . . . . . . . 3.2 Asymptotics o large initial populations . . . . . . . . . 3.2.1 Derivation o the deterministic model . . . . . . 3.2.2 Derivation o the normal probability distribution 3.3 Simulation o population growth . . . . . . . . . . . . . 4 Inectious Disease Modeling 4.1 The SI model . . . . . . . . . . . 4.2 The SIS model . . . . . . . . . . 4.3 The SIR epidemic disease model 4.4 Vaccination . . . . . . . . . . . . 4.5 The SIR endemic disease model . 4.6 Evolution o virulence . . . . . .
. . . . . .
5 Population Genetics 5.1 Haploid genetics . . . . . . . . . . 5.1.1 Spread o a avored allele . 5.1.2 Mutation-selection balance 5.2 Diploid genetics . . . . . . . . . . . 5.2.1 Sexual reproduction . . . . 5.2.2 Spread o a avored allele . 5.2.3 Mutation-selection balance v
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . .
1 1 2 6 8
. . . . . . .
15 15 17 18 21 23 26 29
. . . . .
37 37 41 42 45 47
. . . . . .
51 51 52 54 56 58 59
. . . . . . .
61 62 63 64 66 67 70 72
vi
CONTENTS
5.2.4 Heterosis . . . . . . . . 5.3 Frequency-dependent selection 5.4 Linkage equilibrium . . . . . . 5.5 Random genetic drit . . . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
74 75 78 82
6 Biochemical Reactions 6.1 The law o mass action 6.2 Enzyme kinetics . . . 6.3 Competitive inhibition 6.4 Allosteric inhibition . 6.5 Cooperativity . . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
89 89 92 93 96 99
7 Sequence Alignment 7.1 DNA . . . . . . . . . . 7.2 Brute orce alignment 7.3 Dynamic programming 7.4 Gaps . . . . . . . . . . 7.5 Local alignments . . . 7.6 Sotware . . . . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
103 104 105 108 111 112 113
Chapter 1
Population Dynamics Populations grow in size when the birth rate exceeds the death rate. Thomas Malthus, in An Essay on the Principle o Population (1798), used unchecked population growth to amously predict a global amine unless governments regulated amily size—an idea later echoed by Mainland China’s one-child policy. The reading o Malthus is said by Charles Darwin in his autobiography to have inspired his discovery o what is now the cornerstone o modern biology: the principle o evolution by natural selection. The Malthusian growth model is the granddaddy o all population models, and we begin this chapter with a simple derivation o the amous exponential growth law. Unchecked exponential growth obviously does not occur in nature, and population growth rates may be regulated by limited ood or other environmental resources, and by competition among individuals within a species or across species. We will develop models or three types o regulation. The rst model is the well-known logistic equation, a model that will also make an appearance in subsequent chapters. The second model is an extension o the logistic model to species competition. And the third model is the amous Lotka-Volterra predator-prey equations. Because all these mathematical models are nonlinear diferential equations, mathematical methods to analyze such equations will be developed.
1.1
The Malthusian growth model
Let () be the number o individuals in a population at time , and let and be the average per capita birth rate and death rate, respectively. In a short time Δ, the number o births in the population is Δ , and the number o deaths is Δ . An equation or at time + Δ is then determined to be ( + Δ) = () + Δ ()
− Δ (),
which can be rearranged to ( + Δ) Δ and as Δ
→ 0,
− () = ( − ) ();
= ( 1
− ).
2
CHAPTER 1. POPULATION DYNAMICS
With an initial population size o 0 , and with = or = () grows exponentially:
− positive, the solution
() = 0 . With population size replaced by the amount o money in a bank, the exponential growth law also describes the growth o an account under continuous compounding with interest rate .
1.2
The Logistic equation
The exponential growth law or population size is unrealistic over long times. Eventually, growth will be checked by the over-consumption o resources. We assume that the environment has an intrinsic carrying capacity , and populations larger than this size experience heightened death rates. To model population growth with an environmental carrying capacity , we look or a nonlinear equation o the orm = ( ), where ( ) provides a model or environmental regulation. This unction should satisy (0) = 1 (the population grows exponentially with growth rate when is small), ( ) = 0 (the population stops growing at the carrying capacity), and ( ) < 0 when > (the population decays when it is larger than the carrying capacity). The simplest unction ( ) satisying these conditions is linear and given by ( ) = 1 / . The resulting model is the well-known logistic equation,
−
= (1 / ), (1.1) an important model or many processes besides bounded population growth. Although (1.1) is a nonlinear equation, an analytical solution can be ound by separating the variables. Beore we embark on this algebra, we rst illustrate some basic concepts used in analyzing nonlinear diferential equations. Fixed points, also called equilibria, o a diferential equation such as ( 1.1) are dened as the values o where / = 0. Here, we see that the xed points o (1.1) are = 0 and = . I the initial value o is at one o these xed points, then will remain xed there or all time. Fixed points, however, can be stable or unstable. A xed point is stable i a small perturbation rom the xed point decays to zero so that the solution returns to the xed point. Likewise, a xed point is unstable i a small perturbation grows exponentially so that the solution moves away rom the xed point. Calculation o stability by means o small perturbations is called linear stability analysis . For example, consider the general one-dimensional diferential equation (using the notation ˙ = /) ˙ = (), (1.2)
−
with * a xed point o the equation, that is (* ) = 0. To determine analytically i * is a stable or unstable xed point, we perturb the solution. Let us write our solution = () in the orm () = * + (),
(1.3)
3
1.2. THE LOGISTIC EQUATION
unstable
) x ( f
unstable
stable
0
x
Figure 1.1: Determining one-dimensional stability using a graphical approach. where initially (0) is small but diferent rom zero. Substituting (1.3) into (1.2), we obtain ˙ = (* + ) = (* ) + ′ (* ) + . . . = ′ (* ) + . . . , where the second equality uses a Taylor series expansion o () about * and the third equality uses (* ) = 0. I ′ (* ) = 0, we can neglect higher-order terms in or small times, and integrating we have
̸
′
() = (0) (
*
)
.
The perturbation () to the xed point * goes to zero as ′ (* ) < 0. Thereore, the stability condition on * is *
is
{
a stable xed point i an unstable xed point i
→ ∞ provided
′ (* ) < 0, ′ (* ) > 0.
Another equivalent but sometimes simpler approach to analyzing the stability o the xed points o a one-dimensional nonlinear equation such as (1.2) is to plot () versus . We show a generic example in Fig. 1.1. The xed points are the x-intercepts o the graph. Directional arrows on the x-axis can be drawn based on the sign o (). I () < 0, then the arrow points to the let; i () > 0, then the arrow points to the right. The arrows show the direction o motion or a particle at position satisying ˙ = (). As illustrated in Fig. 1.1, xed points with arrows on both sides pointing in are stable, and xed points with arrows on both sides pointing out are unstable. In the logistic equation (1.1), the xed points are * = 0, . A sketch o ( ) = (1 / ) versus , with , > 0 in Fig. 1.2 immediately shows that * = 0 is an unstable xed point and * = is a stable xed point. The
−
4
CHAPTER 1. POPULATION DYNAMICS
) N ( F
unstable
stable
0
0
K N
Figure 1.2: Determining stability o the xed points o the logistic equation. analytical approach computes ′ ( ) = (1 2/ ), so that ′ (0) = > 0 and ′ ( ) = < 0. Again we conclude that * = 0 is unstable and * = is stable. We now solve the logistic equation analytically. Although this relatively simple equation can be solved as is, we rst nondimensionalize to illustrate this very important technique that will later prove to be most useul. Perhaps here one can guess the appropriate unit o time to be 1/ and the appropriate unit o population size to be . However, we preer to demonstrate a more general technique that may be useully applied to equations or which the appropriate dimensionless variables are dicult to guess. We begin by nondimensionalizing time and population size:
−
−
= /* ,
= / * ,
˙ is computed where * and * are unknown dimensional units. The derivative as ( * ) * = = . * Thereore, the logistic equation (1.1) becomes
−
= * 1
*
,
which assumes the simplest orm with the choices * = 1/ and * = . Thereore, our dimensionless variables are = ,
= /,
and the logistic equation, in dimensionless orm, becomes = (1
− ) ,
(1.4)
5
1.2. THE LOGISTIC EQUATION
with the dimensionless initial condition (0) = 0 = 0 / , where 0 is the initial population size. Note that the dimensionless logistic equation (1.4) has no ree parameters, while the dimensional orm o the equation (1.1) contains and . Reduction in the number o ree parameters (here, two: and ) by the number o independent units (here, also two: time and population size) is a general eature o nondimensionalization. The theoretical result is known as the Buckingham Pi Theorem . Reducing the number o ree parameters in a problem to the absolute minimum is especially important beore proceeding to a numerical solution. The parameter space that must be explored may be substantially reduced. Solving the dimensionless logistic equation (1.4) can proceed by separating the variables. Separating and integrating rom = 0 to and 0 to yields
0
= (1 )
−
.
0
The integral on the let-hand-side can be perormed using the method o partial ractions: 1 (1
−
+ ) 1 + ( ) = ; (1 ) =
− − −
and by equating the coecients o the numerators proportional to 0 and 1 , we nd that = 1 and = 1. Thereore,
0
+ ) 0 0 (1 1 = ln ln 0 1 0 (1 0 ) = ln 0 (1 ) = .
= (1 )
−
− −− − −
−
Solving or , we rst exponentiate both sides and then isolate : (1 0 ) = , or (1 0 (1 )
− −
or (1
− 0) = 0 − 0 ,
− 0 + 0 ) = 0 , or
=
0 0 + (1
−
− 0)
.
Returning to the dimensional variables, we nally have () =
0 . 0 / + (1 0 / )−
−
(1.5)
There are several ways to write the nal result given by (1.5). The presentation o a mathematical result requires a good aesthetic sense and is an important element o mathematical technique. When deciding how to write (1.5), I considered i it was easy to observe the ollowing limiting results: (1) (0) = 0 ; (2) lim→∞ () = ; and (3) lim →∞ () = 0 exp().
6
CHAPTER 1. POPULATION DYNAMICS Logistic Equation 1.4 1.2 1 0.8 η
0.6 0.4 0.2 0 0
2
4
6
8
τ
Figure 1.3: Solutions o the dimensionless logistic equation. In Fig. 1.3, we plot the solution to the dimensionless logistic equation or initial conditions 0 = 0.02, 0.2, 0.5, 0.8, 1.0, and 1.2. The lowest curve is the characteristic ‘S-shape’ usually associated with the solution o the logistic equation. This sigmoidal curve appears in many other types o models. The MATLAB script to produce Fig. 1.3 is shown below. eta0=[0.02 .2 .5 .8 1 1.2]; tau=linspace(0,8); for i=1:length(eta0) eta=eta0(i)./(eta0(i)+(1-eta0(i)).*exp(-tau)); plot(tau,eta);hold on end axis([0 8 0 1.25]); xlabel(’\tau’); ylabel(’\eta’); title(’Logistic Equation’);
1.3
A model o species competition
Suppose that two species compete or the same resources. To build a model, we can start with logistic equations or both species. Diferent species would have diferent growth rates and diferent carrying capacities. I we let 1 and 2 be the number o individuals o species one and species two, then 1 = 1 1 (1 2 = 2 2 (1
− 1/ 1), − 2/ 2).
These are uncoupled equations so that asymptotically, 1 1 and 2 2 . How do we model the competition between species? I 1 is much smaller than 1 , and 2 much smaller than 2 , then resources are plentiul and populations
→
→
7
1.4. THE LOTKA-VOLTERRA PREDATOR-PREY MODEL
grow exponentially with growth rates 1 and 2 . I species one and two compete, then the growth o species one reduces resources available to species two, and vice-versa. Since we do not know the impact species one and two have on each other, we introduce two additional parameters to model the competition. A reasonable modication that couples the two logistic equations is
− −
1 = 1 1 1 2 = 2 2 1
1 + 12 2 1 21 1 + 2 2
,
(1.6a)
,
(1.6b)
where 12 and 21 are dimensionless parameters that model the consumption o species one’s resources by species two, and vice-versa. For example, suppose that both species eat exactly the same ood, but species two consumes twice as much as species one. Since one individual o species two consumes the equivalent o two individuals o species one, the correct model is 12 = 2 and 21 = 1/2. Another example supposes that species one and two occupy the same niche, consume resources at the same rate, but may have diferent growth rates and carrying capacities. Can the species coexist, or does one species eventually drive the other to extinction? It is possible to answer this question without actually solving the diferential equations. With 12 = 21 = 1 as appropriate or this example, the coupled logistic equations (1.6) become
−
1 = 1 1 1
1 + 2 1
,
−
2 = 2 2 1
1 + 2 2
.
(1.7)
For sake o argument, we assume that 1 > 2 . The only xed points other than the trivial one ( 1 , 2 ) = (0, 0) are ( 1 , 2 ) = ( 1 , 0) and ( 1 , 2 ) = (0, 2 ). Stability can be computed analytically by a two-dimensional Taylorseries expansion, but here a simpler argument can suce. We rst consider ( 1 , 2 ) = ( 1 , ), with small. Since 1 > 2 , observe rom (1.7) that ˙ 2 < 0 so that species two goes extinct. Thereore ( 1 , 2 ) = ( 1 , 0) is a stable xed point. Now consider ( 1 , 2 ) = (, 2 ), with small. Again, since ˙ 1 > 0 and species one increases in number. 1 > 2 , observe rom (1.7) that Thereore, ( 1 , 2 ) = (0, 2 ) is an unstable xed point. We have thus ound that, within our coupled logistic model, species that occupy the same niche and consume resources at the same rate cannot coexist and that the species with the largest carrying capacity will survive and drive the other species to extinction. This is the so-called principle o completive exclusion , also called -selection since the species with the largest carrying capacity wins. In act, ecologists also talk about -selection; that is, the species with the largest growth rate wins. Our coupled logistic model does not model -selection, demonstrating the potential limitations o a too simple mathematical model. For some values o 12 and 21 , our model admits a stable equilibrium solution where two species coexist. The calculation o the xed points and their stability is more complicated than the calculation just done, and I present only the results. The stable coexistence o two species within our model is possible only i 12 2 < 1 and 21 1 < 2 .
8
CHAPTER 1. POPULATION DYNAMICS
Figure 1.4: Pelt-trading records o the Hudson Bay Company or the snowshoe hare and its predator the lynx. [From E.P. Odum, Fundamentals o Ecology, 1953.]
1.4
The Lotka-Volterra predator-prey model
Pelt-trading records (Fig. 1.4) o the Hudson Bay company rom over almost a century display a near-periodic oscillation in the number o trapped snowshoe hares and lynxes. With the reasonable assumption that the recorded number o trapped animals is proportional to the animal population, these records suggest that predator-prey populations—as typied by the hare and the lynx—can oscillate over time. Lotka and Volterra independently proposed in the 1920s a mathematical model or the population dynamics o a predator and prey, and these Lotka-Volterra predator-prey equations have since become an iconic model o mathematical biology. To develop these equations, suppose that a predator population eeds on a prey population. We assume that the number o prey grow exponentially in the absence o predators (there is unlimited ood available to the prey), and that the number o predators decay exponentially in the absence o prey (predators must eat prey or starve). Contact between predators and prey increases the number o predators and decreases the number o prey. Let () and () be the number o prey and predators at time . To develop a coupled diferential equation model, we consider population sizes at time +Δ. Exponential growth o prey in the absence o predators and exponential decay o predators in the absence o prey can be modeled by the usual linear terms. The coupling between prey and predator must be modeled with two additional parameters. We write the population sizes at time + Δ as ( + Δ) = () + Δ () Δ () (), ( + Δ) = () + Δ () () Δ ().
−
−
The parameters and are the average per capita birthrate o the prey and the deathrate o the predators, in the absence o the other species. The coupling
9
1.4. THE LOTKA-VOLTERRA PREDATOR-PREY MODEL
terms model contact between predators and prey. The parameter is the raction o prey caught per predator per unit time; the total number o prey caught by predators during time Δ is Δ . The prey eaten is then converted into newborn predators (view this as a conversion o biomass), with conversion actor , so that the number o predators during time Δ increases by Δ . Converting these equations into diferential equations by letting Δ 0, we obtain the well-known Lotka-Volterra predator-prey equations
→
=
− ,
=
− .
(1.8)
Beore analyzing the Lotka-Volterra equations, we rst review xed point and linear stability analysis applied to what is called an autonomous system o diferential equations. For simplicity, we consider a system o only two diferential equations o the orm ˙ = (, ),
˙ = (, ),
(1.9)
though our results can be generalized to larger systems. The system given by (1.9) is said to be autonomous since and do not depend explicitly on the independent variable . Fixed points o this system are determined by setting ˙ = ˙ = 0 and solving or and . Suppose that one xed point is (* , * ). To determine its linear stability, we consider initial conditions or (, ) near the xed point with small independent perturbations in both directions, i.e., (0) = * + (0), (0) = * + (0). I the initial perturbation grows in time, we say that the xed point is unstable; i it decays, we say that the xed point is stable. Accordingly, we let () = * + (),
() = * + (),
(1.10)
and substitute (1.10) into (1.9) to determine the time-dependence o and . Since * and * are constants, we have ˙ = (* + , * + ),
˙ = (* + , * + ).
The linear stability analysis proceeds by assuming that the initial perturbations (0) and (0) are small enough to truncate the two-dimensional Taylor-series expansion o and about = = 0 to rst-order in and . Note that in general, the two-dimensional Taylor series o a unction (, ) about the origin is given by (, ) = (0, 0) + (0, 0) + (0, 0) 1 2 + (0, 0) + 2 (0, 0) + 2 (0, 0) + . . . , 2
[
]
where the terms in the expansion can be remembered by requiring that all o the partial derivatives o the series agree with that o (, ) at the origin. We now Taylor-series expand (* + , * + ) and (* + , * + ) about (, ) = (0, 0). The constant terms vanish since (* , * ) is a xed point, and we neglect all terms with higher orders than and . Thereore, ˙ = (* , * ) + (* , * ),
˙ = (* , * ) + (* , * ),
10
CHAPTER 1. POPULATION DYNAMICS
which may be written in matrix orm as
=
* *
* *
,
(1.11)
where * = (* , * ), etc. Equation (1.11) is a system o linear ode’s, and its solution proceeds by assuming the orm
= v.
(1.12)
Upon substitution o (1.12) into (1.11), and canceling , we obtain the linear algebra eigenvalue problem *
*
J v = v, with J =
* *
* , *
where is the eigenvalue, v the corresponding eigenvector, and J* the Jacobian matrix evaluated at the xed point. The eigenvalue is determined rom the characteristic equation det(J* I) = 0,
−
which or a two-by-two Jacobian matrix results in a quadratic equation or . From the orm o the solution (1.12), the xed point is stable i or all eigenvalues , Re < 0, and unstable i or at least one , Re > 0. Here Re means the real part o the (possibly) complex eigenvalue . We now reconsider the Lotka-Volterra equations (1.8). Fixed point solutions ˙ = ˙ = 0, and the two solutions are are ound by solving
{}
{}
( * , * ) = (0, 0) or (
{}
, ).
The trivial xed point (0, 0) is unstable since the prey population grows exponentially i it is initially small. To determine the stability o the second xed point, we write the Lotka-Volterra equation in the orm = (, ),
= (, ),
with (, ) =
− ,
(, ) =
The partial derivatives are then computed to be = , = ,
− .
= = .
−
−
The Jacobian at the xed point ( * , * ) = (/,/ ) is *
J =
0
−/
⃒⃒− ⃒
−/ −
0
;
and det(J*
−
I) =
= 2 + =0
⃒⃒ ⃒
−
1.4. THE LOTKA-VOLTERRA PREDATOR-PREY MODEL
11
± √
has the solutions ± = , which are pure imaginary. When the eigenvalues o the two-by-two Jacobian are pure imaginary, the xed point is called a center and the perturbation neither grows nor decays, but oscillates. Here, the angular requency o oscillation is = , and the period o the oscillation is 2/. We plot and versus (time series plot), and versus (phase space diagram) to see how the solutions behave. For a nonlinear system o equations such as (1.8), a numerical solution is required. The Lotka-Volterra system has our ree parameters , , and . The relevant units here are time, the number o prey, and the number o predators. The Buckingham Pi Theorem predicts that nondimensionalizing the equations can reduce the number o ree parameters by three to a manageable single dimensionless grouping o parameters. We choose to nondimensionalize time using the angular requency o oscillation and the number o prey and predators using their xed point values. With carets denoting the dimensionless variables, we let
√
^ =
√
,
^ = / * = ,
^ = / * = .
(1.13)
Substitution o (1.13) into the Lotka-Volterra equations (1.8) results in the dimensionless equations ^ ^ = ( ^
^ − ^ ),
^ 1 ^^ = ( ^
√
^ − ),
with single dimensionless grouping = / . A numerical solution uses MATLAB’s ode45.m built-in unction to integrate the diferential equations. The code below produces Fig. 1.5. Notice how the predator population lags the prey population: an increase in prey numbers results in a delayed increase in predator numbers as the predators eat more prey. The phase space diagrams clearly show the periodicity o the oscillation. Note that the curves move counterclockwise: prey numbers increase when predator numbers are minimal, and prey numbers decrease when predator numbers are maximal.
12
CHAPTER 1. POPULATION DYNAMICS
function lotka_volterra % plots time series and phase space diagrams clear all; close all; t0=0; tf=6*pi; eps=0.1; delta=0; r=[1/2, 1, 2]; options = odeset(’RelTol’,1e-6,’AbsTol’,[1e-6 1e-6]); %time series plots for i=1:length(r); [t,UV]=ode45(@lv_eq,[t0,tf],[1+eps 1+delta],options,r(i)); U=UV(:,1); V=UV(:,2); subplot(3,1,i); plot(t,U,t,V,’ --’); axis([0 6*pi,0.8 1.25]); ylabel(’predator,prey’); text(3,1.15,[’r=’,num2str(r(i))]); end xlabel(’t’); subplot(3,1,1); legend(’prey’, ’predator’); %phase space plot xpos=[2.5 2.5 2.5]; ypos=[3.5 3.5 3.5];%for annotating graph for i=1:length(r); for eps=0.1:0.1:1.0; [t,UV]=ode45(@lv_eq,[t0,tf],[1+eps 1+eps],options,r(i)); U=UV(:,1); V=UV(:,2); figure(2);subplot(1,3,i); plot(U,V); hold on; end axis equal; axis([0 4 0 4]); text(xpos(i),ypos(i),[’r=’,num2str(r(i))]); if i==1; ylabel(’predator’); end; xlabel(’prey’); end function dUV=lv_eq(t,UV,r) dUV=zeros(2,1); dUV(1) = r*(UV(1)-UV(1)*UV(2)); dUV(2) = (1/r)*(UV(1)*UV(2)-UV(2));
13
1.4. THE LOTKA-VOLTERRA PREDATOR-PREY MODEL
y e r p , r o t a d e r p
y e r p , r
o t a d e r p
y e r p , r o t a d e r p
prey predator
1.2
r=0.5
1 0.8
0
2
1.2
4
6
8
10
12
14
16
18
6
8
10
12
14
16
18
6
8
10
12
14
16
18
r=1
1 0.8
0
2
1.2
4 r=2
1 0.8
0
2
4
t
4
4
4
r=0.5 r o t a d e r p
r=1
r=2
3
3
3
2
2
2
1
1
1
0
0
2 prey
4
0
0
2 prey
4
0
0
2 prey
4
Figure 1.5: Solutions o the dimensionless Lotka-Volterra equations. Upper plots: time-series solutions; lower plots: phase space diagrams.
14
CHAPTER 1. POPULATION DYNAMICS
Chapter 2
Age-structured Populations Determining the age-structure o a population helps governments plan economic development. Age-structure theory can also help evolutionary biologists better understand a species’s lie-history. An age-structured population occurs because ofspring are born to mothers at diferent ages. I average per capita birth and death rates at diferent ages are constant, then a stable age-structure arises. However, a rapid change in birth or death rates can cause the age-structure to shit distributions. In this section, we develop the theory o age-structured populations using both discrete- and continuous-time models. We also present two interesting applications: (1) modeling age-structure changes in China and other countries as these populations age, and; (2) modeling the lie cycle o a hermaphroditic worm. We begin this section, however, with one o the oldest problems in mathematical biology: Fibonacci’s rabbits. This will lead us to a brie digression about the golden mean, rational approximations and ower development, beore returning to our main topic.
2.1
Fibonacci’s rabbits
In 1202, Fibonacci proposed the ollowing puzzle, which we paraphrase here: A man put a male-emale pair o newly born rabbits in a eld. Rabbits take a month to mature beore mating. One month ater mating, emales give birth to one male-emale pair and then mate again. Assume no rabbits die. How many rabbit pairs are there ater one year? Let be the number o rabbit pairs at the start o the th month, counted ater the birth o any new rabbits. There are twelve months in a year, so the number o rabbits at the start o the 13th month will be the solution to Fibonacci’s puzzle. Now 1 = 1 because the population consists o the initial newborn rabbit pair at the start o the rst month; 2 = 1 because the population still consists o the initial rabbit pair, now mature; and 3 = 2 + 1 because the population consists o 2 rabbit pairs present in the previous month (here equal to one), and 1 newborn rabbit pairs born to the 1 emale rabbits that are now old enough to give birth (here, also equal to one). In general +1 = + −1 . 15
(2.1)
16
CHAPTER 2. AGE-STRUCTURED POPULATIONS
For only 12 months we can simply count using (2.1) and 1 = 2 = 1, so the Fibonacci numbers are 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, . . . where 13 = 233 is the solution to Fibonacci’s puzzle. Let us solve (2.1) or all the ’s. This equation is a second-order linear diference equation, and to solve it, we look or a solution o the orm = . Substitution into (2.1) yields +1 = + −1 , or ater division by −1 : 2 with solution
− − 1 = 0, √ 1± 5
± = Dene
.
2
√
1+ 5 Φ= = 1.61803 . . . , 2 and =
√ 5 − 1 2
=Φ
− 1 = 0.61803 . . . .
Then + = Φ and − = . Also, notice that since Φ 2 by Φ yields 1/Φ = Φ 1, so that
−
−
=
− Φ − 1 = 0, division
1 . Φ
As in the solution o linear homogeneous diferential equations, the two values o can be used to construct a general solution to the linear diference equation using the principle o linear superposition: = 1 Φ + 2 ( 1) .
−
Extending the Fibonacci sequence to 0 = 0 (since 0 = 2 conditions 0 = 0 and 1 = 1:
− 1), we satisy the
1 + 2 = 0, 1 Φ 2 = 1.
−
√
−1/√ 5. We can
1 + ( 1)+1 Φ−2 .
(2.2)
Thereore, 2 = 1 , and 1 (Φ + ) = 1, or 1 = 1/ 5, 2 = rewrite the solution as
−
=
√ 15 Φ
In this orm, we see that, as
[
−
]
→ ∞, → Φ/√ 5, and +1/ → Φ.
2.2. THE GOLDEN RATIO Φ
2.2
The golden ratio
17 Φ
The number Φ is known as the golden ratio. Two positive numbers and , with > , are said to be in the golden ratio i the ratio between the sum o those numbers and the larger one is the same as the ratio between the larger one and the smaller; that is, + = . (2.3) Solution o (2.3) yields / = Φ. In some well-dened way, Φ can also be called the most irrational o the irrational numbers. To understand why Φ has this distinction as the most irrational number, we need rst to understand an algorithm—continued raction expansions—or approximating irrational numbers by rational numbers. Let be positive integers. To construct a continued raction expansion to the positive irrational number , we choose the largest possible ’s that satisy the ollowing inequalities: > 1 = 1 , < 2 = 1 + > 3 = 1 +
< 4 = 1 +
1 , 2 1 1 2 + 3 1 2 +
,
,
1 3 +
etc.
1 4
There is a simple algorithm or calculating the ’s. Firstly, 1 = 1 is the integer part o . One computes the remainder, 1 , and takes its reciprocal 1/( 1 ). The integer part o this reciprocal is 2 . One then takes the remainder 1/( 1 ) 2 , and its reciprocal. The integer part o this reciprocal is 3 . We ollow this algorithm to nd successively better rational approximations o = 3.141592654 . . . :
− −
−
−
’s
remainder 1/remainder
1 = 3 2 = 3 + 3 = 3 +
4 = 3 +
0.14159265 7.06251329 1 7
0.06251329 15.99659848 1
0.99659848 1.00341313
1 7+ 15 1 7+
1 15 +
1 1
18
CHAPTER 2. AGE-STRUCTURED POPULATIONS
Follow how this works. The rst rational approximation o is 1 = 3, which is less than . The remainder is 0.141592654. The remainder’s reciprocal is 7.06251329. The integer part o this reciprocal is 7, so our next rational approximation o is 2 = 3 + 1/7, which is larger than since 1/7 is greater than 0.141592654. The next remainder is 7.06251329 7 = 0.06251329, and its reciprocal is 15.99659848. Thereore, our next rational approximation is 3 = 3 + 1/(7 + (1/15)), and since 1/15 is greater than 0.06251329, 3 is less than . Rationalizing the continued ractions gives us the ollowing successive rational number approximation o : 3, 22/7, 333/106, 355/113, . . . , or in decimal orm: 3, 3.1428 . . . , 3.141509 . . . , 3.14159292 . . . . Now consider a rational approximation or Φ = 1.61803399 . . . :
−
{
{
’s
}
remainder
1/remainder
Φ
Φ
1 = 1 2 = 1 +
}
1 1
and so on. All the ’s are one, and the sequence results in the slowest possible convergence or the ’s. Formerly, the limiting expression or Φ is 1
Φ= 1+
.
1 1+
1 ...
Another derivation o the continued raction expansion or Φ can be obtained rom Φ2 Φ 1 = 0, written as Φ = 1 + 1/Φ. Then inserting this expression or Φ back into the right-hand-side gives Φ = 1 + 1/(1 + 1/Φ), and iterating innitum shows that all the ’s are one. (As a side note, another interesting expression or Φ comes rom iterating Φ = 1 + Φ.) The sequence o rational approximations just constructed or Φ is 1, 2, 3/2, 5/3, 8/5, 13/8, etc., which are just the ratios o consecutive Fibonacci numbers. We have already shown that this ratio converges to Φ. Because the golden ratio is the most irrational number, it has a way o appearing unexpectedly in nature. One well-known example is the orets in a sunower head, which we discuss in the next section.
− −
√
2.3
The Fibonacci numbers in a sunower
Consider the photo o a sunower shown in Fig. 2.1, and notice the apparent spirals in the orets radiating out rom the center to the edge. These spirals appear to rotate both clockwise and counterclockwise. By counting them, one nds 34 clockwise spirals and 21 counterclockwise spirals. The numbers 21 and 34 are notable as they are consecutive numbers in the Fibonacci sequence. Why do Fibonacci numbers appear in the sunower head? To answer this question, we construct a very simple model or the way that the orets develop. Suppose that during development, orets rst appear close to the center o the head and subsequently move radially out with constant speed as the sunower head grows. To ll the circular sunower head, we suppose that as each new oret is created at the center, it is rotated through a constant angle beore
2.3. THE FIBONACCI NUMBERS IN A SUNFLOWER
19
Figure 2.1: The owering head o a sunower. moving radially. We will urther assume that the angle o rotation is optimum in the sense that the resulting sunower head consists o orets that are wellspaced. We denote the rotation angle by 2, and we can assume that 0 < < 1. We rst consider the possibility that is a rational number, say /, where and are positive integers with no common actors, and < . Since ater rotations the orets return to the line on which they started, the resulting sunower head would consist o orets lying along straight lines. Such a sunower head or = 3/5 is shown in Fig. 2.2a, where ve straight lines are apparent. This gure and subsequent ones were produced using MATLAB, and the simulation code is shown at the end o this subsection. We next consider the possibility that is an irrational number. Now, no matter how many rotations, the orets will never return to their starting line. Nevertheless, the resulting sunower head may not have orets that are wellspaced. For example, i = 3, then the resulting sunower head would look like Fig. 2.2b. Intriguingly, seven counterclockwise spirals are evident. We have already shown that a good rational approximation or —obtained rom a continued raction expansion—is 3 + 1/7, which is slightly larger than . So that on every seventh rotation, the new oret alls just short o the radial line traced by the oret rom seven rotations ago. Ater every seven rotations as these orets move out along the growing sunower head, they will take on the appearance o a counterclockwise spiral; and the entire ower head will seem to consist o seven o these spirals, as evident in Fig. 2.2b. The irrational number that is most likely to result in a sunower head with well-spaced orets is the most irrational o the irrational numbers, that is , which is the ractional part o the golden ratio. The rational approximations
−
20
CHAPTER 2. AGE-STRUCTURED POPULATIONS
(a)
(b)
Figure 2.2: Simulation o a sunower head or (a) = 3/5; (b) =
(a)
− 3.
(b)
Figure 2.3: Simulation o a sunower head or = . (a) 0 = 1/2; (b) 0 = 1/4.
2.4. RABBITS ARE AN AGE-STRUCTURED POPULATION
21
to obtained rom a continued raction expansion are just the ratio o two consecutive Fibonacci numbers, e.g., 13/21 or 21/34. Depending on the density and size o the sunower head, diferent consecutive Fibonacci numbers may be observed. For the sunower shown in Fig. 2.1, the relevant Fibonacci numbers are 13, 21 and 34. Since 13/21 > , the 21 spirals are expected to rotate counterclockwise; and since 21/34 < , the 34 spirals are expected to rotate clockwise. Furthermore, the 34 clockwise spirals should be easier to observe than the 21 counterclockwise spirals (since 21/34 is a better approximation to than is 13/21), and indeed, this is exactly what we nd. Two simulations o the sunower head with = are shown in Fig. 2.3. These simulations difer only by the choice o radial velocity, with 0 = 1/2 in Fig. 2.3a and 0 = 1/4 in Fig. 2.3b, where the precise denition o 0 can be obtained rom the MATLAB code. In Fig. 2.3a, there can be counted 13 counterclockwise spirals and 21 clockwise spirals; in Fig. 2.3b, there are 21 clockwise spirals and 34 counterclockwise spirals, similar to the real sunower in Fig. 2.1. As expected, these number o spirals and their orientation derive rom the ratio o consecutive Fibonacci numbers satisying 8/13 < , 13/21 > and 21/34 < . function flower(alpha) %draws points moving radially but rotated an angle 2*pi*num figure angle=2*pi*alpha; r0=0; v0=1.0; %for higher Fibonacci numbers try v0=0.5, 0.25, 0.1 npts=100/v0; theta=(0:npts-1)*angle; %all the angles of the npts r=r0*ones(1,npts); %starting radius of the npts for i=1:npts % i is total number of points to be plotted (also time) for j=1:i %find coordinates of all i points; j=i is newest point r(j)=r(j)+v0; %*(i-j); x(j)=r(j)*cos(theta(j)); y(j)=r(j)*sin(theta(j)); end plot(x,y,’.’,’MarkerSize’,5.0); axis equal; axis([-npts*v0-r0 npts*v0+r0 -npts*v0-r0 npts*v0+r0]); pause(0.1); end %draw a circle at the center hold on; phi=linspace(0,2*pi); x=r0*cos(phi); y=r0*sin(phi); fill(x,y,’r’)
2.4
Rabbits are an age-structured population
Fibonacci’s rabbits orm an age-structured population and we can use this simple case to illustrate the more general approach. Fibonacci’s rabbits can be categorized into two meaningul age classes: juveniles and adults. Here, juveniles are the newborn rabbits that can not yet mate; adults are those rabbits at least one month old. Beginning with a newborn pair at the beginning o the rst month, we census the population at the beginning o each subsequent month
22
CHAPTER 2. AGE-STRUCTURED POPULATIONS
ater mated emales have given birth. At the start o the th month, let 1, be the number o newborn rabbit pairs, and let 2, be the number o rabbit pairs at least one month old. Since each adult pair gives birth to a juvenile pair, the number o juvenile pairs at the start o the ( + 1)-st month is equal to the number o adult pairs at the start o the -th month. And since the number o adult pairs at the start o the ( + 1)-st month is equal to the sum o adult and juvenile pairs at the start o the -th month, we have 1,+1 = 2, , 2,+1 = 1, + 2, ; or written in matrix orm
1,+1 2,+1
0 1 1 1
=
1, . 2,
(2.4)
Rewritten in vector orm, we have u+1 = Lu ,
(2.5)
where the denitions o the vector u and the matrix L are obvious. The initial conditions, with one juvenile pair and no adults, are given by
1,1 2,1
=
1 . 0
The solution o this system o coupled, rst-order, linear, diference equations, (2.5), proceeds similarly to that o coupled, rst-order, linear, diferential equations. With the ansatz, u = v, we obtain upon substitution into (2.5) the eigenvalue problem Lv = v, whose solution yields two eigenvalues 1 and 2 , with corresponding eigenvectors v1 and v2 . The general solution to (2.5) is then u = 1 1 v1 + 2 2 v2 ,
(2.6)
with 1 and 2 determined rom the initial conditions. Now suppose that 1 > 2 . I we rewrite (2.6) in the orm
| |
| |
u =
1
1 v1 + 2
2 1
v2 ,
then because 2 /1 < 1, u 1 1 v1 as . The long-time asymptotics o the population, thereore, depends only on 1 and the corresponding eigenvector v1 . For our Fibonacci’s rabbits, the eigenvalues are obtained by solving det(L I) = 0, and we nd
|
|
→
→∞
−
det
−
1 1 1
−
=
−(1 − ) − 1
= 0,
23
2.5. DISCRETE AGE-STRUCTURED POPULATIONS
, = 1 = 0 =
···
ℛ
number o emales in age class at census raction o emales surviving rom age class 1 to expected number o emale ofspring rom a emale in age class raction o emales surviving rom birth to age class
−
∑
basic reproductive ratio
Table 2.1: Denitions needed in an age-structured, discrete-time population model or 2 1 = 0, with solutions Φ and . Since Φ > , the eigenvalue Φ and its corresponding eigenvector determine the long-time asymptotic population age-structure. The eigenvector may be ound by solving
− −
−
(L or
−
− ΦI)v1 = 0,
Φ 1 1 1 Φ
−
11 12
=
0 . 0
The rst equation is just Φ times the second equation (use Φ2 so that 12 = Φ11 . Taking 11 = 1, we have
−
v1 =
− Φ − 1 = 0),
1 . Φ
The asymptotic age-structure obtained rom v1 shows that the ratio o adults to juveniles approaches the golden mean; that is, 2, = 12 /11 →∞ 1, lim
= Φ.
2.5
Discrete age-structured populations
In a discrete model, population censuses occur at discrete times and individuals are assigned to age classes, spanning a range o ages. For model simplicity, we assume that the time between censuses is equal to the age span o all age classes. An example is a country that censuses its population every ve years, and assigns individuals to age classes spanning ve years (e.g., 0-4 years old, 5-9 years old, etc.). Although country censuses commonly count both emales and males separately, we will only count emales and ignore males. There are several new denitions in this Section and I place these in Table 2.1 or easy reerence. We dene , to be the number o emales in age class at census . We assume that = 1 represents the rst age class and = the last. No emale survives past the last age class. We also assume that the rst census takes place when = 1. We dene as the raction o emales that survive rom age class 1 to age class (with 1 the raction o newborns that survive to their rst census), and dene as the expected number o emale births per emale in age class .
−
24
CHAPTER 2. AGE-STRUCTURED POPULATIONS
We construct diference equations or ,+1 in terms o , . First, newborns at census + 1 were born between census and + 1 to diferent aged emales, with difering ertilities. Also, only a action o these newborns survive to their rst census. Second, only a raction o emales in age class that were counted in census survive to be counted in age class + 1 in census + 1. Putting these two ideas together with the appropriately dened parameters, the diference equations or ,+1 are determined to be
{
{
}
{ }
}
1,+1 = 1 (1 1, + 2 2, + 2,+1 = 2 1, , 3,+1 = 3 2, , .. . ,+1 = −1, ,
··· + , ) ,
which can be rewritten as the matrix equation
⎛ ⎜⎜ ⎜⎝
1,+1 2,+1 3,+1 .. .
, +1
⎞ ⎟⎟ ⎟⎠
=
⎛ ⎜⎜ ⎜⎝
1 1 2 0 .. .
1 2 0 3 .. .
. .. ... ... .. .
1 −1 0 0 .. .
1 0 0 .. .
0
0
...
0
or in compact vector orm as
⎞⎛ ⎟⎟ ⎜⎜ ⎟⎠ ⎜⎝
1, 2, 3, .. .
,
⎞ ⎟⎟ ⎟⎠
u+1 = Lu ,
;
(2.7)
where L is called the Leslie Matrix. This system o linear equations can be solved by determining the eigenvalues and associated eigenvectors o the Leslie Matrix. One can solve directly the characteristic equation, det (L I) = 0, or reduce the system o rst-order diference equations (2.7) to a single high-order equation or the number o emales in the rst age class. Following the latter approach, and beginning with the second row o (2.7), we have
−
2,+1 = 2 1, , 3,+1 = 3 2, = 3 2 1,−1 , .. . ,+1 = −1, = −1 −2,−1 .. . = −1
··· 21,
− +2
.
I we dene = 1 2 to be the raction o emales that survive rom birth to age class , and = , then the rst row o (2.7) becomes
···
1,+1 = 1 1, + 2 1,−1 + 3 1,−2 +
··· + 1,
− +1
.
(2.8)
25
2.5. DISCRETE AGE-STRUCTURED POPULATIONS
Here, we have made the simpliying assumption that so that all the emales counted in the + 1 census were born ater the rst census. The high-order linear diference equation (2.8) may be solved using the ansatz 1, = . Direct substitution and division by +1 results in the discrete Euler-Lotka equation
≥
− = 1,
(2.9)
=1
which may have both real and complex-conjugate roots. Once an eigenvalue is determined rom (2.9), the corresponding eigenvector v can be computed using the Leslie matrix. We have
⎛ ⎜⎜ ⎜⎝
1 1 2 0 .. . 0
−
1 2 3 .. .
−
0
. . . 1 −1 ... 0 ... 0 .. .. . . ...
1 0 0 .. .
−
⎞⎛ ⎟⎟ ⎜⎜ ⎟⎠ ⎜⎝
1 2 3 .. .
⎞ ⎟⎟ ⎟⎠
=
⎛⎞ ⎜⎜ ⎟⎟ ⎜⎝ ⎟⎠ 0 0 0 .. .
.
0
Taking = / , and beginning with the last row and working backwards, we have: −1 = −1 /−1 , −2 = −2 /−2 , .. . 1 = 1 /, so that
= / ,
or = 1, 2, . . . , .
We can obtain an interesting implication o this result by orming the ratio o two consecutive age classes. I is the dominant eigenvalue (and is real and positive, as is the case or human populations), then asymptotically, +1, /,
∼ +1/
= +1 /. With the survival ractions xed, a smaller positive implies a larger ratio: a slower growing (or decreasing) population has relatively more older people than a aster growing population. In act, we are now living through a time when developed countries, particularly Japan and those in Western Europe, as well as Hong Kong and Singapore, have substantially lowered their population growth rates and are increasing the average age o their citizens. I we want to simply determine i a population grows or decays, we can calculate the basic reproduction ratio 0 , dened as the net expectation o emale ofspring to a newborn emale. Stasis is obtained i the emale only replaces hersel beore dying. I 0 > 1, then the population grows, and i 0 < 1 then the population decays. 0 is equal to the number o emale ofspring expected rom a newborn when she is in age class , summed over all age classes, or
{ }
ℛ
ℛ
ℛ
ℛ
ℛ0 =
=1
.
26
CHAPTER 2. AGE-STRUCTURED POPULATIONS
For a population with approximately equal numbers o males and emales, 0 = 1 means a newborn emale must produce on average two children over her lietime. News stories in the western press requently state that or zero population growth, women need to have 2.1 children. The term women used in these stories presumably means women o child-bearing age. Since girls who die young have no children, the statistic o 2.1 children implies that 0.1/2.1, or about 5% o children die beore reaching adulthood. A useul application o the mathematical model developed in this Section is to predict the uture age structure within various countries. This can be important or economic planning—or instance, determining the tax revenues that can pay or the rising costs o health care as a population ages. For accurate predictions on the uture age-structure o a given country, immigration and migration must also be modeled. An interesting website to browse is at
ℛ
http://www.census.gov/ipc/www/idb. This website, created by the US census bureau, provides access to the International Data Base (IDB), a computerized source o demographic and socioeconomic statistics or 227 countries and areas o the world. In class, we will look at and discuss the dynamic output o some o the population pyramids, including those or Hong Kong and China.
2.6
Continuous age-structured populations
We can derive a continuous-time model by considering the discrete model in the limit as the age span Δ o an age class (also equal to the time between censuses) goes to zero. For > , (2.8) can be rewritten as
1, =
1,− .
(2.10)
=1
The rst age class in the discrete model consists o emales born between two consecutive censuses. The corresponding unction in the continuous model is the emale birth rate o the population as a whole, (), satisying 1, = ( )Δ. I we assume that the th census takes place at a time = Δ, we also have 1,− = (− )Δ = ( )Δ.
−
To determine the continuous analogue o the parameter = , we dene the age-specic survival unction () to be the raction o newborn emales that survive to age , and dene the age-specic maternity unction (), multiplied by Δ, to be the average number o emales born to a emale between the ages o and + Δ. With the denition o the age-specic net maternity unction, () = ()(), and = Δ, we have = ( )Δ.
2.6. CONTINUOUS AGE-STRUCTURED POPULATIONS
27
With these new denitions, (2.10) becomes
( )Δ =
( )(
=1
− )(Δ)2.
Canceling one actor o Δ, and using = , the right-hand side becomes a Riemann sum. Taking = and assigning () = 0 when is greater than the maximum age o emale ertility, the limit Δ 0 transorms (2.10) to
→
∞
() =
(
0
− ) ().
(2.11)
Equation (2.11) states that the population-wide emale birth rate at time has contributions rom emales o all ages, and that the contribution to this birth rate rom emales between the ages o and + is determined rom the population-wide emale birth rate at the earlier time times the raction o emales that survive to age times the number o emale births to emales between the ages o and + . Equation (2.11) is a linear homogeneous integral equation, valid or greater than the maximum age o emale ertility. A more complete but inhomogeneous equation valid or smaller can also be derived. Equation (2.11) can be solved by the ansatz () = . Direct substitution yields
−
∞
=
()(−) ,
0
which upon canceling results in the continuous Euler-Lotka equation ∞
()− = 1.
(2.12)
0
Equation (2.12) is an integral equation or given the age-specic net maternity unction (). It is possible to prove that or () a continuous non-negative unction, (2.13) has exactly one real root * , and that the population grows (* > 0) or decays (* < 0) asymptotically as . The population growth rate * has been called the intrinsic rate o increase, the intrinsic growth rate, or the Malthusian parameter. Typically, (2.12) must be solved or * using Newton’s method. We briey digress here to explain Newton’s method, which is an ecient root-nding algorithm that solves () = 0 or . Newton’s method can be derived graphically by plotting () versus , approximating () by the tangential line at = with slope ′ ( ), and letting +1 be the intercept o the tangential line with the x-axis. An alternative derivation starts with the Taylor series (+1 ) = ( ) + (+1 ) ′ ( ) + . . . . We set (+1 ) = 0, drop higher-order terms in the Taylor series, and solve or +1 : *
−
+1 =
− (()) , ′
which is to be solved iteratively, starting with a wise choice o initial guess 0 , until convergence.
28
CHAPTER 2. AGE-STRUCTURED POPULATIONS
Solution o (2.12) by Newton’s method requires dening ∞
() =
()−
0
− 1,
(2.13)
rom which one derives ∞
′
() =
−
()− .
0
Newton’s method then begins with an initial guess 0 and computes the iteration ∞
+1 = +
()− 1 ∞ ()− 0
∫ ∫ 0
−
until converges suciently close to * . Ater asymptotically attaining a stable age structure, the population grows like , and our previous discussion o the Malthusian growth model suggests that * may be ound rom the constant per capita birth rate and death rate . By determining expressions or and , we will indeed show that * = . Because emales that have survived to age at time were born earlier at time , ( )() represents the number o emales at time that are between the ages o and + . The total number o emales () at time is thereore given by *
−
−
−
∞
() =
(
0
− )().
(2.14)
The per capita birth rate () equals the population-wide birth rate () divided by the population size (), and using (2.11) and (2.14), () = ()/ () ∞
=
∫ ∫ 0
∞
0
( (
− ) () . − )()
Similarly, the per capita death rate () equals the population-wide death rate () divided by (). To derive (), we rst dene the age-specic mortality unction (), multiplied by Δ, to be the raction o emales o age that die beore attaining the age + Δ. The relationship between the agespecic mortality unction () and the age-specic survival unction () may be obtained by computing the raction o emales that survive to age + Δ. This raction is equal to the raction o emales that survive to age times the raction o emales 1 ()Δ that do not die in the next small interval o time Δ; that is, ( + Δ) = ()(1 ()Δ),
−
−
or
and as Δ
( + Δ) Δ
→ 0,
− () = −()();
′ () =
−()().
(2.15)
The age-specic mortality unction () is analogous to the age-specic maternity unction (), and we dene the age-specic net mortality unction () =
29
2.7. THE BROOD SIZE OF A HERMAPHRODITIC WORM
()() in analogy to the age-specic net maternity unction () = ()(). The population-wide birth rate () is determined rom () using (2.11), and in analogy, the population-wide death rate () is determined rom () using ∞
() =
(
0
− )(),
(2.16)
where the integrand represents the contribution to the death rate rom emales that die between the ages o and + . The per capita death rate is thereore () = ()/ () ∞
∫ ∫ − − ∫ ∫ − 0
=
∞
0
( (
)() ; )()
and the diference between the per capita birth and death rates is calculated rom ∞ ( ) [ () ()] 0 () () = . (2.17) ∞ ( )() 0
−
−
−
Asymptotically, a stable age structure is established and the population-wide birth rate grows as () . Substitution o this expression or () into (2.17) and cancelation o results in
∼
*
*
∞
−= =
[ () ()] − ∞ ()− 0
∫ ∫ ∫ ∫ 0
−
*
*
∞ ′
()− , ∞ ()− 0
1+
*
0
*
where use has been made o (2.12) and (2.15). Simpliying the numerator using integration by parts, ∞
∞
′
−*
()
−* ∞
= ()
|0
0
+ *
∞
=
−1 +
*
()− *
0
()− , *
0
produces the desired result, * =
− .
It is usually supposed that evolution by natural selection will result in populations with the largest value o the Malthusian parameter * , and that natural selection would avor those emales that constitute such a population. We will exploit this idea in the next section to compute the brood size o the selertilizing hermaphroditic worm o the species Caenorhabditis elegans .
2.7
The brood size o a hermaphroditic worm
Caenorhabditis elegans , a soil-dwelling nematode worm about 1 mm in length, is a widely studied model organism in biology. With a body made up o approximately 1000 cells, it is one o the simplest multicellular organisms under study. Advances in understanding the development o this multicellular organism led
30
CHAPTER 2. AGE-STRUCTURED POPULATIONS
Figure 2.4: Caenorhabditis elegans , a nematode worm used by biologists as a simple animal model o a multicellular organism. (Photo by Amy Pasquinelli.)
sperm-
0
juvenile
-
-
eggs
+
++
adult
-
Figure 2.5: A simplied timeline o a hermaphrodite’s lie. to the awarding o the 2002 Nobel prize in Physiology or Medicine to the three C. elegans biologists Sydney Brenner, H. Robert Horvitz and John E. Sulston. The worm C. elegans has two sexes: hermaphrodites, which are essentially emales that can produce internal sperm and sel-ertilize their own eggs, and males, which must mate with hermaphrodites to produce ofspring. In laboratory cultures, males are rare and worms generally propagate by sel-ertilization. Typically, a hermaphrodite lays about 250-350 sel-ertilized eggs beore becoming inertile. It is reasonable to assume that the orces o natural selection have shaped the lie-history o C. elegans , and that the number o ofspring produced by a selng hermaphrodite must be in some sense optimal. Here, we show how an age-structured model applied to C. elegans yields theoretical insights into the brood size o a selng hermaphrodite. To develop a mathematical model or C. elegans , we need to know some details o its lie history. As a rst approximation (Barker, 1992), a simplied timeline o a hermaphrodite’s lie is shown in Fig. 2.5. The ertilized egg is laid at time = 0. During a juvenile growth period, the immature worm develops through our larval stages (L1-L4). Towards the end o L4 and or a short while ater its nal molt to adulthood, the hermaphrodite produces sperm, which is then stored or later use. Then the hermaphrodite produces eggs, sel-ertilizes them using her internally stored sperm, and lays them. In the absence o males, egg production ceases ater all the sperm are utilized. We assume that the juvenile growth period occurs during 0 < < , spermatogenesis occurs during < < +, and egg-production, sel-ertilization, and egg laying occurs during
2.7. THE BROOD SIZE OF A HERMAPHRODITIC WORM
growth period sperm production period egg production period sperm production rate egg production rate brood size
31
72 h 11.9 h 65 h 24 h−1 4.4 h−1 286
Table 2.2: Parameters in a lie-history model o C. elegans , with experimental estimates. + < < + + . Here, we want to understand why hermaphrodites limit their sperm production. Biologists dene males and emales rom the size and metabolic cost o their gametes: sperm are cheap and eggs are expensive. So on rst look, it is puzzling why the total number o ofspring produced by a hermaphrodite is limited by the number o sperm produced, rather than by the number o eggs. There must be a hidden cost to the hermaphrodite o producing additional sperm other than metabolic. To understand the basic biology, it is instructive to consider two limiting cases: (1) no sperm production; (2) innite sperm production. In both cases, the hermaphrodite produces no ofspring—in the rst case because there are no sperm, and in the second case because there are no eggs. The number o sperm produced by a hermaphrodite beore laying eggs is thereore a compromise; although more sperm means more ofspring, more sperm also means delayed egg production. Our main theoretical assumption is that natural selection will avor worms with the ability to establish populations with the largest Malthusian parameter . Worms containing a genetic mutation resulting in a larger value or will eventually outnumber all other worms. The parameters we need or our mathematical model are listed in Table 2.2, together with estimated experimental values (Cutter, 2004). In addition to the growth period , sperm production period , and egg production period (all in units o hours), we need the sperm production rate and the egg production rate (both in units o inverse hours). We also dene the brood size as the total number o ertilized eggs laid by a selng hermaphrodite. The brood size is equal to the number o sperm produced, and also equal to the number o eggs laid, so that = = . (2.18) We may use (2.18) to eliminate and in avor o : = /,
= /.
(2.19)
The continuous Euler-Lotka equation (2.12) or requires a model or () = ()(), where () is the age-specic maternity unction and () is the agespecic survival unction. The unction () satises the diferential equation (2.15), and here we make the simpliying assumption that the age-specic mortality unction () = , where is the age-independent per capita death rate. Implicitly, we are assuming that worms do not die o old age during egg-laying, but rather die o predation, starvation, disease, or other age-independent causes. Such an assumption is reasonable since worms can live in the laboratory or
32
CHAPTER 2. AGE-STRUCTURED AGE-STRUCTURED POPULATIONS POPULATIONS
several several weeks weeks ater sperm depletion. depletion. Solving Solving (2.15 ( 2.15)) with the initial condition (0) = 1 results in () = exp exp ( ). (2.20)
− ·
The age-specic maternity unction () is dened such that ()Δ )Δ is the expected number o ofspring produced over the age interval interval Δ Δ. We assume that a hermaphrodite lays eggs at a constant rate over the ages + < < + +; thereore, or + < < + + , () = (2.21) 0 othe otherrwis wise.
Using (2.19 (2.19), ), (2.20 2.20)) and (2.21), (2.21), the continuous Euler-Lotka equation (2.12 ( 2.12)) or the Malthusian parameter becomes +/ +/
exp
+/
Integrating,
[−
]
( + ) = 1.
(2.22)
+/ +/
1=
exp
+/
[−
] ] − [− ] − [− ] − [−
( + )
exp ( + /)( /)( + ) exp ( + / + /)( / )( + ) + = exp ( + /)( / )( + ) 1 exp (/)( / )( + ) , + =
[− [−
which may be rewritten as
[
( + )exp ( + /)( / )( + ) = 1
exp
]}
]}
(/)( / )( + )
]}
.
(2.23)
With the parameters parameters , , , and xed, the integrated Euler-Lotka equation (2.23 2.23)) is an implicit equation or = ( ). To demonstrate that = ( ) has a maximum at some value o , we numerically solve (2.23 ( 2.23)) or with the parameter values o , and obtained rom Table 2.2. Since + is maximum at the same value o that is maximum, and only enters (2.23 (2.23)) in the orm + , without loss o generality we can take = 0. To solve (2.23 (2.23), ), it is best to make use o Newton’s method. We let
[
( () = ( ( + )exp ( + /)( / )( + )
] − − [− 1
exp
]}
(/)( / )( + )
and diferentiate with respect to to obtain
′ () = 1 + ( + /)( / )( + ) exp ( + /)( / )( + )
[
] [
For a given , we then solve ( () = 0 by iterating +1 =
]−
exp
[−
,
]
(/)( / )( + ) .
( ) . − ( ( ) ′
Using appropriate starting values or , the unction = ( ) can be computed and is presented in Fig. 2.6. 2.6. Evidently, Evidently, is maximum near the value = 152, which is 53% o the experimental value or shown in Table 2.2.
2.7. THE BROOD BROOD SIZE OF A HERMAPHRO HERMAPHRODITI DITIC C WORM
33
0.056
0.054
0.052
0.05 r
0.048
0.046
0.044
0.042 0
100
200
300
400
500
B
Figure 2.6: A plot o o = ( ), which demonstrates that the Malthusian growth rate is maximum near a brood size o = 152. We can also directly determine a single equation or the value o at which is maximum. We implicitly diferentiate (2.23 ( 2.23)) with respect to —with the only parameter that depends on —and apply the condition / = 0. We nd ( + )exp ( + /)( (2.24) / )( + ) = exp (/)( / )( + ) .
[
]
[− ]− }
Taking the ratio o (2.23 ( 2.23)) to (2.24 (2.24)) results in 1= exp (/)( / )( + ) rom which we can nd
[
1 ,
ln(1 + / + /)). Substituting (2.26) (2.26) back into either (2.23 (2.23)) or (2.24 (2.24)) results in +=
1+
+
= . +
ln 1 +
]
(2.25)
(2.26)
(2.27)
Equation (2.27 (2.27)) contains the our parameters , , and , which can be urther urther reduced reduced to three parameters parameters by dimensiona dimensionall analysis. analysis. The brood size is already dimensionless. The rate at which eggs are produced and laid may be multiplied by the sperm production period = / to orm the dimensionless parameter = /. /. The parameter represents the number o eggs orgone because o the adult sperm production period and is a measure o the cost o producing sperm. Similarly, may be multiplied by the larval growth period to orm the dimensionless parameter = . . The parameter parameter represents the number o eggs orgone because o the juvenile growth period and is a measure o the cost o development. With , and as our three dimensionless parameters, (2.27 2.27)) becomes 1
+
1+
ln 1 +
=
. +
(2.28)
34
CHAPTER 2. AGE-STRUCTURED AGE-STRUCTURED POPULATIONS POPULATIONS
700
600
500
y 400
300
200
5
10
15
20
25
30
35
40
45
50
55
x
Figure Figure 2.7: Solution Solution curve o versus with brood size = 286, obtained by solving (2.28 (2.28). ). Values or or , , and are taken rom Table 2.2. The cross, crossed crossed open circle circle and open circle correspond to = and = / with = 1, 1/ 1/3 and 1/ 1/8, respectively.
sperm-
0
juvenile
−
-
eggs
+
+ +
adult
-
Figure 2.8: A more rened timeline o a hermaphrodite’s lie. Given values or two o the three dimensionless parameters , and , (2.28 2.28)) may be solved or the remaining parameter, either explicitly or the case o = (, ), or by Newton’s method. The values o and obtained rom Table 2.2 are = 52 52..5 and = 317. With = 286, the solution = () is shown in Fig. 2.7, 2.7, with the experimental value o (, (, ) plotted as a cross. The seemingly large disagreement between the theoretical result and the experimental data leads us to question the underlying assumptions o the model. Indeed, Cutter (2004) rst suggested that the sperm produced precociously as a juvenile does not delay egg production and should be considered cost ree. One possibility is to x the absolute number o sperm produced precociously and to optimize optimize the number o sperm produced produced as an adult. Another Another possibility possibility is to x the raction o sperm produced precociously and to optimize the total number number o sperm produced. produced. This latter assumption assumption was made by Cutter (2004) and seems to best improve the model agreement with the experimental data. We thereore split the total sperm production period into juvenile and adult sperm production periods and , with = + . The The revi revise sed d
35
2.7. THE BROOD SIZE OF A HERMAPHRODITIC WORM
timeline o the hermaphrodite’s lie is now shown in Fig. 2.8. With the raction o sperm produced as an adult denoted by , and the raction produced as a juvenile by 1 , we have
−
= (1
− ),
= .
(2.29)
Equation (2.21) or the age-specic maternity unction becomes () =
or + < < + + , 0 otherwise.
(2.30)
with = /.
(2.31)
The Euler-Lotka equation (2.22) is then changed by the substitution / . Following this substitution to the nal result given by (2.28) shows that this equation still holds (which is in act equivalent to (12) o Chasnov (2011)), but with the now changed denition
→
= /.
(2.32)
A close examination o the results shown in Fig. 2.7 demonstrates that nearperect agreement can be made between the results o the theoretical model and the experimental data i = 1/8 (shown as the open circle in Fig. 2.7). Cutter (2004) suggested the value o = 1/3, and this result is shown as a crossed open circle in Fig. 2.7, still in much better agreement with the experimental data than the open circle corresponding to = 1. The added modeling o precocious sperm production through the parameter thus seems to improve the verisimilitude o the model to the underlying biology.
Reerences Barker, D.M. Evolution o sperm shortage in a selng hermaphrodite. Evolution (1992) 46, 1951-1955. Cutter, A.D. Sperm-limited ecundity in nematodes: How many sperm are enough? Evolution (2004) 58, 651-655. Chasnov, J.R. Evolution o increased sel-sperm production in postdauer hermaphroditic nematodes. Evolution (2011) 65, 2117-2122. Hodgkin, J. & Barnes, T.M. More is not better: brood size and population growth in a sel-ertilizing nematode. Proc. R. Soc. Lond. B. (1991) 246, 19-24. Charlesworth, B. Evolution in age-structured populations. (1980) Cambridge University Press.
36
CHAPTER 2. AGE-STRUCTURED POPULATIONS
Chapter 3
Stochastic Modeling o Population Growth Our derivation o the Malthusian growth model implicitly assumed a large population size. Smaller populations exhibit stochastic efects and these can considerably complicate the modeling. Since in general, modeling stochastic processes in biology is an important yet dicult topic, we will spend some time here analyzing the simplest model o births in nite populations.
3.1
A stochastic model o population growth
The size o the population is now considered to be a discrete random variable. We dene the time-dependent probability mass unction () o to be the probability that the population is o size at time . Since must take on one o the values rom zero to innity, we have ∞
() = 1,
=0
or all 0. Again, let be the average per capita birth rate. We make the simpliying approximations that all births are singlets, and that the probability o an individual giving birth is independent o past birthing history. We can then interpret probabilistically by supposing that as Δ 0, the probability that an individual gives birth during the time Δ is given by Δ. For example, i the average per capita birthrate is one ofspring every 365-day year, then the probability that a given individual gives birth on a given day is 1 /365. As we will be considering the limit as Δ 0, we neglect probabilities o more than one birth in the population in the time interval Δ since they are o order (Δ)2 or higher. Furthermore, we will suppose that at = 0, the population size is known to be 0 , so that 0 (0) = 1, with all other ’s at = 0 equal to zero. We can determine a system o diferential equations or the probability mass unction () as ollows. For a population to be o size > 0 at a time + Δ, either it was o size 1 at time and one birth occurred, or it was o size at time and there were no births; that is
≥
→
→
−
( + Δ) = −1 ()(
− 1)Δ + ()(1 − Δ).
37
38
CHAPTER 3. STOCHASTIC POPULATION GROWTH
Subtracting () rom both sides, dividing by Δ, and taking the limit Δ results in the orward Kolmogorov diferential equations , = (
[
]
− 1) 1 − , −
= 1, 2, . . . ,
→0
(3.1)
where 0 () = 0 (0) since a population o zero size remains zero. This system o coupled, rst-order, linear diferential equations can be solved iteratively. We rst review how to solve a rst-order linear diferential equation o the orm + = (), (0) = 0 , (3.2) where = () and is constant. First, we look or an integrating actor such that () = + .
Diferentiating the let-hand-side and multiplying out the right-hand-side results in + = + ; and canceling terms yields = . We may integrate this equation with an arbitrary initial condition, so or simplicity we take (0) = 1. Thereore, () = . Hence, = ().
( )
Integrating this equation rom 0 to t yields
()
− (0) =
().
0
Thereore, the solution is −
() =
(0) +
() .
0
(3.3)
The orward Kolmogorov diferential equation (3.1) is o the orm (3.2) with = and () = ( 1) −1 . With the population size known to be 0 at = 0, the initial conditions can be written succinctly as (0) = , 0 , where is the Kronecker delta, dened as
−
=
0, 1,
i = ; i = .
̸
Thereore, ormal integration o (3.1) using (3.3) results in
() = − , 0 + (
− 1)
0
−1 () .
(3.4)
39
3.1. A STOCHASTIC MODEL OF POPULATION GROWTH
The rst ew solutions o (3.4) can now be obtained by successive integrations: 0, i < 0 ; − 0 , i = 0 ; − 0 − () = 0 [1 ], i = 0 + 1; 1 − 0 − 2 [1 ] , i = 0 + 2; 2 0 ( 0 + 1) ..., i . . . .
⎧⎪ ⎪⎨ ⎪⎪⎪ ⎩
−
−
Although we will not need this, or completeness I give the complete solution. By dening the binomial coecient as the number o ways one can select k objects rom a set o n identical objects, where the order o selection is immaterial, we have ! = , !( )!
−
(read as “ choose ”). The general solution or (), () =
≥ 0, is known to be
− [− ] − 0
1 − 0 1 1
−
− 0
,
which statisticians call a shited negative binomial distribution . The determination o the time-evolution o the probability mass unction o completely solves this stochastic problem. O usual main interest is the mean and variance o the population size, and although both could in principle be computed rom the probability mass unction, we will compute them directly rom the diferential equation or . The denitions o the mean population size and its variance 2 are
⟨ ⟩
∞
⟨ ⟩ =
∞
2
,
=
=0
(
=0
− ⟨ ⟩ 2 ,
)
(3.5)
and we will make use o the equality 2 = 2
⟨ ⟩ − ⟨ ⟩2.
(3.6)
Multiplying the diferential equation (3.1) by the constant , summing over , and using = 0 or < 0 , we obtain =
⟨ ⟩
∞
∞
(
= 0 +1
− 1)
−1
−
2 .
= 0
Now ( 1) = ( 1)2 +( 1), so that the rst term on the right-hand-side is
−
−
−
∞
= 0 +1
∞
(
− 1)
−1
=
∞
(
= 0 +1
− 1)
∞
=
= 0
2
−1
+
+
= 0
(
= 0 +1
∞
2
,
− 1)
−1
40
CHAPTER 3. STOCHASTIC POPULATION GROWTH
where the second equality was obtained by shiting the summation index downward by one. Thereore, we nd the amiliar Malthusian growth equation = .
⟨ ⟩
⟨ ⟩
Together with the initial condition (0) = 0 , we can nd the solution
⟨ ⟩ ⟨ ⟩() = 0.
(3.7)
We proceed similarly to nd 2 by rst determining the diferential equation or 2 . Multiplying the diferential equation or , (3.1), by 2 and summing over results in
⟨ ⟩
2 =
⟨ ⟩
∞
∞
2 (
= 0 +1
− 1)
−1
−
3 .
= 0
Here, we use the equality 2 ( 1) = ( 1)3 +2( 1)2 +( 1). Proceeding in the same way as above by shiting the index downward, we obtain
−
−
−
−
2
⟨ ⟩ − 2⟨ 2⟩ = ⟨ ⟩.
(3.8)
Since is known, (3.8) is a rst-order, linear, inhomogeneous equation or 2 , which can be solved using an integrating actor. The solution obtained using (3.3) is
⟨ ⟩
⟨ ⟩
2
⟨ ⟩ =
2
02
+
−2
0
⟨ ⟩()
,
with given by (3.7). Perorming the integration, we obtain
⟨ ⟩
⟨ 2⟩ = 2 02 + 0(1 − ) . Finally, using 2 = ⟨ 2 ⟩ −⟨ ⟩2 , we obtain the variance. Thus we arrive at our −
[
]
nal result or the population mean and variance:
⟨ ⟩ = 0,
2 = 0 2 1
− .
(− )
(3.9)
The coecient o variation measures the standard deviation relative to the mean, and is here given by = / =
⟨ −⟩ 1
− . 0
√
For large , the coecient o variation thereore goes like 1 / 0 , and is small when 0 is large. In the next section, we will determine the limiting orm o the probability distribution or large 0 , recovering both the deterministic model and a Gaussian model approximation.
41
3.2. ASYMPTOTICS OF LARGE INITIAL POPULATIONS
3.2
Asymptotics o large initial populations
Our goal here is to solve or an expansion o the distribution in powers o 1/ 0 to leading-orders; notice that 1/ 0 is small i 0 is large. To zeroth-order, that is in the limit 0 , we will show that the deterministic model o population growth is recovered. To rst-order in 1/ 0 , we will show that the probability distribution is normal. The latter result will be seen to be a consequence o the well-known Central Limit Theorem in probability theory. We develop our expansion by working directly with the diferential equation or (). Now, when the population size is a discrete random variable (taking only the nonnegative integer values o 0, 1, 2, . . . ), () is the probability mass unction or . I 0 is large, then the discrete nature o is inconsequential, and it is preerable to work with a continuous random variable and its probability density unction. Accordingly, we dene the random variable = / 0 , and treat as a continuous random variable, with 0 < . Now, () is the probability that the population is o size at time , and the probability density unction o , (, ), is dened such that (, ) is the probability that . The relationship between and can be determined by considering how to approximate a discrete probability distribution by a continuous distribution, that is, by dening such that
→∞
≤
∞
∫
≤
≤
( + 12 )/ 0
() =
(, )
( − 12 )/ 0
= (/ 0 , )/ 0 where the last equality becomes exact as 0 denition or (, ) is given by (, ) = 0 (),
→ ∞. Thereore, the appropriate = / 0 ,
(3.10)
which satises ∞
∞
(, ) =
0
(/ 0 , )(1/ 0 )
=0 ∞
=
()
=0
= 1,
the rst equality (exact only when 0 ) being a Reimann sum approximation o the integral. We now transorm the innite set o ordinary diferential equations ( 3.1) or () into a single partial diferential equation or (, ). We multiply (3.1) by 0 and substitute = 0 , () = (, )/ 0 , and −1 () = ( 10 , )/ 0 to obtain
→∞
−
(, ) = ( 0
− 1) ( −
We next Taylor series expand (
1 , ) 0
− 0 (, )
.
(3.11)
− 1/ 0, ) around , treating 1/ 0 as a small
42
CHAPTER 3. STOCHASTIC POPULATION GROWTH
parameter. That is, we make use o (
− 10 , ) = (, ) − 10 (, ) + 2 1 2 (, ) − . . . 0 (−1) = . ∞
! 0
=0
The two leading terms proportional to 0 on the right-hand-side o (3.11) cancel exactly, and i we group the remaining terms in powers o 1 / 0 , we obtain or the rst three leading terms in the expansion = =
− −
1 ( + ) + + 2 2! 1! 0 1 1 ( ) ( ) + 2 ( ) 0 2! 0 3!
−
1 0
−
+ 3! 2!
− ...
− ...
; (3.12)
and higher-order terms can be obtained by ollowing the evident pattern. Equation (3.12) may be urther analyzed by a perturbation expansion o the probability density unction in powers o 1/ 0 : (, ) = (0) (, ) +
1 (1) 1 (, ) + 2 (2) (, ) + . . . . 0 0
(3.13)
Here, the unknown unctions (0) (, ), (1) (, ), (2) (, ), etc. are to be determined by substituting the expansion (3.13) into (3.12) and equating coecients o powers o 1/ 0 . We thus obtain or the coecients o (1/ 0 )0 and (1/ 0 )1 , (0)
=
(1)
=
3.2.1
− − − (0)
(1)
,
1 (0) 2
(3.14) .
(3.15)
Derivation o the deterministic model
The zeroth-order term in the perturbation expansion (3.13), (0) (, ) = lim (, ), 0 →∞
satises (3.14). Equation (3.14) is a rst-order linear partial diferential equation with a variable coecient. One way to solve this equation is to try the ansatz (0) (, ) = ℎ() (),
= (, ),
(3.16)
together with the initial condition (0) (, 0) = (), since at = 0, the probability distribution is assumed to be a known unction. This initial condition imposes urther useul constraints on the unctions ℎ() and (, ): ℎ(0) = 1, (, 0) = .
43
3.2. ASYMPTOTICS OF LARGE INITIAL POPULATIONS
The partial derivatives o (3.16) are (0)
= ℎ′ + ℎ ′ ,
(0) = ℎ ′ ,
which ater substitution into (3.14) results in (ℎ′ + ℎ) + ( + ) ℎ ′ = 0. This equation can be satised or any provided that ℎ() and (, ) satisy ℎ′ + ℎ = 0,
+ = 0.
(3.17)
The rst equation or ℎ(), together with the initial condition ℎ(0) = 1, is easily solved to yield ℎ() = − . (3.18) To determine a solution or (, ), we try the technique o separation o variables. We write (, ) = () (), and upon substitution into the diferential equation or (, ), we obtain ′ + ′ = 0; and division by XT and separation yields ′ =
′
. −
(3.19)
Since the let-hand-side is independent o , and the right-hand-side is independent o , both the let-hand-side and right-hand-side must be constant, independent o both and . Now, our initial condition is (, 0) = , so that () (0) = . Without lose o generality, we can take (0) = 1, so that () = . The right-hand-side o (3.19) is thereore equal to the constant , and we obtain the diferential equation and initial condition
−
′ + = 0,
(0) = 1,
which we can solve to yield () = − . Thereore, our solution or (, ) is (, ) = − .
(3.20)
By putting our solutions (3.18) and (3.20) together into our ansatz (3.16), we have obtained the general solution to the pde: (0) (, ) = − (− ). To determine , we apply the initial condition o the probability mass unction, (0) = , 0 . From (3.10), the corresponding initial condition on the probability distribution unction is (, 0) =
0 0
1 i 1 2 0 otherwise.
−
≤ ≤ 1 + 2 1 , 0
In the limit 0 , (, 0) (0) (, 0) = ( 1), where ( 1) is the Dirac delta-unction, centered around 1. The delta-unction is widely used in
→∞
→
−
−
44
CHAPTER 3. STOCHASTIC POPULATION GROWTH
quantum physics and was introduced by Dirac or that purpose. It now nds many uses in applied mathematics. It can be dened by requiring that, or any unction (), +∞
() () = (0).
−∞
The usual view o the delta-unction ( ) is that it is zero everywhere except at = at which it is innite, and its integral is one. It is not really a unction, but it is what mathematicians call a distribution. Now, since (0) (, 0) = () = ( 1), our solution becomes
−
−
(0) (, ) = − (−
− 1). This can be rewritten by noting that (letting = − ), +∞
() (
−∞
−
(3.21)
1 +∞ ) = ( + )/ () −∞ 1 = (/),
(
yielding the identity (
)
− ) = 1 ( − ).
From this, we can rewrite our solution (3.21) in the more intuitive orm (0) (, ) = (
− ).
(3.22)
Using (3.22), the zeroth-order expected value o is ∞
⟨0⟩ = =
0
(0) (, )
∞
(
0
− )
= ; while the zeroth-order variance is 20 = 20
⟨ ⟩ − ⟨0⟩2 ∞
=
2 (0) (, )
0
− 2
∞
=
0 2
= = 0.
2 (
− ) − 2
− 2
Thus, in the innite population limit, the random variable has zero variance, and is thereore no longer random, but ollows = deterministically. We say that the probability distribution o becomes sharp in the limit o large population sizes. The general principle o modeling large populations deterministically can simpliy mathematical models when stochastic efects are unimportant.
45
3.2. ASYMPTOTICS OF LARGE INITIAL POPULATIONS
3.2.2
Derivation o the normal probability distribution
We now consider the rst-order term in the perturbation expansion (3.13), which satises (3.15). We do not know how to solve (3.15) directly, so we will attempt to nd a solution ollowing a more circuitous route. First, we proceed by computing the moments o the probability distribution. We have ∞
⟨ ⟩ =
(, )
0
∞
=
(0)
0
= 0
⟨ ⟩
∞
1 (, ) + 0
(1) (, ) + . . .
0
1 + +..., 0 1
⟨ ⟩
where the last equality denes 0 , etc. Now, using (3.22),
⟨ ⟩
∞
0
⟨ ⟩= =
=
0
(0) (, )
∞
(
0
(3.23)
− )
.
To determine 1 , we use (3.15). Multiplying by and integrating, we have
⟨ ⟩
1 =
⟨ ⟩ −
∞
(1)
0
−
1 2
∞
(0)
0
.
(3.24)
We integrate by parts to remove the derivatives o , assuming that and all its derivatives vanish on the boundaries o integration, where is equal to zero or innity. We have ∞
(1)
0
∞
=
−
(1)
0
=
1 ,
−⟨ ⟩
and 1 2
∞
= 2 ( =
(0)
0
−
=
(
∞
0
− 1)
2
−1 (0) ∞
−1 (0)
0
− 1) ⟨ 1⟩. 0 2 −
Thereore, ater integration by parts, (3.24) becomes 1 =
⟨ ⟩ − ⟨⟩ + ( − 1) ⟨ 1⟩ 1 0 2 −
.
(3.25)
Equation (3.25) is a rst-order linear inhomogeneous diferential equation and can be solved using an integrating actor (see (3.3) and the preceding discussion).
46
CHAPTER 3. STOCHASTIC POPULATION GROWTH
Solving this diferential equation using (3.23) and the initial condition 1 (0) = 0, we obtain ( 1) 1 = 1 − . (3.26) 2 The probability distribution unction, accurate to order 1/ 0 , may be obtained by making use o the so-called moment generating unction Ψ(), dened as
⟨ ⟩
−
⟨ ⟩
(− )
Ψ() =
⟨ ⟩
2 2 3 3 = 1+ + + + ... 2! 3! ∞ = . ! =0
⟨⟩ ⟨ ⟩ ⟨ ⟩
⟨ ⟩
To order 1/ 0 , we have ∞
∞
0 1 Ψ() = + ! 0 =0
1 + O(1/ 02 ). ! =0
⟨ ⟩
⟨ ⟩
(3.27)
Now, using (3.23), ∞
∞
0 = ! ! =0 =0
( )
⟨ ⟩
= ,
and using (3.26), ∞
∞
1 1 = 1 ! 2 =0
⟨ ⟩
=
(− (− (− (− (− (−
1 1 2
1 = 1 2 1 = 1 2 1 1 2 1 = 1 2 =
Thereore,
Ψ() =
−
1 ( ! =0
) − ( ) ) − ( ) ) ( ) ) ( ) ) ) (− ) 1)
∞
− 2 2
1 ( ! =0
1)
−2
∞
−
−
2
1 2 2 ! =0
2 2
∞
!
2
− 2
2
2
=0
− 2 2 .
1+
1 1 2 0
− 2 2 + . . .
.
(3.28)
We can recognize the parenthetical term o (3.28) as a Taylor-series expansion o an exponential unction truncated to rst-order, i.e., exp
(− ) 1 1 2 0
− 2 2
= 1+
1 1 2 0
− 2 2 + O(1/ 02 ).
(− )
47
3.3. SIMULATION OF POPULATION GROWTH
Thereore, to rst-order in 1/ 0 , we have
(− )
1 Ψ() = exp + 1 2 0
− 2 2
+ O(1/ 02 ).
(3.29)
Standard books on probability theory (e.g., A rst course in probability by Sheldon Ross, pg. 365) detail the derivation o the moment generating unction o a normal random variable: Ψ() = exp
⟨ ⟩
1 + 2 2 , 2
or a normal random variable;
(3.30)
and comparing (3.30) with (3.29) shows us that the probability distribution (, ) to rst-order in 1/ 0 is normal with the mean and variance given by
⟨⟩ = ,
2 =
1 2 1 0
− .
(− )
(3.31)
The mean and variance o = / 0 is equivalent to those derived or in (3.9), but now we learn that is approximately normal or large populations. The appearance o a normal probability distribution (also called a Gaussian probability distribution) in a rst-order expansion is in act a particular case o the Central Limit Theorem, one o the most important and useul theorems in probability and statistics. We state here a simple version o this theorem without proo:
Central Limit Theorem: Suppose that 1 , 2 , . . . , are independent and 2 . identically distributed (iid) random variables with mean and variance Then or suciently large , the probability distribution o the average o the ’s, denoted as the random variable = 1 =1 , is well approximated by 2 /. a Gaussian with mean and variance 2 =
⟨ ⟩
⟨ ⟩
∑
The central limit theorem can be applied directly to our problem. Consider that our population consists o 0 ounders. I () denotes the number o individuals descendent rom ounder at time (including the still living ounder), 0 then the total number o individuals at time is () = =1 (); and the average number o descendants o a single ounder is () = ()/ 0 . I the 2 , then mean number o descendants o a single ounder is , with variance by applying the central limit theorem or large 0 , the probability distribution unction o is well approximated by a Gaussian with mean = and 2 / . Comparing with our results (3.31), we nd = variance 2 = 0 2 = 2 (1 and − ).
⟨ ⟩
∑
−
3.3
⟨⟩ ⟨ ⟩ ⟨ ⟩
Simulation o population growth
As we have seen, stochastic modeling is signicantly more complicated than deterministic modeling. As the modeling becomes more detailed, it may become necessary to solve a stochastic model numerically. Here, or illustration, we show how to simulate individual realizations o population growth. A straightorward approach would make use o the birth rate directly. During the time Δ, each individual o a population has probability Δ o
48
CHAPTER 3. STOCHASTIC POPULATION GROWTH
giving birth, accurate to order Δ. I we have individuals at time , we can compute the number o individuals at time + Δ by computing random deviates (random numbers between zero and unity) and calculating the number o new births by the number o random deviates less than Δ. In principle, this approach will be accurate provided that Δ is suciently small. There is, however, a more ecient way to simulate population growth. We rst determine the probability density unction o the interevent time , dened as the time required or the population to grow rom size to size +1 because o a single birth. A simulation rom population size 0 to size would then simply require computing 0 diferent random values o , a relatively easy and quick computation. Let ( ) be the probability that the interevent time or a population o size lies in the interval (, + ), ( ) = 0 ( ′ ) ′ the probability that the interevent time is less than , and ( ) = 1 ( ) the probability that the interevent time is greater than . ( ) is called the probability density unction o , and ( ) the cumulative distribution unction o . They satisy the relation ( ) = ′ ( ). Since the probability that the interevent time is greater than + Δ is given by the probability that it is greater than times the probability that no new births occur in the time interval ( , + Δ ), the distribution ( + Δ ) is given or small Δ by
−
∫
( + Δ ) = ( )(1 Diferencing and taking the limit Δ =
−
− Δ ).
→ 0 yields the diferential equation
−,
which can be integrated using the initial condition (0) = 1 to urther yield ( ) = − . Thereore, ( ) = 1
−
−
,
( ) = − ;
(3.32)
( ) has the orm o an exponential distribution. We now compute random values or using its known cumulative distribution unction ( ). We do this in reerence to Figure 3.1, where we plot = ( ) versus . The range o ( ) is [0, 1), and or each , we can compute by inverting ( ) = 1 exp( ); that is,
−
−
() =
− ln(1 − ) .
Our claim is that with a random number sampled uniormly on (0, 1), the random numbers () have probability density unction ( ). Now, the raction o random numbers in the interval (1 , 2 ) is equal to the corresponding raction o random numbers in the interval ( 1 , 2 ), where 1 = (1 ) and 2 = (2 ); and or sampled uniormly on (0, 1), this raction is equal to 2 1 . We thus have the ollowing string o equalities: (the raction o random numbers in (1 , 2 )) =(the raction o random numbers in ( 1 , 2 ))= 2 1 = ( 2 ) ( 1 ) = 12 ( ) . Thereore, (the raction o random
−
−
−
∫
49
3.3. SIMULATION OF POPULATION GROWTH
1
y = F(τ) = 1 − e
−bNτ
y =F(τ ) 2
y
2
=F(τ )
1
1
0 τ
1
τ
2
Figure 3.1: How to generate a random number with a given probability density unction.
∫
numbers in ( 1 , 2 )) = 12 ( ) . This is exactly the denition o having probability density unction ( ). Below, we illustrate a simple MATLAB unction that simulates one realization o population growth rom initial size 0 to nal size , with birth rate . function [t, N] = population_growth_simulation(b,N0,Nf) % simulates population growth from N0 to Nf with birth rate b N=N0:Nf; t=0*N; y=rand(1,Nf-N0); tau=-log(1-y)./(b*N(1:Nf-N0)); % interevent times t=[0 cumsum(tau)]; % cumulative sum of interevent times
The unction population growth simulation.m can be driven by a MATLAB script to compute realizations o population growth. For instance, the ollowing script computes 100 realizations or a population growth rom 10 to 100 with = 1 and plots all the realizations: % calculate nreal realizations and plot b=1; N0=10; Nf=100; nreal=100; for i=1:nreal [t,N]=population_growth_simulation(b,N0,Nf); plot(t,N); hold on; end xlabel(’t’); ylabel(’N’);
Figure 3.2 presents three graphs, showing 25 realizations o population growth starting with population sizes o 10, 100, and 1000, and ending with population
50
CHAPTER 3. STOCHASTIC POPULATION GROWTH
Figure 3.2: Twenty-ve realizations o population growth with initial population sizes o 10, 100, and 1000, in (a), (b), and (c), respectively. sizes a actor o 10 larger. Observe that the variance, relative to the initial population size, decreases as the initial population size increases, ollowing our analytical result (3.31).
Chapter 4
Inectious Disease Modeling In the late 1320’s, an outbreak o the bubonic plague occurred in China. This disease is caused by the bacteria Yersinia pestis and is transmitted rom rats to humans by eas. The outbreak in China spread west, and the rst major outbreak in Europe occurred in 1347. During a ve year period, 25 million people in Europe, approximately 1/3 o the population, died o the black death. Other more recent epidemics include the inuenza pandemic known as the Spanish u killing 50-100 million people worldwide during the year 1918-1919, and the present AIDS epidemic, originating in Arica and rst recognized in the USA in 1981, killing more than 25 million people. For comparison, the SARS epidemic or which Hong Kong was the global epicenter resulted in 8096 known SARS cases and 774 deaths. Yet, we know well that this relatively small epidemic caused local social and economic turmoil. Here, we introduce the most basic mathematical models o inectious disease epidemics and endemics. These models orm the basis o the necessarily more detailed models currently used by world health organizations, both to predict the uture spread o a disease and to develop strategies or containment and eradication.
4.1
The SI model
The simplest model o an inectious disease categorizes people as either susceptible or inective ( ). One can imagine that susceptible people are healthy and inective people are sick. A susceptible person can become inective by contact with an inective. Here, and in all subsequent models, we assume that the population under study is well mixed so that every person has equal probability o coming into contact with every other person. This is a major approximation. For example, while the population o Amoy Gardens could be considered well mixed during the SARS epidemic because o shared water pipes and elevators, the population o Hong Kong as a whole could not because o the larger geographical distances, and the limited travel o many people outside the neighborhoods where they live. We derive the governing diferential equation or the SI model by considering the number o people that become inective during time Δ. Let Δ be the probability that a random inective person inects a random susceptible 51
52
CHAPTER 4. INFECTIOUS DISEASE MODELING
person during time Δ. Then with susceptible and inective people, the expected number o newly inected people in the total population during time Δ is Δ . Thus, ( + Δ) = () + Δ () (), and in the limit Δ
→ 0,
= .
(4.1)
We diagram (4.1) as
−−→ .
Later, diagrams will make it easier to construct more complicated systems o equations. We now assume a constant population size , neglecting births and deaths, so that + = . We can eliminate rom (4.1) and rewrite the equation as = 1 ,
−
which can be recognized as a logistic equation, with growth rate and carrying capacity . Thereore as and the entire population will become inective.
→
4.2
→∞
The SIS model
The SI model may be extended to the SIS model, where an inective can recover and become susceptible again. We assume that the probability that an inective recovers during time Δ is given by Δ. Then the total number o inective people that recover during time Δ is given by Δ, and
× ( + Δ) = () + Δ () () − Δ (),
or as Δ
→ 0,
=
− ,
(4.2)
which we diagram as −− −−
↽ ⇀ . Using + = , we eliminate rom (4.2) and dene the basic reproductive ratio as . (4.3) 0 =
ℛ
Equation(4.2) may then be rewritten as = (
ℛ0 − 1)
− 1
(1 1/
− ℛ0)
,
which is again a logistic equation, but now with growth rate ( 0 1) and carrying capacity (1 1/ 0 ). The disease will disappear i the growth rate is negative, that is, 0 < 1, and it will become endemic i the growth rate is
ℛ
− ℛ
ℛ −
53
4.2. THE SIS MODEL
positive, that is, 0 > 1. For an endemic disease with 0 > 1, the number o inected people approaches the carrying capacity: (1 1/ 0 ) as . We can give a biological interpretation to the basic reproductive ratio 0 . Let () be the probability that an individual initially inected at = 0 is still inective at time . Since the probability o being inective at time +Δ is equal to the probability o being inective at time multiplied by the probability o not recovering during time Δ, we have
ℛ
→
( + Δ) = ()(1 or as Δ
→ 0,
=
ℛ
− ℛ
→∞ ℛ
− Δ),
−.
With initial condition (0) = 1, () = − .
(4.4)
Now, the expected number o secondary inections produced by a single primary inective over the time period (, + ) is given by the probability that the primary inective is still inectious at time multiplied by the expected number o secondary inections produced by a single inective during time ; that is, () (). We assume that the total number o secondary inections rom a single inective individual is small relative to the population size . Thereore, the expected number o secondary inectives produced by a single primary inective introduced into a completely susceptible population is
×
∞
∞
() ()
0
≈
()
0
∞
=
−
0
= = 0,
ℛ
where we have approximated () during the time period in which the inective remains inectious. I a single inected individual introduced into a completely susceptible population produces more than one secondary inection beore recovering, then 0 > 1 and the disease becomes endemic. We have seen previously two other analogous calculations with results similar to (4.4): the rst when computing the probability distribution o a worm’s lietime ( 2.7); and the second when computing the probability distribution o the interevent time between births ( 3.3). We have also seen an analogous denition o the basic reproductive ratio in our previous discussion o age-structured populations ( 2.5). There, the basic reproductive ratio was the number o emale ofspring expected rom a new born emale over her lietime; the population size would grow i this value was greater than unity. Both () and 0 are good examples o how the same mathematical concept can be applied to seemingly unrelated biological problems.
≈
ℛ
S
S
S
ℛ
54
4.3
CHAPTER 4. INFECTIOUS DISEASE MODELING
The SIR epidemic disease model
The SIR model, rst published by Kermack and McKendrick in 1927, is undoubtedly the most amous mathematical model or the spread o an inectious disease. Here, people are characterized into three classes: susceptible , inective and removed . Removed individuals are no longer susceptible nor inective or whatever reason; or example, they have recovered rom the disease and are now immune, or they have been vaccinated, or they have been isolated rom the rest o the population, or perhaps they have died rom the disease. As in the SIS model, we assume that inectives leave the class with constant rate , but in the SIR model they move directly into the class. The model may be diagrammed as ,
−→ −→
and the corresponding coupled diferential equations are =
−,
=
− ,
= ,
(4.5)
with the constant population constraint + + = . For convenience, we nondimensionalize (4.5) using or population size and −1 or time; that is, let ^ = /, ^ = /, ^ = /, ^ = , and dene the dimensionless basic reproductive ratio as . ℛ0 =
(4.6)
The dimensionless SIR equations are then given by ^ = ^
−ℛ0 ^,^
^ = ^
^ ^ = , ^
ℛ0 ^ ^ − ,^
(4.7)
^ ^+ ^ = 1. with dimensionless constraint + We will use the SIR model to address two undamental questions: (1) Under what condition does an epidemic occur? (2) I an epidemic occurs, what raction o a well-mixed population gets sick? ^* , ^* , ^ * ) be the xed points o (4.7). Setting / ^ ^ = / ^ ^ = Let ( ^ ^ = 0, we immediately observe rom the equation or / ^ ^ that ^ = 0, / ^ ^ and this value orces all the time-derivatives to vanish or any and . Since ^ = 0, we have ^ = 1 ^, evidently all the xed points o (4.7) are given with ^* , ^* , ^ * ) = ( ^* , 0, 1 ^* ). by the one parameter amily ( An epidemic occurs when a small number o inectives introduced into a susceptible population results in an increasing number o inectives. We can assume an initial population at a xed point o (4.7), perturb this xed point by introducing a small number o inectives, and determine the xed point’s stability. An epidemic occurs when the xed point is unstable. The linear ^ ^ in stability problem may be solved by considering only the equation or / ^ ^ ^ (4.7). With 1 and 0 , we have
−
≪
−
≈
^ = ^
ℛ
^ 0 0
−1
^ ,
55
4.3. THE SIR EPIDEMIC DISEASE MODEL fraction of population that get sick 0.8 0.7 0.6 0.5 ∞ 0.4 R
ˆ
0.3 0.2 0.1 0 0
0.5
1
1.5
2
R0
Figure 4.1: The raction o the population that gets sick in the SIR model as a unction o the basic reproduction ratio 0 .
ℛ
^0 1 > 0. With the basic reproductive ratio so that an epidemic occurs i 0 ^0 = 0 / , where 0 is the number o initial susceptible given by (4.6), and individuals, an epidemic occurs i
ℛ
−
ℛ0 ^0 = 0 > 1,
(4.8)
which could have been guessed. An epidemic occurs i an inective individual introduced into a population o 0 susceptible individuals inects on average more than one other person. We now address the second question: I an epidemic occurs, what raction o the population gets sick? For simplicity, we assume that the entire initial ^0 = 1. We expect the solution population is susceptible to the disease, so that o the governing equations (4.7) to approach a xed point asymptotically in time (so that the nal number o inectives will be zero), and we dene this ^ , ^ ) ^ = (1 ^ ∞ , 0, ^ ∞ ), with ^ ∞ equal to the raction xed point to be (, ^ ∞ , it is simpler to work with a o the population that gets sick. To compute ^ ^ = (/ ^ )( ^ / ^ ^), so transormed version o (4.7). By the chain rule, / that
−
^ ^ ^ / = ^ ^ ^ / ^ = 0 ,
−ℛ
which is separable. Separating and integrating rom the initial to nal conditions, ^ ^ ^ 1− ^ ^ = , 0 ^ 1 0
∞
−ℛ
∞
56
CHAPTER 4. INFECTIOUS DISEASE MODELING
which upon integration and simplication, results in the ollowing transcenden^∞: tal equation or ^ ∞ −ℛ0 ^ = 0, 1 (4.9)
−
−
∞
an equation that can be solved numerically using Newton’s method. We have ^ ∞) = 1 ( ^∞ ′ (
^∞ −ℛ0
− ^ − ) = −1 + ℛ0 ∞
^∞ −ℛ0
,
;
^ ∞ ) = 0 iterates and Newton’s method or solving ( ^ (∞+1) = ^ ∞
()
−
^∞ ) ( () ^∞ ′ ( )
(0) or xed 0 and a suitable initial condition or ∞ , which we take to be unity. My code or computing ∞ as a unction o 0 is given below, and the result is shown in Fig. 4.1. There is an explosion in the number o inections as 0 increases rom unity, and this rapid increase is a classic example o what is known more generally as a threshold phenomenon .
ℛ
ℛ
ℛ
function [R0, R_inf] = sir_rinf % computes solution of R_inf using Newton’s method from SIR model nmax=10; numpts=1000; R0 = linspace(0,2,numpts); R_inf = ones(1,numpts); for i=1:nmax R_inf = R_inf - F(R_inf,R0)./Fp(R_inf,R0); end plot(R0,R_inf); axis([0 2 -0.02 0.8]) xlabel(’$\mathcal{R}_0$’, ’Interpreter’, ’latex’, ’FontSize’,16) ylabel(’$\hat R_\infty$’, ’Interpreter’, ’latex’, ’FontSize’,16); title(’fraction of population that get sick’) %subfunctions function y = F(R_inf,R0) y = 1 - R_inf - exp(-R0.*R_inf); function y = Fp(R_inf,R0) y = -1 + R0.*exp(-R0.*R_inf);
4.4
Vaccination
Table 4.1 lists the diseases or which vaccines exist and are widely administered to children. Health care authorities must determine the raction o a population that must be vaccinated to prevent epidemics. We address this problem within the SIR epidemic disease model. Let be the raction o the population that is vaccinated and * the minimum raction required to prevent an epidemic. When > * , an epidemic can not occur. Since even non-vaccinated people are protected by the absence o epidemics, we say that the population has acquired herd immunity . We assume that individuals are susceptible unless vaccinated, and vaccinated individuals are in the removed class. The initial population is then modeled as
57
4.4. VACCINATION
Disease Diphtheria Haemophilus inuenzae type b (Hib) Hepatitis A
Description A bacterial respiratory disease A bacterial inection occurring primarily in inants A viral liver disease
Hepatitis B
Same as Hepatitis A
Measles
A viral respiratory disease
Mumps
A viral lymph node disease
Pertussis (whooping cough) Pneumococcal disease
A bacterial respiratory disease
Polio
A viral lymphatic and nervous system disease
Rubella (German measles) Tetanus (lock jaw)
A viral respiratory disease A bacterial nervous system disease
Varicella (chickenpox)
A viral disease in the Herpes amily A viral skin and mucous membrane disease
Human papillomavirus
A bacterial disease
Symptoms Sore throat and lowgrade ever Skin and throat inections, meningitis, pneumonia, sepsis, and arthritis Potentially none; yellow skin or eyes, tiredness, stomach ache, loss o appetite, or nausea Same as Hepatitis A
Complications Airway obstruction, coma, and death Death in one out o 20 children, and permanent brain damage in 10% 30% o the survivors usually none
Lie-long liver problems, such as scarring o the liver and liver cancer Rash, high ever, Diarrhea, ear inections, cough, runny nose, pneumonia, encephalitis, and red, watery eyes seizures, and death Fever, headache, Meningitis, inammamuscle ache, and tion o the testicles or swelling o the ovaries, inammation lymph nodes close o the pancreas and to the jaw deaness Severe spasms o Pneumonia, encephalitis, coughing and death, especially in inants High ever, cough, death and stabbing chest pains, bacteremia, and meningitis Fever, sore throat, Paralysis that can lead to nausea, headaches, permanent disability and stomach aches, stif- death ness in the neck, back, and legs Rash and ever or Birth deects i acquired two to three days by a pregnant woman Lockjaw, stifness in Death in one third o the the neck and ab- cases, especially people domen, and di- over age 50 culty swallowing A skin rash o Bacterial inection o the blister-like lesions skin, swelling o the brain, and pneumonia Warts, cervical can- The 5-year survival rate cer rom all diagnoses o cervical cancer is 72%
Table 4.1: Previously common diseases or which vaccines have been developed.
58
CHAPTER 4. INFECTIOUS DISEASE MODELING
^ , ^ ) ^ = (1 , 0, ). We have already determined the stability o this xed (, point to perturbation by a small number o inectives. The condition or an ^0 = 1 , an epidemic occurs i epidemic to occur is given by (4.8), and with
−
−
ℛ0(1 − ) > 1. Thereore, the minimum raction o the population that must be vaccinated to prevent an epidemic is 1 * = 1 .
− ℛ0
Diseases with smaller values o 0 are easier to eradicate than diseases with larger values 0 since a population can acquire herd immunity with a smaller raction o the population vaccinated. For example, smallpox with 0 4 has been eradicated throughout the world whereas measles with 0 17 still has occasional outbreaks.
ℛ
ℛ
ℛ ≈
4.5
ℛ ≈
The SIR endemic disease model
A disease that is constantly present in a population is said to be endemic . For example, malaria is endemic to Sub-Saharan Arica, where about 90% o malaria-related deaths occur. Endemic diseases prevail over long time scales: babies are born, old people die. Let be the birth rate and the diseaseunrelated death rate. We separately dene to be the disease-related death rate; is now the immune class. We may diagram a SIR model o an endemic disease as
6
-
?
-
?
?
and the governing diferential equations are =
− − ,
=
− ( + + ),
=
− ,
(4.10)
with = + + . In our endemic disease model, separately satises the diferential equation / = ( ) , (4.11)
−
−
and is not necessarily constant. A disease can become endemic in a population i / stays nonnegative; that is, () 1. + +
≥
For a disease to become endemic, newborns must introduce an endless supply o new susceptibles into a population.
59
4.6. EVOLUTION OF VIRULENCE
4.6
Evolution o virulence
Microorganisms continuously evolve due to selection pressures in their environments. Antibiotics are a common source o selection pressure on pathogenic bacteria, and the development o antibiotic-resistant strains presents a major health challenge to medical science. Bacteria and viruses also compete directly with each other or reproductive success resulting in the evolution o virulence. Here, using the SIR endemic disease model, we study how virulence may evolve. For the sake o argument, we will assume that a population is initially in equilibrium with an endemic disease caused by a wildtype virus; that is, , and are assumed to be nonzero and at equilibrium values. Now suppose that some virus particles mutate by a random, undirected process that occurs naturally. We want to determine the conditions under which the mutant virus will replace the wildtype virus in the population. In mathematical terms, we want to determine the linear stability o the endemic disease equilibrium to the introduction o a mutant viral strain. We assume that the original wildtype virus has inection rate , removal rate , and disease-related death rate , and that the mutant virus has corresponding rates ′ , ′ and ′ . We urther assume that an individual inected with either a wildtype or mutant virus gains immunity to subsequent inection rom both wildtype and mutant viral orms. Our model thus has a single susceptible class , two distinct inective classes and ′ depending on which virus causes the inection, and a single recovered class . The appropriate diagram is
6 ( + )
-
- - ′
?
′ ′
′ ′
( + ′ ) ′
?
?
with corresponding diferential equations = ( + ′ ′ ), = ( + + ), ′ = ′ ′ ( ′ + + ′ ) ′ , = + ′ ′ .
− − − − −
(4.12) (4.13) (4.14) (4.15)
I the population is initially in equilibrium with the wildtype virus, then we have ˙ = 0 with = 0, and the equilibrium value or is determined rom (4.13) to be + + * = , (4.16)
̸
60
CHAPTER 4. INFECTIOUS DISEASE MODELING
which corresponds to a basic reproductive ratio * /( + + ) o unity. We perturb this endemic disease equilibrium by introducing a small number o inectives carrying the mutated virus, that is, by letting ′ be small. Rather than solve the stability problem by means o a Jacobian analysis, we can directly examine the equation or ′ / given by (4.14). Here, with = * given by (4.16), we have ′ ′ ( + + ) =
′
′
− ( + + )
′ ;
and ′ increases exponentially i ′ ( + + )
′
′
− ( + + ) > 0,
or ater some elementary algebra, ′ > . ′ + + ′ + +
(4.17)
Our result (4.17) suggests that endemic viruses (or other microorganisms) will tend to evolve (i) to be more easily transmitted between people ( ′ > ); (ii) to make people sick longer ( ′ < ), and; (iii) to be less deadly ′ < . In other words, viruses evolve to increase their basic reproductive ratios. For instance, our model suggests that viruses evolve to be less deadly because the dead do not spread disease. Our result would not be applicable, however, i the dead in act did spread disease, a possibility i disposal o the dead was not done with sucient care, perhaps because o certain cultural traditions such as amily washing o the dead body.
Chapter 5
Population Genetics Deoxyribonucleic acid, or DNA—a large double-stranded, helical molecule, with rungs made rom the our base pairs adenine (A), cytosine (C), thymine (T) and guanine (G)—carries inherited genetic inormation. The ordering o the base pairs A, C, T and G determines the DNA sequence. A gene is a particular DNA sequence that is the undamental unit o heredity or a particular trait. Some species develop as diploids, carrying two copies o every gene, one rom each parent, and some species develop as haploids with only one copy. There are even species that develop as both diploids and haploids. When we say there is a gene or pea color, say, we mean there is a particular DNA sequence that may vary in a pea plant population, and that there are at least two subtypes, called alleles , where plants with two copies o the yellow-color allele have yellow peas, those with two copies o the green-color allele, green peas. A plant with two copies o the same allele is homozygous or that particular gene (or a homozygote ), while a plant carrying two diferent alleles is heterozygous (or a heterozygote ) . For the pea color gene, a plant carrying both a yellow- and green-color allele has yellow peas. We say that the green color is a recessive trait (or the green-color allele is recessive ), and the yellow color is a dominant trait (or the yellow-color allele is dominant ). The combination o alleles carried by the plant is called its genotype , while the actual trait (green or yellow peas) is called its phenotype . A gene that has more than one allele in a population is called polymorphic , and we say the population has a polymorphism or that particular gene. Population genetics can be dened as the mathematical modeling o the evolution and maintenance o polymorphism in a population. Population genetics together with Charles Darwin’s theory o evolution by natural selection and Gregor Mendel’s theory o biological inheritance orms the modern evolutionary synthesis (sometimes called the modern synthesis, the evolutionary synthesis, the neo-Darwinian synthesis, or neo-Darwinism). The primary ounders in the early twentieth century o population genetics were Sewall Wright, J. B. S. Haldane and Ronald Fisher. Allele requencies in a population can change due to the inuence o our primary evolutionary orces: natural selection, genetic drit, mutation, and migration. Here, we mainly ocus on natural selection and mutation. Genetic drit is the study o stochastic efects, and it is important in small populations. Migration requires consideration o the spatial distribution o a population, and it is usually modeled mathematically by partial diferential equations. 61
62
CHAPTER 5. POPULATION GENETICS
genotype number viability tness ertility tness
Table 5.1: Haploid genetics using population size, absolute viability, and ertility tnesses.
The simplied models we will consider assume innite population sizes (neglecting stochastic efects except in 5.5), well-mixed populations (neglecting any spatial distribution), and discrete generations (neglecting any age-structure). Our main purpose is to illustrate the undamental ways that a genetic polymorphism can be maintained in a population.
S
5.1
Haploid genetics
We rst consider the modeling o selection in a population o haploid organisms. Selection is modeled by tness coecients, with diferent genotypes having diferent tnesses. We begin with a simple model that counts the number o individuals in the next generation, and then show how this model can be reormulated in terms o allele requencies and relative tness coecients. Table 5.1 ormulates the basic model. We assume that there are two alleles and or a particular haploid gene. These alleles are carried in the population by and individuals, respectively. A raction ( ) o individuals carrying allele () is assumed to survive to reproduction age, and those that survive contribute ( ) ofspring to the next generation. These are o course average values, but under the assumption o an (almost) innite population, our model is deterministic. Accordingly, with () (( ) ) representing the number o individuals carrying allele () in the th generation, and ormulating a discrete generation model, we have
(+1) = () ,
(+1) = () .
(5.1)
It is mathematically easier and more transparent to work with allele requencies rather than individual numbers. We denote the requency (or more accurately, proportion) o allele () in the th generation by ( ); that is, ()
=
, () + ( )
()
=
() + ( )
,
where evidently + = 1. Now, rom (5.1), (+1)
()
+ (+1) = + () ,
(5.2)
63
5.1. HAPLOID GENETICS
genotype req. o gamete relative tness req ater selection normalization
1+ 1 (1 + ) / / = (1 + ) +
Table 5.2: Haploid genetic model o the spread o a avored allele. so that dividing the rst equation in (5.1) by (5.2) yields
()
+1 =
() + ( ) = + =
(5.3)
,
+
where the second equality comes rom dividing the numerator and denominator () () by + , and the third equality rom dividing the numerator and denominator by . Similarly, +1 =
,
(5.4)
+
which could also be derived using +1 = 1 +1 . We observe rom the evolution equations or the allele requencies, (5.3) and (5.4), that only the relative tness / o the alleles matters. Accordingly, in our models, we will consider only relative tnesses, and we will arbitrarily set one tness to unity to simpliy the algebra and make the nal result more transparent.
−
5.1.1
Spread o a avored allele
We consider a simple model or the spread o a avored allele in Table 5.2, with > 0. Denoting ′ by the requency o in the next generation (not (!) the derivative o ), the evolution equation is given by (1 + ) (1 + ) = , 1 +
′ =
(5.5)
where we have used (1 + ) + = 1 + , since + = 1. Note that (5.5) is the same as (5.3) with ′ = +1 , = , and / = 1 + . Fixed points o (5.5) are determined rom ′ = . We nd two xed points: * = 0, corresponding to a population in which allele is absent; and * = 1, corresponding to a population in which allele is xed . Intuitively, * = 0 is unstable while * = 1 is stable.
64
CHAPTER CHAPTER 5. POPULA POPULATION TION GENETIC GENETICS S
To illustrate how a stability analysis is perormed analytically or a diference equation (instead o a diferential equation), consider the general diference equation ′ = ( ( ) ). (5.6) With = * a xed point such that * = ( ( * ), we write = * + so that (5.6) 5.6) becomes * + ′ = ( ( * + ) = ( ( * ) + ′ ( * ) + . . . = * + ′ ( * ) + . . . , where ′ ( * ) denotes the derivative o evaluated at *. Thereore, to leadingorder in ′ / = ′ ( * ) ,
| | |
|
and the xed point is stable provided that ′ ( * ) < 1. For our haploid model,
|
( ( ) ) =
(1 + ) , 1 +
|
′ ( ) ) =
1+ , (1 + ) )2
so that ′ ( * = 0) = 1 + > 1, and ′ ( * = 1) = 1/ 1/(1 + ) < 1, conrming that * = 0 is unstable and * = 1 is stable. I the selection coecient is small, the model model equati equation on (5.5 ( 5.5)) simplies urther. We have (1 + ) 1 + = (1 + ) (1 (1
′ =
= + ( (
O(2 )) − + O(
− 2) + O( O(2 ),
so that to leading-order in , ′
(1 − ) ). − = (1 − ≪ 1, which is valid or ≪ 1, we can approximate this diference
I ′ equation by the diferential equation
/ = (1 (1
− ) ),
which shows that the requency o allele satises the now very amiliar logistic equation. Although a polymorphism or this gene exists in the population as the new allele spreads, eventually becomes xed in the population and the polymorphism phism is lost. lost. In the next section, section, we consider consider how a polymorph polymorphism ism can be maintained in a haploid population by a balance between mutation and selection.
5.1.2
MutationMutation-sele selectio ction n balance balance
We consider a gene with two alleles: alleles: a wildtype allele and a mutant allele . We view view the mutant mutant allele allele as a deect deectiv ivee genot genotype, ype, which which coner conerss on the
65
5.1. HAPLOID HAPLOID GENETICS GENETICS
genotype req. o gamete relative tness req ater selection req req ate aterr muta mutati tion on normalization
1 1 / (1 )/ (1-u) (1-u)p/ p/w w [(1 [(1 ) + + ] ]/ = + (1 )
−
−
−
−
Table 5.3: A haploid haploid genetic genetic model o mutation-s mutation-selec election tion balance. carrier a lowered tness 1 relative relative to the wildtype. wildtype. Although Although all mutant mutant alleles may not have identical DNA sequences, we assume that they share in common commo n the same phenotype o reduced tness. We model the opposing efects o two evolutionary orces: natural selection, which avors the wildtype allele over the mutant allele , and mutation, which coners a small probability that allele mutates mutates to allele allele in each newborn individual. Schematically, Schematically,
−
−−
⇀ ↽ where represents mutation and represents selection. The model is shown in Table 5.3. The equations or and in the next generation are
− ) (1 − ) = , 1 − (1 − ) )
′ =
and
(1
(5.7)
− ) + + (1 − − ) + + = , 1 −
′ =
(1
(5.8)
where we have used + = 1 to eliminate rom the equation or ′ and rom the equation or ′ . The equat equation ionss or ′ and ′ are linearly dependent since ′ ′ + = 1, and we need solve only one o them. Considering (5.7 (5.7), ), the xed points determined rom ′ = are * = 0, or which the mutant allele is xed in the population and there is no polymorphism, and the solution to 1
− (1 − ) = 1 − , *
which is * = 1 /, /, and there there is a polymor polymorphi phism. sm. The stabili stability ty o these these ′ two xed points is determined by considering = ( ( ), ), with ( ( ) ) given by the right-hand-side o (5.7 (5.7). ). Taking the derivative o , ,
−
′ ( ) ) = so that
(1 [1
− )(1 − ) , )]2 − (1 − )]
1 1 , ′ ( * = 1 /) /) = . 1 1 Applying the criterion ′ ( * ) < 1 or stability, * = 0 is stable or < and * = 1 / is stable or > . A polymorphism polymorphism is thereore possible possible under mutation-selection balance when > > 0. ′ ( * = 0) =
−
|
− − |
−
− −
66
CHAPTER CHAPTER 5. POPULA POPULATION TION GENETIC GENETICS S
genotype reerr reerred ed to as requency
wildt wildtype ype homozy homozygot gotee
hetero heterozyg zygote ote
mutan mutantt homozy homozygot gotee
Table 5.4: The terminology o diploidy.
5.2 5.2
Dipl Diploi oid d gene geneti tics cs
Most sexually reproducing species are diploid. In particular, our species Homo diploid with with two two except exception ions: s: we are haploid haploid at the gamete gamete stage sapiens is diploid (sperm and unertilized egg); and males are haploid or most genes on the unmatched matched X and Y sex chromoso chromosomes mes (emales (emales are XX and diploid). This latter seemingly innocent act is o great signicance to males sufering rom genetic diseases due to an X-linked recessive mutation inherited rom their mother. Females inheriting this mutation are most probably disease-ree because o the unctional unctional gene inherited inherited rom their ather. A polymorphic gene with alleles and can appear in a diploid gene as three distinct distinct genotypes: genotypes: ,, and . . Conventionally Conventionally,, we denote to be the wildtype allele and the mutant allele. Table 5.4 presents the terminology o diploidy. As or haploi haploid d geneti genetics, cs, we will will determ determine ine evolu evolutio tion n equati equations ons or allele allele and/or and/or genotype genotype requencies. requencies. To develop the appropriate appropriate denitions and relations, we initially assume a population o size (which we will later take to be innite), innite), and assume assume that the number o individuals with genotypes genotypes ,, and are , and . Now, = + + . Dene genotype requencies , , and as =
,
=
,
=
,
so that + + + = 1. It will also be useul to dene allele requencies. Let and be the number o alleles and in the population, with = + the total number o alleles. Since the population is o size and diploidy, = 2 ; and since each homozygote contains two identical alleles, and each heterozygote contains one o each allele, = 2 2 + and = 2 + . Dening the allele requencies and as previously, = / 2 + = 2 1 = + ; 2 and similarly, = / 2 + = 2 1 = + . 2
67
5.2. DIPLOID GENETICS
With ve requencies, , , , , , and our constraints + + = 1, + = 1, = + /2, = + /2, how many independent requencies are there? In act, there are two because one o the our constraints is linearly dependent. We may choose any two requencies other than the choice , as our linearly independent set. For instance, one choice is , ; then, = 1
− ,
= 2(
{ }
{ } = 1 + − 2 .
− ),
Similarly, another choice is , ; then
{
=1
5.2.1
− − ,
}
1 = + , 2
= 1
− − 12 .
Sexual reproduction
Diploid reproduction may be sexual or asexual, and sexual reproduction may be o varying types (e.g., random mating, selng, brother-sister mating, and various other types o assortative mating). The two simplest types to model exactly are random mating and selng. These mating systems are useul or contrasting the biology o both outbreeding and inbreeding. Random mating The simplest type o sexual reproduction to model mathematically is random mating . Here, we assume a well-mixed population o individuals that have equal probability o mating with every other individual. We will determine the genotype requencies o the zygotes (ertilized eggs) in terms o the allele requencies using two approaches: (1) the gene pool approach, and (2) the mating table approach. The gene pool approach models sexual reproduction by assuming that males and emales release their gametes into pools. Ofspring genotypes are determined by randomly combining one gamete rom the male pool and one gamete rom the emale pool. As the probability o a random gamete containing allele or is equal to the allele’s population requency or , respectively, the probability o an ofspring being is 2 , o being is 2 (male emale + emale male ), and o being is 2 . Thereore, ater a single generation o random mating, the genotype requencies can be given in terms o the allele requencies by = 2 , = 2 , = 2 . This is the celebrated Hardy-Weinberg law. Notice that under the assumption o random mating, there is now only a single independent requency, greatly simpliying the mathematical modeling. For example, i is taken as the independent requency, then = 1
− ,
= 2 ,
= 2 (1
− ),
= (1
− )2.
Most modeling is done assuming random mating unless the biology under study is inuenced by inbreeding. The second approach uses a mating table (see Table 5.5). This approach to modeling sexual reproduction is more general and can be applied to other mating systems. We explain this approach by considering the mating . The
×
68
CHAPTER 5. POPULATION GENETICS
mating Totals
× × × × × ×
requency 2 2 2 2 2 2 ( + + )2 =1
AA 2 0 1 2 4 0 0 ( + 12 )2 = 2
progeny requency Aa 0 2 1 2 2 0 2( + 12 )( + 12 ) = 2
aa 0 0 0 1 2 4 2 ( + 12 )2 = 2
Table 5.5: Random mating table.
mating Totals
⊗ ⊗ ⊗
requency 1
progeny requency AA Aa aa 0 0 1 1 1 4 2 4 0 0 1 1 + 4 2 + 14
Table 5.6: Selng mating table. genotypes and have requencies and , respectively. The requency o males mating with emales is and is the same as emales mating with males, so the sum is 2 . Hal o the ofspring will be and hal , and the requencies are denoted under progeny requency . The sums o all the progeny requencies are given in the Totals row, and the random mating results are recovered upon use o the relationship between the genotype and allele requencies. Selfng Perhaps the next simplest type o mating system is sel-ertilization, or selng. Here, an individual reproduces sexually (passing through a haploid gamete stage in its lie-cycle), but provides both o the gametes. For example, the nematode worm C. elegans can reproduce by selng. The mating table or selng is given in Table 5.6. The selng requency o a particular genotype is just the requency o the genotype itsel. For a selng population, disregarding selection or any other evolutionary orces, the genotype requencies evolve as 1 ′ = + , 4
′ =
1 , 2
1 ′ = + . 4
(5.9)
Assuming an initially heterozygous population, we solve (5.9) with the initial conditions 0 = 1 and 0 = 0 = 0. In the worm lab, this type o initial population is commonly created by crossing wildtype homozygous C. elegans males with mutant homozygous C. elegans hermaphrodites, where the mutant allele is recessive. Wildtype hermaphrodite ofspring, which are necessarily heterozygous, are then picked to separate worm plates and allowed to sel-ertilize. (Do you see why the experiment is not done with wildtype hermaphrodites and
69
5.2. DIPLOID GENETICS
mutant males?) From the equation or ′ in (5.9), we have = (1/2) , and rom symmetry, = . Then, since + + = 1, we obtain the complete solution 1 = 2
− 1
1 2
,
=
1 2
,
1 = 2
− 1
1 2
.
The main result to be emphasized here is that the heterozygosity o the population decreases by a actor o two in each generation. Selng populations rapidly become homozygous. For homework, I will ask you to determine the genotype requencies or an initially random mating population that transitions to selng. Constancy o allele requencies We can show that both random mating and selng do not by themselves change the allele requencies o a population, but only reshues alleles into diferent genotypes. For random mating, 1 ′ = ′ + ′ 2 1 = 2 + (2 ) 2 = ( + ) = ; and or selng, 1 ′ = ′ + ′ 2 1 1 = + + 4 2 1 = + 2 = .
1 2
These results are not general, however, and it is possible to construct other mating systems or which allele requencies do change. Nevertheless, the conservation o allele requencies by random mating is an important element o neo-Darwinism. In Darwin’s time, most biologists believed in blending inheritance , where the genetic material rom parents with diferent traits actually blended in their ofspring, rather like the mixing o paints o diferent colors. I blending inheritance occurred, then genetic variation, or polymorphism, would eventually be lost over several generations as the “genetic paints” became well-mixed. Mendel’s work on peas, published in 1866, suggested a particulate theory o inheritance, where the genetic material, later called genes, maintain their integrity across generations. Sadly, Mendel’s paper was not read by Darwin (who published The Origin o Species in 1859 and died in 1882) or other inuential biologists during Mendel’s lietime (Mendel died in 1884). Ater being rediscovered in 1900, Mendel and his work eventually became widely celebrated.
70
CHAPTER 5. POPULATION GENETICS
Figure 5.1: Evolution o peppered moths in industrializing England.
5.2.2
Spread o a avored allele
We consider the spread o a avored allele in a diploid population. The classic example – widely repeated in biology textbooks as a modern example o natural selection – is the change in the requencies o the dark and light phenotypes o the peppered moth during England’s industrial revolution. The evolutionary story begins with the observation that pollution killed the light colored lichen on trees during industrialization o the cities. Light colored peppered moths camouage well on light colored lichens, but are exposed to birds on plain tree bark. On the other hand, dark colored peppered moths camouage well on plain tree bark, but are exposed on light colored lichens (see Fig. 5.1). Natural selection thereore avored the light-colored allele in preindustrialized England and the dark-colored allele during industrialization. The dark-colored allele, presumably kept at low requency by mutation-selection balance in pre-industrialized England, increased rapidly under natural selection in industrializing England. We present our model in Table 5.7. Here, we consider as the wildtype genotype and normalize its tness to unity. The allele is the mutant whose requency increases in the population. In our example o the peppered moth, the phenotype is light colored and the phenotype is dark colored. The color o the phenotype depends on the relative dominance o and . Usually, no pigment results in light color and is a consequence o nonunctioning pigmentproducing genes. One unctioning pigment-producing allele is usually sucient to result in a dark-colored moth. With a unctioning pigment-producing allele and the mutated nonunctioning allele, is most likely recessive, is most likely dominant, and the phenotype o is most likely dark, so ℎ 1. For the moment, though, we leave ℎ as a ree parameter.
≈
We assume random mating, and this simplication is used to write the genotype requencies as = 2 , = 2 , and = 2 . Since = 1 , we reduce our problem to determining an equation or ′ in terms o . Using ′ = +(1/2) , where ′ is the allele requency in the next generation’s zygotes, and and are the and genotype requencies, respectively, in the present generation
−
71
5.2. DIPLOID GENETICS
genotype req. o zygote relative tness req ater selection normalization
2 2 2 1+ 1 + ℎ 1 2 2 (1 + ) / 2(1 + ℎ) / / = (1 + ) 2 + 2(1 + ℎ) + 2
Table 5.7: A diploid genetic model o the spread o a avored allele assuming random mating. ater selection, ′ = where = 1
(1 + ) 2 + (1 + ℎ) ,
− , and = (1 + ) 2 + 2(1 + ℎ) + 2 = 1 + ( 2 + 2ℎ ).
Ater some algebra, the nal evolution equation written solely in terms o is ′ =
(1 + ℎ) + (1 ℎ) 2 . 1 + 2ℎ + (1 2ℎ) 2
− −
(5.10)
The expected xed points o this equation are * = 0 (unstable) and * = 1 (stable), where our assignment o stability assumes positive selection coecients. The evolution equation (5.10) in this orm is not particularly illuminating. In general, a numerical solution would require speciying numerical values or and ℎ, as well as an initial value or . Here, to determine how the spread o depends on the dominance coecient ℎ, we investigate analytically the increase o assuming 1. We Taylor-series expand the right-hand-side o (5.10) in powers o , keeping terms to order :
≪
(1 + ℎ) + (1 ℎ) 2 1 + 2ℎ + (1 2ℎ) 2 + ℎ + (1 ℎ) 2 = 1 + 2ℎ + (1 2ℎ) 2
− − − − = + ℎ + (1 − ℎ) 2 1 − 2ℎ + (1 − 2ℎ) 2 = + ℎ + (1 − 3ℎ) − (1 − 2ℎ) 2 + O(2 ).
′ =
(( ( (
)) ) (
)
(5.11)
)
+ O(2 )
I 1, we expect a small change in allele requency in each generation, so we can approximate ′ /, where denotes the generation number, and = (). The approximate diferential equation obtained rom (5.11) is
≪
− ≈
= ℎ + (1
− 3ℎ) − (1 − 2ℎ) 2 . (5.12) I is partially dominant so that ℎ = ̸ 0 (e.g., the heterozygous moth is
(
)
darker than the homozgygous mutant moth), then the solution to (5.12) behaves similarly to the solution o a logistic equation: initially grows exponentially as () = 0 exp(ℎ), and asymptotes to one or large . I is recessive so
72
CHAPTER 5. POPULATION GENETICS
disease Thalassemia Sickle cell anemia Haemophilia Cystic Fibrosis Tay-Sachs disease Fragile X syndrome Huntington’s disease
mutation haemoglobin haemoglobin blood clotting actor chloride ion channel Hexosaminidase A enzyme FMR1 gene HD gene
symptoms anemia anemia uncontrolled bleeding thick lung mucous nerve cell damage mental retardation brain degeneration
Table 5.8: Seven common monogenic diseases. that ℎ = 0 (e.g., the heterozygous moth is as light-colored as the homozygous mutant moth), then (5.12) reduces to = 2 (1
− ) ,
or ℎ = 0.
(5.13)
O main interest is the initial growth o when (0) = 0 1, so that / 2 . This diferential equation may be integrated by separating variables to yield
≪
() =
≈
≈
0 1 0 0 (1 + 0 ).
−
The requency o a recessive avored allele increases only linearly across generations, a consequence o the heterozygote being hidden rom natural selection. Most likely, the peppered-moth heterozygote is signicantly darker than the light-colored homozygote since the dark colored moth rapidly increased in requency over a short period o time. As a nal comment, linear growth in the requency o when ℎ = 0 is sensitive to our assumption o random mating. I selng occurred, or another type o close amily mating, then a recessive avored allele may still increase exponentially. In this circumstance, the production o homozygous ofspring rom more requent heterozygote pairings allows selection to act more efectively.
5.2.3
Mutation-selection balance
By virtue o sel-knowledge, the species with the most known mutant phenotypes is Homo sapiens . There are thousands o known genetic diseases in humans, many o them caused by mutation o a single gene (called a monogenic disease). For an easy-to-read overview o genetic disease in humans, see the website http://www.who.int/genomics/public/geneticdiseases.
Table 5.8 lists seven common monogenic diseases. The rst two diseases are maintained at signicant requencies in some human populations by heterosis . We will discuss in 5.2.4 the maintenance o a polymorphism by heterosis, or which the heterozygote has higher tness than either homozygote. It is postulated that Tay-Sachs disease, prevalent among ancestors o Eastern European Jews, and cystic brosis may also have been maintained by heterosis acting in
S
73
5.2. DIPLOID GENETICS
genotype req. o zygote relative tness
2 2 2 1 1 ℎ 1 2 / 2(1 ℎ) / (1 ) 2 / = 2 + 2(1 ℎ) + (1 ) 2
− − −
freq after selection
normalization
− −
−
Table 5.9: A diploid genetic model o mutation-selection balance assuming random mating. the past. (Note that the cystic brosis gene was identied in 1989 by a Toronto group led by Lap Chee Tsui, who later became President o the University o Hong Kong.) The other disease genes listed may be maintained by mutationselection balance. Our model or diploid mutation selection-balance is given in Table 5.9. We urther assume that mutations o type occur in gamete production with requency . Back-mutation is neglected. The gametic requency o and ater selection but beore mutation is given by ^ = + /2 and ^ = + /2, and the gametic requency o ater mutation is given by ′ = ^ + ^ . Thereore,
→
′
((
2
= + (1 where
) (
− ℎ) +
(1
− ) + (1 − ℎ)
= 2 + 2(1 ℎ) + (1 = 1 (2ℎ + ).
−
−
)
2
/,
− ) 2
Using = 1 , we write the evolution equation or ′ in terms o alone. Ater some algebra that could be acilitated using a computer algebra sotware such as Mathematica, we obtain
−
− ℎ(1 + ) − 1 − ℎ(1 + ) 2 = . (5.14) 1 − 2ℎ − (1 − 2ℎ) 2 To determine the equilibrium solutions o (5.14), we set ≡ = to obtain ′
(−
+ 1
) (
)
*
′
a cubic equation or * . Because o the neglect o back mutation in our model, one solution readily ound is * = 1, in which all the alleles have mutated to . The * = 1 solution may be actored out o the cubic equation resulting in a quadratic equation, with two solutions. Rather than show the exact result here, we determine equilibrium solutions under two approximations: (i) 0 < ℎ, , and; (ii) 0 = ℎ < < . First, when 0 < ℎ, , we look or a solution o the orm * = +O(2 ), with constant, and Taylor series expand in (assuming , ℎ = O(0 )). I such a solution exists, then (5.14) will determine the unknown coecient . We have
≪
≪
+ O(2 ) =
+ (1 ℎ) + O(2 ) 1 2ℎ + O(2 )
−
−
− ℎ) + O(2); and equating powers o , we nd = 1 + − ℎ, or = 1/ℎ. Thereore, = /ℎ + O(2 ), or 0 < ≪ ℎ,. = (1 +
*
74
CHAPTER 5. POPULATION GENETICS
genotype req: 0 < , ℎ req: 0 = ℎ < <
≪
1 + O() 1 + O( )
2/ℎ + O(2 ) 2 / + O()
√
√
+ O(3 ) /
2 /(ℎ)2
Table 5.10: Equilibrium requencies o the genotypes at the diploid mutationselection balance. Second, when 0 = ℎ < < , we substitute ℎ = 0 directly into (5.14), * =
+ (1 ) * 1 *2
− −
− 2 , *
which we then write as a cubic equation * , *3
− 2 − + = 0. *
*
By actoring this cubic equation, we nd ( *
− 1)( 2 − /) = 0; *
and the polymorphic equilibrium solution is * =
√
/,
or 0 = ℎ < < .
Because * < 1 only i > , this solution does not exist i < . Table 5.10 summarizes our results or the equilibrium requencies o the genotypes at mutation-selection balance. The rst row o requencies, 0 < , ℎ, corresponds to a dominant (ℎ = 1) or partially-dominant ( ℎ < 1) mutation, where the heterozygote is o reduced tness and shows symptoms o the genetic disease. The second row o requencies, 0 = ℎ < < , corresponds to a recessive mutation, where the heterozygote is symptom-ree. Notice that individuals carrying a dominant mutation are twice as prevalent in the population as individuals homozygous or a recessive mutation (with the same and ). A heterozygote carrying a dominant mutation most commonly arises either de novo (by direct mutation o allele ) or by the mating o a heterozygote with a wildtype. The latter is more common or 1, while the ormer must occur or = 1 (a heterozygote with an = ℎ = 1 mutation by denition does not reproduce). One o the most common autosomal dominant genetic diseases is Huntington’s disease, resulting in brain deterioration during middle age. Because individuals with Huntington’s disease have children beore disease symptoms appear, is small and the disease is usually passed to ofspring by the mating o a (heterozygote) with a wildtype homozygote. For a recessive mutation, a mutant homozygote usually occurs by the mating o two heterozygotes. I both parents carry a single recessive disease allele, then their child has a 1/4 chance o getting the disease.
≪
≪
≪
5.2.4
Heterosis
Heterosis , also called overdominance or heterozygote advantage , occurs when the heterozygote has higher tness than either homozygote. The best-known examples are sickle-cell anemia and thalassemia, diseases that both afect hemoglobin,
75
5.3. FREQUENCY-DEPENDENT SELECTION
genotype req. o zygote relative tness req ater selection normalization
2 2 2 1 1 1 2 (1 ) / 2 / (1 ) 2 / = (1 ) 2 + 2 + (1 ) 2
−
−
− −
−
−
Table 5.11: A diploid genetic model o heterosis assuming random mating. the oxygen-carrier protein o red blood cells. The sickle-cell mutations are most common in people o West Arican descent, while the thalassemia mutations are most common in people rom the Mediterranean and Asia. In Hong Kong, the television stations occasionally play public service announcements concerning thalassemia. The heterozygote carrier o the sickle-cell or thalassemia gene is healthy and resistant to malaria; the wildtype homozygote is healthy, but susceptible to malaria; the mutant homozygote is sick with anemia. In class, we will watch the short video, A Mutation Story , about the sickle cell gene, which you can also view at your convenience: http://www.pbs.org/wgbh/evolution/library/01/2/l_012_02.html.
Table 5.11 presents our model o heterosis. Both homozygotes are o lower tness than the heterozygote, whose relative tness we arbitrarily set to unity. Writing the equation or ′ , we have
− ) 2 + − 2 − 2 − 2 = . 1 − + 2 − ( + ) 2
′ =
At equilibrium, *
′
≡
(1 1
= , and we obtain a cubic equation or * : ( + ) 3*
− ( + 2) 2 + *
*
= 0.
(5.15)
Evidently, * = 0 and * = 1 are xed points, and (5.15) can be actored as (1
− ) ( − ( + ) ) = 0.
The p olymorphic solution is thereore * =
, +
* =
, +
valid when , > 0. Since the value o * can be large, recessive mutations that cause disease, yet are highly prevalent in a population, are suspected to provide some benet to the heterozygote. However, only a ew genes are unequivocally known to exhibit heterosis.
5.3
Frequency-dependent selection
A polymorphism may also result rom requency-dependent selection. A wellknown model o requency-dependent selection is the Hawk-Dove game. Most
76
CHAPTER 5. POPULATION GENETICS
player H D
∖ opponent
H = 2 = 0
−
D = 2 = 1
Table 5.12: General payof matrix or the Hawk-Dove game, and the usually assumed values. The payofs are payed to the player (rst column) when playing against the opponent (rst row). genotype req o zygote relative tness req ater selection normalization
H p
D q
+ + + + ( + + ) / ( + + )/ = ( + + ) + ( + + )
Table 5.13: Haploid genetic model or requency-selection in the Hawk-Dove game commonly, requency-dependent selection is studied using game theory, and ollowing John Maynard Smith, one looks or an evolutionarily stable strategy (ESS). Here, we study requency-dependent selection using a populationgenetics model. We consider two phenotypes: Hawk and Dove, with no mating between diferent phenotypes (or example, diferent phenotypes may correspond to dierent species, such as hawks and doves). We describe the Hawk-Dove game as ollows: (i) when Hawk meets Dove, Hawk gets the resource and Dove retreats beore injury; (ii) when two Hawks meet, they engage in an escalating ght, seriously risking injury, and; (iii) when two Doves meet, they share the resource. The Hawk-Dove game is modeled by a payof matrix, as shown in Table 5.12. The player in the rst column receives the payof when playing the opponent in the rst row. For instance, Hawk playing Dove gets the payof . The numerical values are commonly chosen such that < < < , that is, Hawk playing Dove does better than Dove playing Dove does better than Dove playing Hawk does better than Hawk playing Hawk. Frequency-dependent selection occurs because the tness o Hawk or Dove depends on the requency o Hawks and Doves in the population. For example, a Hawk in a population o Doves does well, but a Hawk in a population o Hawks does poorly. We model the Hawk-Dove game using a haploid genetic model, with and the population requencies o Hawks and Doves, respectively. We assume that the times a player phenotype plays an opponent phenotype is proportional to the population requency o the opponent phenotype. For example, i the population is composed o 1/4 Hawks and 3/4 Doves, then Hawk or Dove plays Hawk 1/4 o the time and Dove 3/4 o the time. Our haploid genetic model assuming requency-dependent selection is ormulated in Table 5.13. The tness parameter is arbitrary, but assumed to be large and positive so that the tness o Hawk or Dove is always positive. A negative tness in our haploid model has no meaning (see 5.1). From Table 5.13, the equation or ′ is
S
′ = ( + + ) /,
(5.16)
77
5.3. FREQUENCY-DEPENDENT SELECTION
with = ( + + ) + ( + + ).
(5.17)
Both * = 0 and * = 1 are xed points; a stable polymorphism exists when both these xed points are unstable. First, consider the xed point * = 0. With perturbation 1 and = 1 , to leading-order in
≪
−
′
=
+ . +
Thereore, * = 0 is an unstable xed point when > —a population o Doves is unstable i Hawk playing Dove does better than Dove playing Dove. By symmetry, * = 1 is an unstable xed point when > —a population o Hawks is unstable i Dove playing Hawk does better than Hawk playing Hawk. Both homogeneous populations o either all Hawks or all Doves are unstable or the numerical values shown in Table 5.12. The polymorphic solution may be determined rom (5.16) and (5.17). Assuming * ′ = , * = 1 * , and * , * = 0, we have
≡
−
̸
( + * + * ) * + ( + * + * ) * = + * + * ; and solving or * , we have * =
(
−
) + (
−
− ) .
Since we have assumed > and > , the polymorphic solution satises 0 < * < 1. With the numerical values in Table 5.12, * =
(2 = 1/3;
−
2 1 1) + (0 + 2)
−
the stable polymorphic population, maintained by requency-dependent selection, consists o 1/3 Hawks and 2/3 Doves. Interestingly, this stable polymorphism is not the best outcome or the population as a whole. The normalization actor can be interpreted to be proportional to the average tness o the entire population (the relative tness o doves times the requency o doves plus the relative tness o hawks times the requency o hawks). At the stable polymorphism o 1/3 Hawks and 2/3 Doves, the value o is calculated to be = + 1/3. However, the maximum possible value o —computed by determining the value o such that / = 0—can be shown to occur or a population consisting o all Doves, and this value is calculated to be = + 1. I only all the Hawks will act like Doves, the entire population would benet as a whole. Yet a single Hawk that acts like its true sel has a selective advantage, and natural selection results in an all Dove population becoming unstable. In this situation, natural selection acts to move the population as a whole away rom its maximum average tness by avoring the evolution o Hawks.
78
CHAPTER 5. POPULATION GENETICS
5.4
Recombination and the approach to linkage equilibrium
When considering a polymorphism at a single genetic locus, we assumed two distinct alleles, and . The diploid then occurs as one o three types: , and . We now consider a polymorphism at two genetic loci, each with two distinct alleles. I the alleles at the rst genetic loci are and , and those at the second and , then our distinct haploid gametes are possible, namely , , and . Ten distinct diplotypes are possible, obtained by orming pairs o all possible haplotypes. We can write these ten diplotypes as /, /, /, /, /, /, /, /, /, and /, where the numerator represents the haplotype rom one parent, the denominator represents the haplotype rom the other parent. We do not distinguish here which haplotype came rom which parent. To proceed urther, we dene the allelic and gametic requencies or our two loci problem in Table 5.14. I the probability that a gamete contains allele or does not depend on whether the gamete contains allele or , then the two loci are said to be independent. Under the assumption o independence, the gametic requencies are the products o the allelic requencies, i.e., = , = , etc. Oten, the two loci are not independent. This can be due to epistatic selection, or epistasis . As an example, suppose that two loci in humans inuence height, and that the most t genotype is the one resulting in an average height. Selection that avors the average population value o a trait is called normalizing or stabilizing. Suppose that and are hypothetical tall alleles, and are short alleles, and a person with two tall and two short alleles obtains average height. Then selection may avor the specic genotypes /, /, /, and /. Selection may act against both the genotypes yielding above average heights, /, /, and /, and those yielding below average heights, /, / and /. Epistatic selection occurs because the tness o the , loci depends on which alleles are present at the , loci. Here, has higher tness when paired with than when paired with . The two loci may also not be independent because o a nite population size (i.e., stochastic efects). For instance, suppose a mutation occurs only once in a nite population (in an innite population, any possible mutation occurs an innite number o times), and that is strongly avored by natural selection. The requency o may then increase. I a nearby polymorphic locus on the same chromosome as happens to be (say, with a polymorphism in the population), then gametes may substantially increase in requency, with absent. We say that the allele hitchhikes with the avored allele . When the two loci are not independent, we say that the loci are in gametic phase disequilibrium, or more commonly linkage disequilibrium , sometimes abbreviated as LD. When the loci are independent, we say they are in linkage
→
allele or gamete genotype requency
Table 5.14: Denitions o allelic and gametic requencies or two genetic loci each with two alleles.
5.4. LINKAGE EQUILIBRIUM
79
Figure 5.2: A schematic o crossing-over and recombination during meiosis (gure rom Access Excellence @ the National Health Museum)
equilibrium . Here, we will model how two loci, initially in linkage disequilibrium, approach linkage equilibrium through the process o recombination. To begin, we need a rudimentary understanding o meiosis . During meiosis, a diploid cell’s DNA, arranged in very long molecules called chromosomes, is replicated once and separated twice, producing our haploid cells, each containing hal o the original cell’s chromosomes. Sexual reproduction results in syngamy , the using o a haploid egg and sperm cell to orm a diploid zygote cell. Fig. 5.2 presents a schematic o meiosis and the process o crossing-over resulting in recombination. In a diploid, each chromosome has a corresponding sister chromosome, one chromosome originating rom the egg, one rom the sperm. These sibling chromosomes have the same genes, but possibly diferent alleles. In Fig. 5.2, we schematically show the alleles ,, on the light chromosome, and the alleles ,, on its sister’s dark chromosome. In the rst step o meiosis, each chromosome replicates itsel exactly. In the second step, sister chromosomes exchange genetic material by the process o crossing-over. All our chromosomes then separate into haploid cells. Notice rom the schematic that the process o crossing-over can result in genetic recombination. Suppose that the schematic o Fig. 5.2 represents the production o sperm by a male. I the chromosome rom the male’s ather contains the alleles and that rom the male’s mother , recombination can result in the sperm containing a chromosome with alleles (the third gamete in Fig. 5.2). We say this chromosome is a recombinant ; it contains alleles rom both its paternal grandather and paternal grandmother. It is likely that the precise combination o alleles on this recombinant chromosome has never existed beore in a single person. Recombination is the reason why everybody, with the exception o identical twins, is genetically unique.
80
CHAPTER 5. POPULATION GENETICS
Genes that occur on the same chromosome are said to be linked. The closer the genes are to each other on the chromosome, the tighter the linkage, and the less likely recombination will separate them. Tightly linked genes are likely to be inherited rom the same grandparent. Genes on diferent chromosomes are by denition unlinked; independent assortment o chromosomes results in a 50% chance o a gamete receiving either grandparents’ genes. To dene and model the evolution o linkage disequilibrium, we rst obtain allele requencies rom gametic requencies by = + ,
= + ,
= + ,
= + .
(5.18)
Since the requencies sum to unity, + = 1,
+ = 1,
+ + + = 1.
(5.19)
There are three independent gametic requencies and only two independent allelic requencies, so in general it is not possible to obtain the gametic requencies rom the allelic requencies without assuming an additional constraint such as linkage equilibrium. We can, however, introduce an additional variable , called the coecient o linkage disequilibrium, and dene to be the diference between the gametic requency and what this gametic requency would be i the loci were in linkage equilibrium: = + .
(5.20a)
=
− ;
(5.20b)
=
− ;
(5.20c)
Using + = , using + = , and using + = ,
= + .
(5.20d)
With our denition, positive linkage disequilibrium ( > 0) implies excessive and gametes and decient and gametes; negative linkage disequilibrium ( < 0) implies the opposite. An equality obtainable rom (5.20) that we will later nd useul is
− = ( + )( + ) − ( − )( − ) = ( + + + ) = .
(5.21)
Without selection and mutation, evolves only because o recombination. With primes representing the values in the next generation, and using ′ = and ′ = because sexual reproduction by itsel does not change allele requencies, ′ = ′ ′ ′ = ′ = ′ ( ) = + ( ′ ) ,
− − −
−
−
81
5.4. LINKAGE EQUILIBRIUM
diploid / / / / / / / / / /
dip req 2 2 2 2 2 2 2 2 2 2
gamete req / diploid req AB Ab aB 1 0 0 1/2 1/2 0 1/2 0 1/2 (1 )/2 /2 /2 (1 0 1 0 (1 )/2 (1 )/2 /2 0 1/2 0 0 0 1 0 0 1/2 0 0 0
−
−
−
ab 0 0 0 )/2 0 /2 1/2 0 1/2 1
−
Table 5.15: Computation o gamete requencies.
where we have used (5.20) to obtain the third equality. The change in is thereore equal to the change in requency o the gametes, ′
′
− = − .
(5.22)
To understand why gametic requencies change across generations, we should rst recognize when they do not change. Without genetic recombination, chromosomes maintain their exact identity across generations. Chromosome requencies without recombination are thereore constant, and or genetic loci with, say, alleles A,a and B,b on the same chromosome, ′ = . In an innite population without selection or mutation, gametic requencies change only or genetic loci in linkage disequilibrium on diferent chromosomes, or or genetic loci in linkage disequilibrium on the same chromosome subjected to genetic recombination. We will compute the requency ′ o gametes in the next generation, given the requency o gametes in the present generation, using two diferent methods. The rst method uses a ‘mating’ table (see Table 5.15). The rst column is the parent diploid genotype beore meiosis. The second column is the diploid genotype requency assuming random mating. The next our columns are the haploid genotype requencies (normalized by the corresponding diploid requencies to simpliy the table presentation). Here, we dene to be the requency at which the gamete arises rom a combination o grandmother and grandather genes. I the A,a and B,b loci occur on the same chromosome, then is the recombination requency due to crossing-over. I the A,a and B,b loci occur on diferent chromosomes, then because o the independent assortment o chromosomes there is an equal probability that the gamete contains all grandather or grandmother genes, or contains a combination o grandmother and grandather genes, so that = 1/2. Notice that crossing-over or independent assortment is o importance only i the grandather’s and grandmother’s contribution to the diploid genotype share no common alleles (i.e., / and / genotypes). The requency ′ in the next generation is given by the sum o the column (ater multiplication by the diploid requencies). There-
82
CHAPTER 5. POPULATION GENETICS
ore, ′ = 2 + + + (1
− ) + = ( + + + ) + ( − ) = − ,
(5.23)
where the nal equality makes use o (5.19) and (5.21). The second method or computing ′ is more direct. An haplotype can arise rom a diploid o general type / without recombination, or a diploid o type / with recombination. Thereore, ′ = (1
− ) + ,
where the rst term is rom non-recombinants and the second term rom recombinants. With = , we have
−
′ = (1 ) + ( = ,
−
−
− )
the same result as (5.23). Using (5.22) and (5.23), we derive ′ = (1
− ),
with the solution = 0 (1
− ).
Recombination decreases linkage disequilibrium in each generation by a actor o (1 ). Tightly linked genes on the same chromosome have small values o ; unlinked genes on diferent chromosomes have = 1/2. For unlinked genes, linkage disequilibrium decreases by a actor o two in each generation. We conclude that very strong selection is required to maintain linkage disequilibrium or genes on diferent chromosomes, while weak selection can maintain linkage disequilibrium or tightly linked genes.
−
5.5
Random genetic drit
Up to now, our simplied genetic models have all assumed an innite population, neglecting stochastic efects. Here, we consider a nite population: the resulting stochastic efects on the allelic requencies is called random genetic drit. The simplest genetic model incorporating random genetic drit assumes a xed-sized population o individuals, and models the evolution o a diallelic haploid genetic locus. There are two widely used genetic models or nite populations: the WrightFisher model and the Moran model. The Wright-Fisher model is most similar to our innite-population discrete-generation model. In the Wright-Fisher model, adult individuals release a very large number o gametes into a gene pool, and the next generation is ormed rom random gametes independently chosen rom the gene pool. The Moran model takes a diferent approach. In the Moran model, a single evolution step consists o one random individual in the population reproducing, and another random individual dying, with the population
83
5.5. RANDOM GENETIC DRIFT
always maintained at the constant size . Because two random events occur every generation in the Moran model, and random events occur every generation in the Wright-Fisher model, a total o /2 evolution steps in the Moran model is comparable, but not exactly identical, to a single discrete generation in the Wright-Fisher model. It has been shown, however, that the two models become identical in the limit o large . For our purposes, the Moran model is mathematically more tractable and we adopt it here. We develop our model analogously to the stochastic population growth model derived in 3.1. We let denote the number o -alleles in the population, and the number o -alleles. With = 0, 1, 2, . . . , a discrete random variable, the probability mass unction () denotes the probability o having -alleles at evolution step . With -alleles in a population o size , the probability o an individual with either reproducing or dying is given by = / ; the corresponding probability or is 1 . There are three ways to obtain a population o -alleles at evolution step + 1. First, there were -alleles at evolution step , and the individual reproducing carried the same allele as the individual dying. Second, there were 1 -alleles at evolution step , and the individual reproducing has and the individual dying has . And third, there were + 1 -alleles at evolution step , and the individual reproducing has and the individual dying has . Multiplying the probabilities and summing the three cases results in
−
S
−
−
( + 1) = 2 + (1
(
− )2 () +
)
−1
(1 −1 ) −1 () + +1 (1 +1 ) +1 ().
−
−
(5.24)
Note that this equation is valid or 0 < < , and that the equations at the boundaries—representing the probabilities that one o the alleles is xed—are 0 ( + 1) = 0 () + 1 (1 1 ) 1 (), ( + 1) = () + −1 (1 −1 ) −1 ().
−
(5.25a) (5.25b)
−
The boundaries are called absorbing and the probability o xation o an allele monotonically increases with each birth and death. Once the probability o xation o an allele is unity, there are no urther changes in allele requencies. We illustrate the solution o (5.24)-(5.25) in Fig. 5.3 or a small population o size = 20, and where the number o -alleles is precisely known in the ounding generation, with either (a) 10 (0) = 1, or; (b) 13 (0) = 1. We plot the probability mass density o the number o -alleles every evolution steps, up to 7 steps, corresponding to approximately ourteen discrete generations o evolution in the Wright-Fisher model. Notice how the probability distribution difuses away rom its initial value, and how the probabilities eventually concentrate on the boundaries, with both 0 and 20 monotonically increasing. In act, ater a large number o generations, 0 approaches the initial requency o the -allele and 20 approaches the initial requency o the -allele (not shown in gures). To better understand this numerical solution, we consider the limit o large (but not innite) populations by expanding (5.24) in powers o 1/ . We rst rewrite (5.24) as ( + 1)
− () = +1(1 − +1) +1() − 2 (1 − ) () + 1 (1 − 1 ) 1 (). −
−
−
(5.26)
84
CHAPTER 5. POPULATION GENETICS
(a) 0.2
0.2
0.1
0.1
0
0 0
5
10
15
20
0.2
0.2
0.1
0.1
0
0
5
10
15
20
0
5
10
15
20
0
5
10
15
20
0
5
10
15
20
0
5
10
15
20
0
5
10
15
20
0
5
10
15
20
0
5
10
15
20
0 0
5
10
15
20
0.2
0.2
0.1
0.1
0
0 0
5
10
15
20
0.2
0.2
0.1
0.1
0
0 0
5
10
15
20
(b) 0.2
0.2
0.1
0.1
0
0 0
5
10
15
20
0.2
0.2
0.1
0.1
0
0 0
5
10
15
20
0.2
0.2
0.1
0.1
0
0 0
5
10
15
20
0.2
0.2
0.1
0.1
0
0 0
5
10
15
20
Figure 5.3: versus with = 20. Evolution steps plotted correspond to = 0, , 2 , . . . , 7 . (a) Initial number o individuals is = 10. (b) Initial number o individuals is = 13.
85
5.5. RANDOM GENETIC DRIFT
We then introduce the continuous random variable = / , with 0 1, and the continuous time = /(/2). The variable corresponds to the requency o the -allele in the population, and a unit o time corresponds to approximately a single discrete generation in the Wright-Fisher model. The probability density unction is dened by
≤ ≤
(, ) = (),
with = /, = 2/.
Furthermore, we dene () = , and note that () = . Similarly, we have +1 = ( + Δ) and −1 = ( Δ), where Δ = 1/ . Then, with Δ = 2/ , (5.26) transorms into
−
(
(, + Δ)
)
− (, ) = ( + Δ) 1 − ( + Δ) ( + Δ, ) − 2 () 1 − () (, ) + ( − Δ) 1 − ( − Δ) ( − Δ, ).
(
)
(
)
(5.27)
To simpliy urther, we use the well-known central-diference approximation to the second-derivative o a unction (), ′′ () =
( + Δ)
− 2 () + ( − Δ) + O(Δ2), Δ2
and recognize the right-hand side o (5.27) to be the numerator o a secondderivative. With Δ = Δ/2 = 1/ 0, and () = , we derive to leadingorder in 1/ the partial diferential equation
→
(, ) 1 2 = () (, ) , 2 2
(
with
)
(5.28)
(1 ) . (5.29) To interpret the meaning o the unction () we will use the ollowing result rom probability theory. For independent trials, each with probability o success and probability o ailure 1 , the number o successes, denoted by , is a binomial random variable with parameters (, ). Well-known results are [ ] = and Var[ ] = (1 ), where [. . . ] is the expected value, and Var[. . . ] is the variance. Now, the number o -alleles chosen when orming the next Wright-Fisher generation is a binomial random variable ′ with parameters (, / ). Thereore, [′ ] = and Var[′ ] = (1 / ) . With ′ = ′ / and = / , we have [′ ] = , and Var[′ ] = (1 )/ . The unction () can thereore be interpreted as the variance o over a single Wright-Fisher generation. Although (5.28) and (5.29) depend explicitly on the population size , the population size can be eliminated by a simple change o variables. I we let
−
() =
−
−
− −
= /,
(5.30)
then the diferential equation (5.28) transorms to (, ) 1 2 = (1 2 2
(
)
− ) (, ) ,
(5.31)
86
CHAPTER 5. POPULATION GENETICS
independent o . The change-o-variables given by (5.30) simply states that a doubling o the population size will lengthen the time scale o evolution by a corresponding actor o two. Remember that here we are already working under the assumption that population sizes are large. Equation (5.31) is a difusion-like equation or the probability distribution unction (, ). A difusion approximation or studying genetic drit was rst introduced into population genetics by its ounders Fisher and Wright, and was later extensively developed in the 1950’s by the Japanese biologist Motoo Kimura. Here, our analysis o this equation relies on a more recent paper by McKane and Waxman (2007), who showed how to construct the analytical solution o (5.31) at the boundaries o . For 0 < < 1, it is suggestive rom the numerical solutions shown in Fig. 5.3 that the probability distribution might asymptotically become independent o . Accordingly, we look or an asymptotic solution to (5.31) at the interior values o satisying (, ) = ( ). Equation (5.31) becomes ( ) 1 ′′ = (1 ) ( ) 2 = ( ),
(
)
−
−
with asymptotic solution ( ) = − ,
(5.32)
and we observe that the total probability over the interior values o decays exponentially. To understand how the solution behaves at the boundaries o , we require boundary conditions or (, ). In act, because (, ) is singular at the boundaries o , appropriate boundary conditions can not be obtained directly rom the diference equations given by (5.25). Rather, boundary conditions are most easily obtained by rst recasting (5.31) into the orm o a continuity equation. We let (, ) denote the so-called probability current. In a small region o size Δ lying in the interval (, +Δ), the time-rate o change o the probability (, )Δ is due to the ow o probability into and out o this region. With an appropriate denition o the probability current, we have in general (, )Δ = (, ) or as Δ
→ 0,
(
)
− ( + Δ),
(, ) (, ) + = 0, (5.33) which is the usual orm o a continuity equation. Identication o (5.33) with (5.31) shows that the probability current o our problem is given by (, ) =
− 12
(
(1
)
− ) (, ) .
Now, since the total probability is unity; that is,
(5.34)
1
0
(, ) = 1,
(5.35)
87
5.5. RANDOM GENETIC DRIFT
probability can not ow in or out o the boundaries o , and we must thereore have (0, ) = 0, (1, ) = 0, (5.36) which are the appropriate boundary conditions or (5.31). We can look or stationary solutions o (5.31). Use o the continuity equation (5.33) together with the boundary conditions (5.36) shows that the stationary solution has zero probability current; that is, () = 0. Integration o () using (5.34) results in (1 ) () = 1 ,
−
where 1 is an integration constant. The readily apparent solution is () = 1 /[(1 )], but there are also two less obvious solutions. Probability distribution unctions are allowed to contain singular solutions corresponding to Dirac delta unctions, and here we consider the possibility o there being Dirac delta unctions at the boundaries. By noticing that both () = 0 a n d (1 ) (1 ) = 0, we can see that the general solution or () can be written as 1 () = + 2 () + 3 (1 ). (1 )
−
−
−
−
−
Requiring () to satisy (5.35) results in 1 = 0 and 2 + 3 = 1. To determine the remaining ree constant, we can compute the mean requency o the -allele in the population. From the continuity equation (5.33), we have
1
1
(, ) =
0
− 0
(, )
1
=
(, )
0
= 0,
where the rst integral on the right-hand-side was done by parts using the boundary conditions given by (5.36), and the second integral was done using (5.34) and the vanishing o (1 ) () on the boundaries. The mean requency o the -allele is thereore a constant—as one would expect or a nondirectional random genetic drit—and we can assume that its initial value is . We thereore obtain or our stationary solution
−
() = (1
− ) () + (1 − ).
The eventual probability o xing the -allele is thereore simply equal to its initial requency. For example, suppose that within a population o individuals homogeneous or , a single neutral mutation occurs so that one individual now carries the -allele. What is the probability that the -allele eventually becomes xed? Our result would yield the probability 1/ , which is the initial requency o the -allele. Intuitively, ater a sucient number o generations has passed, all living individuals should be descendant rom a single ancestral individual living at the time the single mutation occurred. The probability that that single individual carried the -allele is just 1/ . We urther note that Kimura was the rst to nd an analytical solution o the difusion equation (5.31). A solution method using Fourier transorms can be ound in the appendix o McKane and Waxman (2007). In addition to making
88
CHAPTER 5. POPULATION GENETICS
use o Dirac delta unctions, these authors require Heaviside step unctions, Bessel unctions, spherical Bessel unctions, hypergeometric unctions, Legendre polynomials, and Gegenbauer polynomials. The resulting solution is o the orm (, ) = Π0 ( ) () + Π 1 ( ) (1
− ) + (, ).
(5.37)
I the number o -alleles are known at the initial instant, and the requency is , then (, 0) = ( ),
−
and Π0 (0) = Π1 (0) = 0; (, 0) = ( ). As , we have (, ) 0; Π0 ( ) 1 and Π1 ( ) . We can easily demonstrate at least one exact solution o the orm ( 5.37). For an initial uniorm probability distribution given by
→ −
→
−
→∞
→
(, 0) = 1, the probability distribution at the interior points remains uniorm and is given by (5.32) with = 1. I we assume the orm o solution given by ( 5.37) together with the requirement o unit probability given by (5.35), we can obtain the exact result 1 1 − () + (1 ) + − . (, ) = 2
( − )(
−
)
Reerences Kimura, M. Solution o a process o random genetic drit with a continuous model. Proceeding o the National Academy o Sciences, USA (1955) 41, 144150. McKane, A.J. & Waxman, D. Singular solutions o the difusion equation o population genetics. Journal o Theoretical Biology (2007) 247, 849-858.
Chapter 6
Biochemical Reactions Biochemistry is the study o the chemistry o lie. It can be considered a branch o molecular biology, perhaps more ocused on specic molecules and their reactions, or a branch o chemistry ocused on the complex chemical reactions occurring in living organisms. Perhaps the rst application o biochemistry happened about 5000 years ago when bread was made using yeast. Modern biochemistry, however, had a relatively slow start among the sciences, as did modern biology. Isaac Newton’s publication o Principia Mathematica in 1687 preceded Darwin’s Origin o Species in 1859 by almost 200 years. I nd this amazing because the ideas o Darwin are in many ways simpler and easier to understand than the mathematical theory o Newton. Most o the delay must be attributed to a undamental conict between science and religion. The physical sciences experienced this conict early—witness the amous prosecution o Galileo by the Catholic Church in 1633, during which Galileo was orced to recant his heliocentric view—but the conict o religion with evolutionary biology continues even to this day. Advances in biochemistry were initially delayed because it was long believed that lie was not subject to the laws o science the way nonlie was, and that only living things could produce the molecules o lie. Certainly, this was more a religious conviction than a scientic one. Then Friedrich W¨ohler in 1828 published his landmark paper on the synthesis o urea (a waste product neutralizing toxic ammonia beore excretion in the urine), demonstrating or the rst time that organic compounds can be created articially. Here, we present mathematical models or some important biochemical reactions. We begin by introducing a useul model or a chemical reaction: the law o mass action . We then model what may be the most important biochemical reactions, namely those catalyzed by enzymes. Using the mathematical model o enzyme kinetics, we consider three undamental enzymatic properties: competitive inhibition , allosteric inhibition , and cooperativity .
6.1
The law o mass action
The law o mass action describes the rate at which chemicals interact in reactions. It is assumed that diferent chemical molecules come into contact by collision beore reacting, and that the collision rate is directly proportional to 89
90
CHAPTER 6. BIOCHEMICAL REACTIONS
the number o molecules o each reacting species. Suppose that two chemicals and react to orm a product chemical , written as +
→ ,
with the rate constant o the reaction. For simplicity, we will use the same symbol , say, to reer to both the chemical and its concentration. The law o mass action says that / is proportional to the product o the concentrations and , with proportionality constant . That is, = .
(6.1)
Similarly, the law o mass action enables us to write equations or the timederivatives o the reactant concentrations and : =
−,
=
−.
(6.2)
Notice that when using the law o mass action to nd the rate-o-change o a concentration, the chemical that the arrow points towards is increasing in concentration (positive sign), the chemical that the arrow points away rom is decreasing in concentration (negative sign). The product o concentrations on the right-hand-side is always that o the reactants rom which the arrow points away, multiplied by the rate constant that is on top o the arrow. Equation (6.1) can be solved analytically using conservation laws. Each reactant, original and converted to product, is conserved since one molecule o each reactant gets converted into one molecule o product. Thereore, ( + ) = 0 ( + ) = 0
=
⇒ =⇒
+ = 0 , + = 0 ,
where 0 and 0 are the initial concentrations o the reactants, and no product is present initially. Using the conservation laws, (6.1) becomes = (0
− )(0 − ), with (0) = 0,
which may be integrated by separating variables. Ater some algebra, the solution is determined to be (0 −0 ) () = 0 0 0 (0 −0 )
−1 , − 0
which is a complicated expression with the simple limits lim () =
→∞
0 0
i 0 < 0 , i 0 < 0 .
(6.3)
The reaction stops ater one o the reactants is depleted; and the nal concentration o the product is equal to the initial concentration o the depleted reactant.
91
6.1. THE LAW OF MASS ACTION
I we also include the reverse reaction, + +
-
,
− then the time-derivative o the product is given by = +
−
−
.
˙ = 0, and using Notice that + and − have diferent units. At equilibrium, the conservation laws + = 0 , + = 0 , we obtain (0
− )(0 − ) − + = 0, −
rom which we dene the equilibrium constant by = − /+ , which has units o concentration. Thereore, at equilibrium, the concentration o the product is given by the solution o the quadratic equation 2
− (0 + 0 + ) + 00 = 0,
with the extra condition that 0 < < min(0 , 0 ). For instance, i 0 = 0 0 , then at equilibrium, = 0
−
1 2
1 + 40 /
≡
−
1 .
I 0 , then and have a high anity, and the reaction proceeds mainly to , with 0 . Below are two interesting reactions. In reaction (ii), is assumed to be held at a constant concentration. (i) +
≪
→
+
-
2
− (ii) +
→ 2, 1
+
→ 2, 2
→ 3
˙ in reaction (i), and ˙ and ˙ in reaction Can you write down the equations or (ii)? When normalized properly, the equations rom reaction (ii) reduce to the Lotka-Volterra predator-prey equations introduced in 1.4. The chemical concentrations and , thereore, oscillate in time like predators and their prey.
S
92
CHAPTER 6. BIOCHEMICAL REACTIONS
Figure 6.1: A Michaelis-Menten reaction o two substrates converting to one product. (Drawn by User:IMeowbot, released under the GNU Free Documentation License.)
6.2
Enzyme kinetics
Enzymes are catalysts, usually proteins, that help convert other molecules called substrates into products, but are themselves unchanged by the reaction. Each enzyme has high specicity or at least one reaction, and it can accelerate this reaction by millions o times. Without enzymes, most biochemical reactions are too slow or lie to be possible. Enzymes are so important to our lives that a single amino acid mutation in one enzyme out o the more than 2000 enzymes in our bodies can result in a severe or lethal genetic disease. Enzymes do not ollow the law o mass action directly: with substrate, product, and enzyme, the reaction +
→ + ,
is a poor model since the reaction velocity / is known to attain a nite limit with increasing substrate concentration. Rather, Michaelis and Menten (1913) proposed the ollowing reaction scheme with an intermediate molecule: 1 +
-
2
-
+ ,
−1 where is a complex ormed by the enzyme and the substrate. A cartoon o the Michaelis-Menten reaction with an enzyme catalyzing a reaction between two substrates is shown in Fig. 6.1. Commonly, substrate is continuously provided to the reaction and product is continuously removed. The removal o product has been modeled by neglecting the reverse reaction + . A continuous provision o substrate allows us to assume that is held at an approximately constant concentration. The diferential equations or and can be obtained rom the law o mass action:
→
/ = 1 / = 2 .
− (
−1
+ 2 ),
Biochemists usually want to determine the reaction velocity / in terms o the substrate concentration and the total enzyme concentration 0 . We can
93
6.3. COMPETITIVE INHIBITION
eliminate in avor o 0 rom the conservation law that the enzyme, ree and bound, is conserved; that is ( + ) =0
=
⇒
+ = 0
=
⇒
= 0
− ;
and we can rewrite the equation or / eliminating : = 1 ( 0 ) (−1 + 2 ) = 1 0 (−1 + 2 + 1 ).
− − −
(6.4)
Because is assumed to be held constant, the complex is expected to be in equilibrium, with the rate o ormation equal to the rate o dissociation. With ˙ = 0 in this so-called quasi-steady-state approximation, we may assume that (6.4), and we have 1 0 = . −1 + 2 + 1 The reaction velocity is then given by = 2 1 2 0 = −1 + 2 + 1 = , +
(6.5)
where two undamental constants are dened: = (−1 + 2 )/1 ,
= 2 0 .
(6.6)
The Michaelis-Menten constant or the Michaelis constant has units o concentration, and the maximum reaction velocity has units o concentration divided by time. The interpretation o these constants is obtained by considering the ollowing limits: as
→ ∞,
i = ,
→ 0
=
1 0 2
and /
→ ,
and / =
1 . 2
Thereore, is the limiting reaction velocity obtained by saturating the reaction with substrate so that every enzyme is bound; and is the concentration o at which only one-hal o the enzymes are bound and the reaction proceeds at one-hal maximum velocity.
6.3
Competitive inhibition
Competitive inhibition occurs when inhibitor molecules compete with substrate molecules or binding to the same enzyme’s active site. When an inhibitor is bound to the enzyme, no product is produced so competitive inhibition will
94
CHAPTER 6. BIOCHEMICAL REACTIONS
Figure 6.2: Competitive inhibition. (Drawn by G. Andruk, released under the GNU Free Documentation License.) reduce the velocity o the reaction. A cartoon o this process is shown in Fig. 6.2. To model competitive inhibition, we introduce an additional reaction associated with the inhibitor-enzyme binding: 1 +
-
1
2
-
+ ,
−1 3 +
-
2 .
−3 With more complicated enzymatic reactions, the reaction schematic becomes dicult to interpret. Perhaps an easier way to visualize the reaction is rom the ollowing redrawn schematic: 1
-
1
2
- +
6 −1 −3
3
? 2 Here, the substrate and inhibitor are combined with the relevant rate constants, rather than treated separately. It is immediately obvious rom this redrawn schematic that inhibition is accomplished by sequestering enzyme in the orm o 2 and preventing its participation in the catalysis o to .
95
6.3. COMPETITIVE INHIBITION
˙ in terms o the substrate Our goal is to determine the reaction velocity and inhibitor concentrations, and the total concentration o the enzyme (ree and bound). The law o mass action applied to the complex concentrations and the product results in 1 = 1 (−1 + 2 ) 1 , 2 = 3 −3 2 , = 2 1 .
− −
The enzyme, ree and bound, is conserved so that ( + 1 + 2 ) = 0
=
+ 1 + 2 = 0
⇒
=
⇒
= 0
− 1 − 2.
˙ 1 = ˙ 2 = 0, so that Under the quasi-equilibrium approximation, 1 ( 0 1 2 ) 3 ( 0 1
− − − ( 1 + 2) 1 = 0, − − 2) − 3 2 = 0, −
−
which results in the ollowing system o two linear equations and two unknowns ( 1 and 2 ): (−1 + 2 + 1 ) 1 + 1 2 = 1 0 , 3 1 + (−3 + 3 ) 2 = 3 0 .
(6.7) (6.8)
We dene the Michaelis-Menten constant as beore, and an additional constant associated with the inhibitor reaction: =
−1 + 2 , 1
=
−3 . 3
Dividing (6.7) by 1 and (6.8) by 3 yields ( + ) 1 + 2 = 0 ,
(6.9)
1 + ( + ) 2 = 0 .
(6.10)
Since our goal is to obtain the velocity o the reaction, which requires determining 1 , we multiply (6.9) by ( + ) and (6.10) by , and subtract:
−
( + )( + ) 1 + ( + ) 2 = 0 ( + ) 1 + ( + ) 2 = 0
(
( + )( + )
)
− 1 = 0 ;
or ater cancelation and rearrangement 0 + + 0 = . (1 + / ) +
1 =
96
CHAPTER 6. BIOCHEMICAL REACTIONS
Figure 6.3: Allosteric inhibition. (Unknown artist, released under the GNU Free Documentation License.)
Thereore, the reaction velocity is given by (2 0 ) = (1 + / ) + = ′ , +
(6.11)
where = 2 0 ,
′ = (1 + / ).
(6.12)
By comparing the inhibited reaction velocity (6.11) and (6.12) with the uninhibited reaction velocity (6.5) and (6.6), we observe that inhibition increases the Michaelis-Menten constant o the reaction, but leaves unchanged the maximum reaction velocity. Since the Michaelis-Menten constant is dened as the substrate concentration required to attain one-hal o the maximum reaction velocity, addition o an inhibitor with a xed substrate concentration acts to decrease the reaction velocity. However, a reaction saturated with substrate still attains the uninhibited maximum reaction velocity.
6.4
Allosteric inhibition
The term allostery comes rom the Greek word allos , meaning diferent, and stereos , meaning solid, and reers to an enzyme with a regulatory binding site separate rom its active binding site. In our model o allosteric inhibition , an inhibitor molecule is assumed to bind to its own regulatory site on the enzyme, resulting in either a lowered binding anity o the substrate to the enzyme, or a lowered conversion rate o substrate to product. A cartoon o allosteric inhibition due to a lowered binding anity is shown in Fig. 6.3. In general, we need to dene three complexes: 1 is the complex ormed rom substrate and enzyme; 2 rom inhibitor and enzyme, and; 3 rom substrate, inhibitor, and enzyme. We write the chemical reactions as ollows:
97
6.4. ALLOSTERIC INHIBITION
1
-
6 −1
1′
? 2
- +
6 ′ − 3
−3 3′
3
2
1
- ?
2′
3
- + 2
′ − 1
The general model or allosteric inhibition with ten independent rate constants appears too complicated to analyze. We will simpliy this general model to one with ewer rate constants that still exhibits the unique eatures o allosteric inhibition. One possible but uninteresting simplication assumes that i binds to , then does not; however, this reduces allosteric inhibition to competitive inhibition and loses the essence o allostery. Instead, we simpliy by allowing both and to simultaneously bind to , but we assume that the binding o prevents substrate conversion to product. With this simplication, 2′ = 0. To urther reduce the number o independent rate constants, we assume that the binding o to is unafected by the bound presence o , and the binding o to is unafected by the bound presence o . These approximations imply that all the primed rate constants equal the corresponding unprimed rate constants, e.g., 1′ = 1 , etc. With these simplications, the schematic o the chemical reaction simplies to
1
-
6 −1
?
1
2
- +
6
−3 3
3
2
1
−3
- ?
3
−1 and now there are only ve independent rate constants. We write the equations or the complexes using the law o mass action: 1 = 1 + −3 3 (−1 + 2 + 3 ) 1 , 2 = 3 + −1 3 (−3 + 1 ) 2 , 3 = 3 1 + 1 2 (−1 + −3 ) 3 ,
− − −
(6.13) (6.14) (6.15)
98
CHAPTER 6. BIOCHEMICAL REACTIONS
while the reaction velocity is given by = 2 1 .
(6.16)
Again, both ree and bound enzyme is conserved, so that = 0 1 2 3 . ˙ 1 = ˙ 2 = ˙ 3 = 0, we obtain a With the quasi-equilibrium approximation system o three equations and three unknowns: 1 , 2 and 3 . Despite our simplications, the analytical solution or the reaction velocity remains messy (see Keener & Sneyd, reerenced at the chapter’s end) and not especially illuminating. We omit the complete analytical result here and determine only the maximum reaction velocity. The maximum reaction velocity ′ or the allosteric-inhibited reaction is dened as the time-derivative o the product concentration when the reaction is saturated with substrate; that is,
− − −
′ = lim / →∞
= 2 lim 1 . →∞
With substrate saturation, every enzyme will have its substrate binding site occupied. Enzymes are either bound with only substrate in the complex 1 , or bound together with substrate and inhibitor in the complex 3 . Accordingly, the schematic o the chemical reaction with substrate saturation simplies to
2
1
- + 1
6 −3
3
? 3 The equations or 1 and 3 with substrate saturation are thus given by 1 = −3 3 3 = 3 1
− 3 1, − 3 3, −
(6.17) (6.18)
and the quasi-equilibrium approximation yields the single independent equation 3 = (3 /−3 ) 1 = (/ ) 1 ,
(6.19)
with = −3 /3 as beore. The equation expressing the conservation o enzyme is given by 0 = 1 + 3 . This conservation law, together with (6.19), permits us to solve or 1 : 0 1 = . 1 + /
99
6.5. COOPERATIVITY
Figure 6.4: Cooperativity. Thereore, the maximum reaction velocity or the allosteric-inhibited reaction is given by 2 0 1 + / = , 1 + /
′ =
where is the maximum reaction velocity o both the uninhibited and the competitive inhibited reaction. The allosteric inhibitor is thus seen to reduce the maximum velocity o the uninhibited reaction by the actor (1 + / ), which may be large i the concentration o allosteric inhibitor is substantial.
6.5
Cooperativity
Enzymes and other protein complexes may have multiple binding sites, and when a substrate binds to one o these sites, the other sites may become more active. A well-studied example is the binding o the oxygen molecule to the hemoglobin protein. Hemoglobin can bind our molecules o O 2 , and when three molecules are bound, the ourth molecule has an increased anity or binding. We call this cooperativity . We will model cooperativity by assuming that an enzyme has two separated but indistinguishable binding sites or a substrate . For example, the enzyme may be a protein dimer, composed o two identical subproteins with identical binding sites or . A cartoon o this enzyme is shown in Fig. 6.4. Because the two binding sites are indistinguishable, we need consider only two complexes: 1 and 2 , with enzyme bound to one or two substrate molecules, respectively. When the enzyme exhibits cooperativity, the binding o the second substrate molecule has a greater rate constant than the binding o the rst. We thereore consider the ollowing reaction:
100
CHAPTER 6. BIOCHEMICAL REACTIONS
1
-
−1
2
1
- +
6 −3
3
? 2
4
- + 1
where cooperativity supposes that 1 action results in
≪ 3.
Application o the law o mass
1 = 1 + (−3 + 4 ) 2 (−1 + 2 + 3 ) 1 , 2 = 3 1 (−3 + 4 ) 2 .
−
−
˙ 1 = ˙ 2 = 0 and the conservaApplying the quasi-equilibrium approximation tion law 0 = + 1 + 2 results in the ollowing system o two equations and two unknowns:
(
)
−1 + 2 + (1 + 3 ) 1
− (
+ 4 1 ) 2 = 1 0 , 3 1 (−3 + 4 ) 2 = 0. −3
−
−
(6.20) (6.21)
We divide (6.20) by 1 and (6.21) by 3 and dene 1 =
−1 + 2 , 1
2 =
−3 + 4 , 3
= 1 /3
to obtain
(
)
1 + (1 + ) 1
− ( 2 − ) 2 = 0, 1 − 2 2 = 0.
(6.22) (6.23)
We can subtract (6.23) rom (6.22) and cancel to obtain ( 1 + ) 1 + 2 = 0 .
(6.24)
Equations (6.23) and (6.24) can be solved or 1 and 2 : 2 0 , 1 2 + 2 + 2 0 2 2 = , 1 2 + 2 + 2 1 =
(6.25) (6.26)
so that the reaction velocity is given by = 2 1 + 4 2 (2 2 + 4 ) 0 = . 1 2 + 2 + 2
(6.27)
101
6.5. COOPERATIVITY
To illuminate this result, we consider two limiting cases: (i) no cooperativity, where the active sites act independently so that each protein dimer, say, can be considered as two independent protein monomers; (ii) strong cooperativity, where the binding o the second substrate has a much greater rate constant than the binding o the rst. Independent active sites The ree enzyme has two independent binding sites while 1 has only a single binding site. Consulting the reaction schematic: 1 is the rate constant or the binding o to two independent binding sites; −1 and 2 are the rate constants or the dissociation and conversion o a single rom the enzyme; 3 is the rate constant or the binding o to a single ree binding site, and; −3 and 4 are the rate constants or the dissociation and conversion o one o two independent ’s rom the enzyme. Accounting or these actors o two and assuming independence o active sites, we have 1 = 23 ,
−3 = 2−1 ,
4 = 22 .
We dene the Michaelis-Menten constant that is representative o the protein monomer with one binding site; that is, −1 + 2 1 /2 = 2 1 1 = 2 . 2
=
Thereore, or independent active sites, the reaction velocity becomes (22 + 22 ) 0 = 2 + 2 + 2 22 0 = . + The reaction velocity or a dimer protein enzyme composed o independent identical monomers is simply double that o a monomer protein enzyme, an intuitively obvious result. Strong cooperativity We now assume that ater the rst substrate binds to the enzyme, the second substrate binds much more easily, so that 1 3 . The number o enzymes bound to a single substrate molecule should consequently be much less than the number bound to two substrate molecules, resulting in 1 2 . Dividing (6.25) by (6.26), this inequality becomes
≪
≪
1 2 = 2
≪ 1.
Dividing the numerator and denominator o (6.27) by 2 , we have (2 2 / + 4 ) 0 = . ( 1 / )( 2 / ) + ( 2 / ) + 1
102
CHAPTER 6. BIOCHEMICAL REACTIONS
1
n=2 0.9 0.8 0.7
n=1
0.6
t d / 0.5 P d 0.4 0.3 0.2 0.1 0 0
2
4
6
8
10
S
Figure 6.5: The reaction velocity / as a unction o the substrate . Shown are the solutions to the Hill equation with = 1, = 1, and or = 1, 2. To take the limit o this expression as 2 / 0, we set 2 / = 0 everywhere except in the last term in the denominator, since 1 / is inversely proportional to 1 and may go to innity in this limit. Taking the limit and multiplying the numerator and denominator by 2 ,
→
4 0 2 = . 1 2 + 2 Here, the maximum reaction velocity is = 4 0 , and the modied MichaelisMenten constant is = 1 2 , so that
√
2 = 2 . + 2 In biochemistry, this reaction velocity is generalized to = , + known as the Hill equation , and by varying is used to t experimental data. In Fig. 6.5, we have plotted the reaction velocity / versus as obtained rom the Hill equation with = 1 or 2. In drawing the gure, we have taken both and equal to unity. It is evident that with increasing the reaction velocity more rapidly saturates to its maximum value.
Reerences Keener, J. & Sneyd, J. Mathematical Physiology. Springer-Verlag, New York (1998). Pg. 30, Exercise 2.
Chapter 7
Sequence Alignment The sotware program BLAST (Basic Local Alignment Search Tool) uses sequence alignment algorithms to compare a query sequence against a database to identiy other known sequences similar to the query sequence. Oten, the annotations attached to the already known sequences yield important biological inormation about the query sequence. Almost all biologists use BLAST, making sequence alignment one o the most important algorithms o bioinormatics. The sequence under study can be composed o nucleotides (rom the nucleic acids DNA or RNA) or amino acids (rom proteins). Nucleic acids chain together our diferent nucleotides: A,C,T,G or DNA and A,C,U,G or RNA; proteins chain together twenty diferent amino acids. The sequence o a DNA molecule or o a protein is the linear order o nucleotides or amino acids in a specied direction, dened by the chemistry o the molecule. There is no need or us to know the exact details o the chemistry; it is sucient to know that a protein has distinguishable ends called the N-terminus and the C-terminus, and that the usual convention is to read the amino acid sequence rom the N-terminus to the C-terminus. Specication o the direction is more complicated or a DNA molecule than or a protein molecule because o the double helix structure o DNA, and this will be explained in Section 7.1. The basic sequence alignment algorithm aligns two or more sequences to highlight their similarity, inserting a small number o gaps into each sequence (usually denoted by dashes) to align wherever possible identical or similar characters. For instance, Fig 7.1 presents an alignment using the sotware tool ClustalW o the hemoglobin beta-chain rom a human, a chimpanzee, a rat, and a zebrash. The human and chimpanzee sequences are identical, a consequence o our very close evolutionary relationship. The rat sequence difers rom human/chimpanzee at only 27 out o 146 amino acids; we are all mammals. The zebrash sequence, though clearly related, diverges signicantly. Notice the insertion o a gap in each o the mammal sequences at the zebra sh amino acid position 122. This permits the subsequent zebrash sequence to better align with the mammal sequences, and implies either an insertion o a new amino acid in sh, or a deletion o an amino acid in mammals. The insertion or deletion o a character in a sequence is called an indel . Mismatches in sequence, such as that occurring between zebrash and mammals at amino acid positions 2 and 3 is called a mutation . ClustalW places a ‘*’ on the last line to denote exact amino acid matches across all sequences, and a ‘:’ and ‘.’ to 103
104
CHAPTER 7. SEQUENCE ALIGNMENT
CLUSTAL W (1.83) multiple sequence alignment
Human Chimpanzee Rat Zebrafish
VHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESFGDLSTPDAVMGNPKV VHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESFGDLSTPDAVMGNPKV VHLTDAEKAAVNGLWGKVNPDDVGGEALGRLLVVYPWTQRYFDSFGDLSSASAIMGNPKV VEWTDAERTAILGLWGKLNIDEIGPQALSRCLIVYPWTQRYFATFGNLSSPAAIMGNPKV *. * *::*: .****:* *::* :**.* *:*******:* :**:**:. *:******
60 60 60 60
Human Chimpanzee Rat Zebrafish
KAHGKKVLGAFSDGLAHLDNLKGTFATLSELHCDKLHVDPENFRLLGNVLVCVLAHHFGK KAHGKKVLGAFSDGLAHLDNLKGTFATLSELHCDKLHVDPENFRLLGNVLVCVLAHHFGK KAHGKKVINAFNDGLKHLDNLKGTFAHLSELHCDKLHVDPENFRLLGNMIVIVLGHHLGK AAHGRTVMGGLERAIKNMDNVKNTYAALSVMHSEKLHVDPDNFRLLADCITVCAAMKFGQ ***:.*:..:. .: ::**:*.*:* ** :*.:******:*****.: :. . ::*:
120 120 120 120
Human Chimpanzee Rat Zebrafish
E-FTPPVQAAYQKVVAGVANALAHKYH E-FTPPVQAAYQKVVAGVANALAHKYH E -FTPCAQAAFQKV VAGVASALAHKY H AGFNADVQEAWQKFLAVVVSALCRQYH *.. .* *:**.:* *..**.::**
146 146 146 147
Figure 7.1: Multiple alignment o the hemoglobin beta-chain or Human, Chimpanzee, Rat and Zebra sh, obtained using ClustalW. denote chemically similar amino acids across all sequences (each amino acid has characteristic chemical properties, and amino acids can be grouped according to similar properties). In this chapter, we detail the algorithms used to align sequences.
7.1
The minimum you need to know about DNA chemistry and the genetic code
In one o the most important scientic papers ever published, James Watson and Francis Crick, pictured in Fig. 7.2, determined the structure o DNA using a three-dimensional molecular model that makes plain the chemical basis o heredity. The DNA molecule consists o two strands wound around each other to orm the now amous double helix. Arbitrarily, one strand is labeled by the sequencing group to be the positive strand, and the other the negative strand. The two strands o the DNA molecule bind to each other by base pairing: the bases o one strand pair with the bases o the other strand. Adenine (A) always pairs with thymine (T), and guanine (G) always pairs with cytosine (C): A with T, G with C. For RNA, T is replaced by uracil (U). When reading the sequence o nucleotides rom a single strand, the direction o reading must be specied, and this is possible by reerring to the chemical bonds o the DNA backbone. There are o course only two possible directions to read a linear sequence o bases, and these are denoted as 5’-to-3’ and 3’-to-5’. Importantly, the two separate strands o the DNA molecule are oriented in opposite directions. Below is the beginning o the DNA coding sequence or the human hemoglobin beta chain protein discussed earlier: 5’-GTGCACCTGACTCCTGAGGAGAAGTCTGCCGTTACTGCCCTGTGGGGCAAGGTGAACGTG-3’ |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| 3’-CACGTGGACTGAGGACTCCTCTTCAGACGGCAATGACGGGACACCCCGTTCCACTTGCAC-5’
105
7.2. BRUTE FORCE ALIGNMENT
Figure 7.2: James Watson and Francis Crick posing in ront o their DNA model. The original photograph was taken in 1953, the year o discovery, and was recreated in 2003, ty years later. Francis Crick, the man on the right, died in 2004. It is important to realize that there are two unique DNA sequences here, and either one, or even both, can be coding. Reading rom 5’-to-3’, the upper sequence begins ‘GTGCACCTG...’, while the lower sequence ends ‘...CAGGTGCAC’. Here, only the upper sequence codes or the human hemoglobin beta chain, and the lower sequence is non-coding. How is the DNA code read? Enzymes separate the two strands o DNA, and transcription occurs as the DNA sequence is copied into messenger RNA (mRNA). I the upper strand is the coding sequence, then the complementary lower strand serves as the template or constructing the mRNA sequence. The ACUG nucleotides o the mRNA bind to the lower sequence and construct a single stranded mRNA molecule containing the sequence ‘GUGCACCUG...’, which exactly matches the sequence o the upper, coding strand, but with T replaced by U. This mRNA is subsequently translated in the ribosome o the cell, where each nucleotide triplet codes or a single amino acid. The triplet coding o nucleotides or amino acids is the amous genetic code , shown in Fig. 7.3. Here, the translation to amino acid sequence is ‘VHL...’, where we have used the genetic code ‘GUG’ = V, ‘CAC’=H, ‘CUG’=L. The three out o the twenty amino acids used here are V = Valine, H = Histidine, and L = Leucine.
7.2
Sequence alignment by brute orce
One (bad) approach to sequence alignment is to align the two sequences in all possible ways, score the alignments with an assumed scoring system, and determine the highest scoring alignment. The problem with this brute-orce approach is that the number o possible alignments grows exponentially with sequence length; and or sequences o reasonable length, the computation is already impossible. For example, the number o ways to align two sequences o 50 characters each—a rather small alignment problem—is about 1 .5 1037 , already an astonishingly large number. It is inormative to count the number
×
106
CHAPTER 7. SEQUENCE ALIGNMENT
Figure 7.3: The genetic code. o possible alignments between two sequences since a similar algorithm is used or sequence alignment. Suppose we want to align two sequences. Gaps in either sequence are allowed but a gap can not be aligned with a gap. By way o illustration, we demonstrate the three ways that the rst character o the upper-case alphabet and the lowercase alphabet may align: A -A A| , || , || a a-a
and the ve ways in which the rst two characters o the upper-case alphabet can align with the rst character o the lower-case alphabet: AB AB ABA-B -AB || , || , ||| , ||| , ||| . -a a--a -aa--
A recursion relation or the total number o possible alignments o a sequence o characters with a sequence o characters may be derived by considering the alignment o the last character. There are three possibilities that we illustrate by assuming the th character is ‘F’ and the th character is ‘d’: (1) 1 characters o the rst sequence are already aligned with 1 characters o the second sequence, and the th character o the rst sequence aligns exactly with the th character o the second sequence:
−
...F |||| ...d
−
107
7.2. BRUTE FORCE ALIGNMENT
(2) 1 characters o the rst sequence are aligned with characters o the second sequence and the th character o the rst sequence aligns with a gap in the second sequence:
−
...F |||| ...-
(3) characters o the rst sequence are aligned with 1 characters o the second sequence and a gap in the rst sequence aligns with the th character o the second sequence:
−
...|||| ...d
I (, ) is the number o ways to align an character sequence with a character sequence, then, rom our counting, (, ) = (
− 1, − 1) + ( − 1, ) + (, − 1).
(7.1)
This recursion relation requires boundary conditions. Because there is only one way to align an > 0 character sequence against a zero character sequence (i.e., characters against gaps) the boundary conditions are (0, ) = (, 0) = 1 or all , > 0 We may also add the additional boundary condition (0, 0) = 1, obtained rom the known result (1, 1) = 3. Using the recursion relation (7.1), we can construct the ollowing dynamic matrix to count the number o ways to align the two ve-character sequences 1 2 3 4 5 and 1 2 3 4 5 : 1 2 3 4 5
1 1 1 1 1 1
1 1 3 5 7 9 11
2 1 5 13 25 41 61
3 1 7 25 63 129 231
4 1 9 41 129 321 681
5 1 11 61 231 681 1683
The size o this dynamic matrix is 6 6, and or convenience we label the rows and columns starting rom zero (i.e., row 0, row 1, . . . , row 5). This matrix was constructed by rst writing 1 2 3 4 5 to the let o the matrix and 1 2 3 4 5 above the matrix, then lling in ones across the zeroth row and down the zeroth column to satisy the boundary conditions, and nally applying the recursion relation directly by going across the rst row rom let-to-right, the second row rom let-to-right, etc. To demonstrate the lling in o the matrix, we have across the rst row: 1 + 1 + 1 = 3, 1 + 1 + 3 = 5, 1 + 1 + 5 = 7, etc, and across the second row: 1 + 3 + 1 = 5, 3 + 5 + 5 = 13, 5 + 7 + 13 = 25, etc. Finally, the last element entered gives the number o ways to align two ve character sequences: 1683, already a remarkably large number. It is possible to solve analytically the recursion relation (7.1) or (, ) using generating unctions. Although the solution method is interesting—and in act was shown to me by a student—the nal analytical result is messy and we omit it here. In general, computation o (, ) is best done numerically by constructing the dynamic matrix.
×
−
−
108
CHAPTER 7. SEQUENCE ALIGNMENT
7.3
Sequence alignment by dynamic programming
Two reasonably sized sequences cannot be aligned by brute orce. Luckily, there is another algorithm borrowed rom computer science, dynamic programming , that makes use o a dynamic matrix. What is needed is a scoring system to judge the quality o an alignment. The goal is to nd the alignment that has the maximum score. We assume that the alignment o character with character has the score ( , ). For example, when aligning two DNA sequences, a match (A-A, C-C, T-T, G-G) may be scored as +2, and a mismatch (A-C, A-T, A-G, etc.) scored as 1. We also assume that an indel (a nucleotide aligned with a gap) is scored as , with a typical value or DNA alignment being = 2. In the next section, we develop a better and more widely used model or indel scoring that distinguishes gap openings rom gap extensions. Now, let (, ) denote the maximum score or aligning a sequence o length with a sequence o length . We can compute (, ) provided we know ( 1, 1), ( 1, ) and (, 1). Indeed, our logic is similar to that used when counting the total number o alignments. There are again three ways to compute (, ): (1) 1 characters o the rst sequence are aligned with 1 characters o the second sequence with maximum score ( 1, 1), and the th character o the rst sequence aligns with the th character o the second sequence with updated maximum score ( 1, 1)+ ( , ); (2) 1 characters o the rst sequence are aligned with characters o the second sequence with maximum score ( 1, ), and the th character o the rst sequence aligns with a gap in the second sequence with updated maximum score ( 1, ) + , or; (3) characters o the rst sequence are aligned with 1 characters o the second sequence with maximum score (, 1), and a gap in the rst sequence aligns with the th character o the second sequence with updated maximum score (, 1) + . We then compare these three scores and assign (, ) to be the maximum; that is,
−
−
−
− −
−
−
−
− − −
− −
−
−
−
−
−
(, ) = max
⎧⎪⎨ ⎪⎩
( 1, 1) + ( , ), ( 1, ) + , (, 1) + .
− −
−
(7.2)
−
Boundary conditions give the score o aligning a sequence with a null sequence o gaps, so that (, 0) = (0, ) = , > 0, (7.3) with (0, 0) = 0. The recursion (7.2), together with the boundary conditions (7.3), can be used to construct a dynamic matrix. The score o the best alignment is then given by the last lled-in element o the matrix, which or aligning a sequence o length with a sequence o length is (, ). Besides this score, however, we also want to determine the alignment itsel. The alignment can be obtained by tracing back the path in the dynamic matrix that was ollowed to compute each matrix element (, ). There could be more than one path, so that the best alignment may be degenerate.
109
7.3. DYNAMIC PROGRAMMING
Sequence alignment is always done computationally, and there are excellent sotware tools reely available on the web (see 7.6). Just to illustrate the dynamic programming algorithm, we compute by hand the dynamic matrix or aligning two short DNA sequences GGAT and GAATT, scoring a match as +2, a mismatch as 1 and an indel as 2:
S
−
G G A T
0 -2 -4 -6 -8
G -2 2 0 -2 -4
−
A -4 0 1 2 0
A -6 -2 -1 3 1
T -8 -4 -3 1 5
T -10 -6 -5 -1 3
In our hand calculation, the two sequences to be aligned go to the let and above the dynamic matrix, leading with a gap character ‘-’. Row 0 and column 0 are then lled in with the boundary conditions, starting with 0 in position (0, 0) and incrementing by the gap penalty 2 across row 0 and down column 0. The recursion relation (7.2) is then used to ll in the dynamic matrix one row at a time moving rom let-to-right and top-to-bottom. To determine the (, ) matrix element, three numbers must be compared and the maximum taken: (1) inspect the nucleotides to the let o row and above column and add +2 or a match or -1 or a mismatch to the ( 1, 1) matrix element; (2) add 2 to the ( 1, ) matrix element; (3) add 2 to the (, 1) matrix element. For example, the rst computed matrix element 2 at position (1, 1) was determined by taking the maximum o (1) 0+2 = 2, since G-G is a match; (2) 2 2 = 4; (3) 2 2 = 4. You can test your understanding o dynamic programming by computing the other matrix elements. Ater the matrix is constructed, the traceback algorithm that nds the best alignment starts at the bottom-right element o the matrix, here the (4, 5) matrix element with entry 3. The matrix element used to compute 3 was either at (4, 4) (horizontal move) or at (3, 4) (diagonal move). Having two possibilities implies that the best alignment is degenerate. For now, we arbitrarily choose the diagonal move. We build the alignment rom end to beginning with GGAT on top and GAATT on bottom:
−
− − −
−
−−
−
−
−−
−
−
T | T
We illustrate our current position in the dynamic matrix by eliminating all the elements that are not on the traceback path and are no longer accessible:
0 -2 -4 -6
-2 2 0 -2
-4 0 1 2
-6 -2 -1 3
-8 -4 -3 1
3
We start again rom the 1 entry at (3, 4). This value came rom the 3 entry at (3, 3) by a horizontal move. Thereore, the alignment is extended to
110
CHAPTER 7. SEQUENCE ALIGNMENT
-T || TT
where a gap is inserted in the top sequence or a horizontal move. (A gap is inserted in the bottom sequence or a vertical move.) The dynamic matrix now looks like
0 -2 -4 -6
-2 2 0 -2
-4 0 1 2
-6 -2 -1 3
1 3
Starting again rom the 3 entry at (3, 3), this value came rom the 1 entry at (2, 2) in a diagonal move, extending the alignment to A-T ||| ATT
The dynamic matrix now looks like
0 -2 -4
-2 2 0
-4 0 1
3
1
3
Continuing in this ashion (try to do this), the nal alignment is GGA-T : : : , GAATT
where it is customary to represent a matching character with a colon ‘:’. The traceback path in the dynamic matrix is
0
3
1
2 1 3
I the other degenerate path was initially taken, the nal alignment would be GGAT: :: GAATT
111
7.4. GAPS
and the traceback path would be
0
5
3
2 1 3
The score o both alignments is easily recalculated to be the same, with 2 1 + 2 2 + 2 = 3 and 2 1 + 2 + 2 2 = 3. The algorithm or aligning two proteins is similar, except match and mismatch scores depend on the pair o aligning amino acids. With twenty diferent amino acids ound in proteins, the score is represented by a 20 20 substitution matrix . The most commonly used matrices are the PAM series and BLOSUM series o matrices, with BLOSUM62 the commonly used deault matrix.
−
−
−
−
×
7.4
Gap opening and gap extension penalties
Empirical evidence suggests that gaps cluster, in both nucleotide and protein sequences. Clustering is usually modeled by diferent penalties or gap opening ( ) and gap extension ( ), with < < 0. For example, the deault scoring scheme or the widely used BLASTN sotware is +1 or a nucleotide match, 3 or a nucleotide mismatch, 5 or a gap opening, and 2 or a gap extension. Having two types o gaps (opening and extension) complicates the dynamic programming algorithm. When an indel is added to an existing alignment, the scoring increment depends on whether the indel is a gap opening or a gap extension. For example, the extended alignment
−
AB || ab
to
−
−
AB||| abc
adds a gap opening penalty to the score, whereas A|| ab
to
A-||| abc
adds a gap extension penalty to the score. The score increment depends not only on the current aligning pair, but also on the previously aligned pair. The nal aligning pair o a sequence o length with a sequence o length can be one o three possibilities (top:bottom): (1) : ; (2) : ; (3) : . Only or (1) is the score increment ( , ) unambiguous. For (2) or (3), the score increment depends on the presence or absence o indels in the previously aligned characters. For instance, or alignments ending with : , the previously aligned character pair could be one o (i) −1 : , (ii) : , (iii) −1 : . I the previous aligned character pair was (i) or (ii), the score increment would be the gap opening penalty ; i it was (iii), the score increment would be the gap extension penalty . To remove the ambiguity that occurs with a single dynamic matrix, we need to compute three dynamic matrices simultaneously, with matrix elements
−
−
−
−
−
112
CHAPTER 7. SEQUENCE ALIGNMENT
denoted by (, ), − (, ) and − (, ), corresponding to the three types o aligning pairs. The recursion relations are (1) :
(, ) = max (2) :
−
⎧⎪⎨ ⎪⎩
( 1, 1) + ( , ), − ( 1, 1) + ( , ), − ( 1, 1) + ( , );
− − − − − −
− (, ) = max (3)
− : −
(, ) = max
⎧⎪⎨ ⎪⎩ ⎧⎪⎨ ⎪⎩
(7.4)
( 1, ) + , − ( 1, ) + , − ( 1, ) + ;
(7.5)
(, 1) + , − (, 1) + , − (, 1) + ;
(7.6)
− − −
− − −
To align a sequence o length with a sequence o length , the best alignment score is the maximum o the scores obtained rom the three dynamic matrices: opt (, ) = max
⎧⎪⎨ ⎪⎩
(, ), − (, ), − (, ).
(7.7)
The traceback algorithm to nd the best alignment proceeds as beore by starting with the matrix element corresponding to the best alignment score, opt (, ), and tracing back to the matrix element that determined this score. The optimum alignment is then built up rom last-to-rst as beore, but now switching may occur between the three dynamic matrices.
7.5
Local alignments
We have so ar discussed how to align two sequences over their entire length, called a global alignment . Oten, however, it is more useul to align two sequences over only part o their lengths, called a local alignment . In bioinormatics, the algorithm or global alignment is called “Needleman-Wunsch,” and that or local alignment “Smith-Waterman.” Local alignments are useul, or instance, when searching a long genome sequence or alignments to a short DNA segment. They are also useul when aligning two protein sequences since proteins can consist o multiple domains, and only a single domain may align. I or simplicity we consider a constant gap penalty , then a local alignment can be obtained using the rule
(, ) = max
⎧⎪⎪ ⎨ ⎪⎪⎩
0, ( 1, 1) + ( , ), ( 1, ) + , (, 1) + .
− −
−
−
(7.8)
