


Full Text  
PAGE 1 MATHEMATICAL MODELS IN POPULATION DYNAMICS BY ALEXANDER SALISBURY A Thesis Submitted to the Division of Natural Sciences New College of Florida i n partial fulfillment of the requirements for the degree Bachelor of Arts Under the sponsorship of Dr. Necmettin Yildirim Sarasota, FL April 2011 PAGE 2 ii ACKNOWLEDGEMENTS I would like to thank my advisor Dr. Necmettin Yildirim for his support, guidance, and seemingly unlimited supply of patience. Additional thanks to my thesis committee members Dr. Chris Hart and Dr Eirini Poimenidou for their guidance and criticism. Final thanks to family and friends for their love and support. PAGE 3 iii TABLE OF CONTENTS Acknowledgements .............................................................................................................................. ii Table of Contents .................................................................................................................................. iii List of Tables and Figures ................................................................................................................. vi Ab stract .....................................................................................................................................................1 Chapter 1: Background ........................................................................................................................2 1.1 What are Dynamical Systems? ............................................................................................................. 2 1.2 Formulating the Model ........................................................................................................................... 5 1.3 Methods for Analysis of Population Dynamics .............................................................................. 8 Solving Differential Equations ................................................................................................................ 8 Expressing in Dimensionless Form ...................................................................................................... 8 One Dimensional Models: Geometrical Analysis ......................................................................... 10 One Dimensional Models: Local Linearization ............................................................................. 11 Two Dimensional Models: Geometrical Analysis ........................................................................ 13 Two Dimensional Models: Local Linearization ............................................................................ 15 Classification of Equilibria .................................................................................................................... 19 1.4 An Historical Overview of Population Dynamics ....................................................................... 23 Fibonacci ...................................................................................................................................................... 25 Leonhard Euler .......................................................................................................................................... 26 Daniel Bernoul li ........................................................................................................................................ 27 Thomas Robert Malthus ......................................................................................................................... 28 PierreFranois Verhulst ....................................................................................................................... 29 Leland Ossian Howard and William Fuller Fiske ......................................................................... 32 Raymond Pearl .......................................................................................................................................... 33 PAGE 4 iv Alfred Jam es Lotka and Vito Volterra ............................................................................................... 36 Anderson Gray McKendrick and William Ogilvy Kermack ....................................................... 40 Georgy Frantsevich Gause ..................................................................................................................... 43 Chapter 2: Single Species Population Models ........................................................................... 45 2.1 Malthusian Exponential Growth Model ......................................................................................... 47 Analytic Solution ....................................................................................................................................... 47 Geometrical Analysis ............................................................................................................................... 48 Assumptions of the Model ..................................................................................................................... 50 2.2 Classical Logistic Growth Model ...................................................................................................... 51 Analytic Solution ....................................................................................................................................... 52 Obtaining Equilibrium Points .............................................................................................................. 53 Geometrical Analysis ............................................................................................................................... 54 Local Linearization .................................................................................................................................. 56 Assumptions of the Model ..................................................................................................................... 57 2.3 Theta Logistic Growth Model ............................................................................................................ 58 2.4 Logistic Model with Allee Effect ....................................................................................................... 61 Geometrical Analysis ............................................................................................................................... 63 2.5 Growth Model with Multiple Equilibria ........................................................................................ 65 Geometrical Analysis ............................................................................................................................... 66 Chapter 3: Multispecies Population Models .............................................................................. 68 3.1 Interspecific Competition Model ...................................................................................................... 71 Obtaining Equilibrium Points .............................................................................................................. 72 Geometrical Analysis ............................................................................................................................... 73 Local Li nearization .................................................................................................................................. 80 3.2 Facultative Mutualism Model ............................................................................................................ 82 PAGE 5 v Obtaining Equilibrium Points .............................................................................................................. 83 Geometrical Analysis ............................................................................................................................... 84 Local Linearization .................................................................................................................................. 86 3.3 Ob ligate Mutualism Model ................................................................................................................. 88 Geometrical Analysis ............................................................................................................................... 88 3.4 Predator Prey Model ............................................................................................................................ 92 Geometrical Analysis ............................................................................................................................... 93 Local Linearization .................................................................................................................................. 96 Chapter 4: Concluding Remarks .................................................................................................... 98 References .......................................................................................................................................... 101 PAGE 6 vi LIST OF TABLES AND FIGURES Figure 1.1 .........................................................................................................................................4 Figure 1.2 .........................................................................................................................................5 Figure 1.3 .......................................................................................................................................11 Figure 1.4 .......................................................................................................................................15 Figure 1.5 .......................................................................................................................................21 Figure 1.6 .......................................................................................................................................22 Table 1.1 ....................... ................................................................................................................31 Figure 1.7 .......................................................................................................................................35 Figure 1.8 .......................................................................................................................................40 Figure 2.1 .......................................................................................................................................48 Figure 2.2 .......................................................................................................................................48 Figure 2. 3 .......................................................................................................................................49 Table 2.1 .......................................................................................................................................53 Figure 2.4 .......................................................................................................................................54 Figure 2.5 .......................................................................................................................................54 Figure 2.6 .......................................................................................................................................59 Figure 2.7 .......................................................................................................................................59 Figure 2.8 .......................................................................................................................................60 Figure 2.9 .......................................................................................................................................63 Figure 2.10 .......................................................................................................................................63 Figure 2.11 .......................................................................................................................................66 Figure 2.12 .......................................................................................................................................66 Table 3.1 .......................................................................................................................................68 Figure 3.1 .......................................................................................................................................74 Table 3.2 .......................................................................................................................................74 Figure 3.2 .......................................................................................................................................75 Figure 3.3 .......................................................................................................................................75 Figure 3.4 .......................................................................................................................................76 Figure 3.5 .......................................................................................................................................77 Figure 3.6 .......................................................................................................................................7 8 Figure 3.7 .......................................................................................................................................84 Figure 3.8 ....................... ................................................................................................................85 Figure 3.9 .......................................................................................................................................89 Figure 3.10 .......................................................................................................................................90 Figure 3.11 .......................................................................................................................................94 Figure 3.12 .......................................................................................................................................9 5 PAGE 7 1 MATHEMATICAL MODELS IN POPULATION DYNAMICS Alexander Salisbury New College of Florida 2011 ABSTRACT Population dynamics studies the changes in size and composition of populations through time, as well as the biotic and abiotic factors influencing those changes. For the past few centuries, ordinary differential equations (ODEs) have served well as models of both singlespecies and multispecies population dynamics. In this study, we provide a mathematical f ramework for ODE model analysis and an outline of the historical context surrounding mathematical population modeling. Upon this foundation, we pursue a piecemeal construction of ODE models beginning with the simplest onedimensional models and working up in complexity into two dimensional systems. Each modeling step is complimented with mathematical analysis, thereby elucidating the models behaviors, and allowing for biological interpretations to be established. Dr. Necmettin Yildirim Division of Natural Sciences PAGE 8 2 CHAPTER 1 : BACKGROUND The aim of this section is to elaborate on basic concepts and terminology underlying the study of dynamical systems Here, we provide a basic review of the literature to date with the intent of foster ing a better understanding of concepts and analyses that are used in later sections We will begin an introduction to ordinary differential equation (ODE) models and methods of analysis that have been developed over the past several centuries followed by an historical overview of the field of dynamics Applications in population ecology will be of particular emphasis. 1.1 WHAT ARE DYNAMIC AL SYSTEMS? A system may be loosely defined as an assemblage of interacting or interdependent objects that collectively form an integrated whole. D ynamical systems describe the evolution of s ystems in time A dynamical system is said to have a state for every point in time, and the state is subject to an evolution rule which determines what future states may follow from the current or initial, state. Whether the system settles down to a state of equilibrium, becomes fixed into steadily oscillating cycles, or fluctuates chaotically, it is the systems dynamics that describe what is occurring (Strogatz, 1994) A system that appear s steady and stable is in fact, the result of forces acting in cohort to produce a balance of PAGE 9 3 tendencies I n certain instances, only a small perturbation is required to move the system into a completely different state. T his occurrence is called a bifurcation Systems of naturally occurring phenomena are generally constituted by discrete subsystems with their own sets of internal forces. T hus in order to avoid problems of intractable complex ity the system must be simplified via the obs ervers discretion. For instance, we might say that for the microbiologist, the system in question is the cell, and likewise, the organ for the physiologist, the population for the ecologist and so on An apt model thus requires a carefully selected set of variables chosen to represent the corresponding real world phenomenon under investigation. Detailing complex systems requires a language for precise description, and as it turns out mathematical mode ls serve well to describe the systems under consideration. Dynamical systems may be r epresented in a variety of ways. They are m ost commonly represented by continuous ordinary differential equations ( ODEs) or discrete difference equations Other manifestations are frequently found in partial differential equations (PDEs), lattice gas automata (LGA), cellular a utomata (CA), etc The focus of this work lies primarily on systems represented through ODEs The dynamic behavior of a system may be determined by inputs from the environment, but as is often the case, feedback from the system allows it to regulate its own dynamics internally Feedback loops are characterized as positive or negative A typical example of a positive feedback loop is demonstrated by s o called arms races whereby two sovereign powers escalate arms production in response to each other, leading to an explosion of uncontrolled output In contrast, negative feedback is exemplified by the typical household thermostat, whereby perturbations in temperature are regulated by the PAGE 10 4 thermostats response, which maintains temperature constancy by either sending a heated or cooled output. Thus, positive feedback tends to amplify perturbations to the system, or amplify the systems initial state whil e negative feedback tends to dampen disturbances to the system as time progresses As we shall see throughout this work feedback plays an important role in the stability of systems. A system is said to be at equilibrium if opposing forces in the system are balanced, and in turn the state of the system remains constant and unchanged A system is said to be stable if its state returns to a state of equilibrium following some perturbation ( e.g., an environmental disturbance) A system is globally stable if its state returns to equilibrium following a perturbation of any magnitude, whereas, a locally stable system indicates that displacements must occur in a defined neighborhood of the equilibrium in order for the system to retur n to the same state of equilibrium Figure 1. 1. R olling ball analogy for stable and unstable equilibria. The notion of stability is illustrated in Figure 1.1 by means of a ball resting atop a peak (unstable position) and i n the dip of a valley (stable position) Imagining a landscape with multiple peaks and valleys is, by analogy, to imagine a global landscape with multiple local points of stability (valleys) and of instability (peaks). T he peaks in the landscape define the thresholds separating each of the distinct equilibria, and therefore the level of perturbation that the system must undergo is analogous to that of the peaks magnitude. Systems and their stability are considered in greater depth in Section 1.3 PAGE 11 5 1.2 FORMULATING THE MODEL When we refer to dynamical systems, in fact, we are generally referring to an abstracted mathematical model, as opposed to the actual empirical phenomenon whose dynamics we are attempting to describe. We begin by attempting to identify the physical variables that we believe are responsible for the behavior of the phenomenon in question, and then we may formulate an equation, or system of equations, which also reflects the interrelation of our assumed variables. As depicted in Figure 1.2 the model building process involves the repetitive steps of observation, deduction, (re)formulation, and validation (Berryman & Kindlmann, 2008) Figure 1. 2 Flow diagram of general modeling process A basic aim of modeling is to help illuminate the mechanics that underlie some real world phenomenon, whether its nature is biological, chemical, physical, economic, or otherwise. Obtaining results that are consistent with the real world phenomenon is a necessary but not suff icient property of a good model, and as we shall see, there are several criteria by which an apt model should be uphe ld. The first step in formulating a model is to delineate the major factors governing the real world situation that is to be modeled (Berryman & Kindlmann, 2008) Conceptualizing the problem such that all key variables are accounted for, insofar as they reflect the mechanics of the observable phenomenon, is a good method for producing a testable PAGE 12 6 mode l. Initial sketches of a model may be done using a flow chart diagram or pseudocode, which illustrates state variables and the nature of their connections. G ilpin & Ayala ( 1973) propose the following criteria by which a good model should uphold : 1. Simplicity By virtue of Occams razor simple models are favorable ove r complicated models "because their empirical content is greater; and because they are better testable" (Popper, 1992) Incorporating the minimum possible number of parameters to account for the observed results is always favorable As Albert Einstein famo usly stated, Everything should be made as simple as possible, but not simpler. 2. Reality All of the models parameters should have biological relevance and attempt to reflect the mechanics of the biological system in question. The modeler should therefore hold a sound understanding of the real world phenomenon in question. Models that explain a phenomenon f rom first principles or from the bottom up are said to be mechanistic Th ey acknowledge that a biological phenomen on is the sum of multiple distinct, yet intertwined, processes, and therefore they attempt to describe the phenomenon in terms of its primary mechanisms ( in ecology, often at the level of the individual ). By contrast, models that describe a phenomenon are said to be phenomenological. The structure of a phenomenological model is empirically determined top down from a populations characteristics and therefore cannot predict behaviors independent of the original data. The parameters used in phenomenological models are therefore conglomerate sums of numerous lower level mechanisms ; (Schoener ) calls them megaparameters (1986). PAGE 13 7 Schoener provides a worthwhile summary of the mechanistic approach in eco logical modeling ( 1986) ultimately favoring it over the phenomenological approach by imagining a mechanistic ecologists utopia. However, both modeling approaches, m echanistic and phenomenological, have their advantages and disadvantages in different scenarios. Nearly all of the models we have chosen to consider herein are phenomenological because they offer a comparatively convenient mathematical form and flexibility in terms of analysis 3. Generality Using dimensionless variables allows magnitudes to take on a general significance, in turn providing scalability. Then, for specific scenarios the general m odel may take on specificity to account for the particular case. 4. Accuracy The model should vary from the observed data as little as possible. Hence a model w ith little to no predictive or explanatory utility should undergo further revision. Prior to embarking on the step by step procedures used to formulate and analyze continuous population models, we will consider population dynamics modeling from an historical perspective, providing insights into the key figures associated with the field of populati on ecology in addition to the methods they developed in order to understand population systems from a mathematical perspective. Additionally, we shall take this as an opportunity to introduce new terms and concepts. PAGE 14 8 1.3 METHODS FOR ANALYSIS OF POPULATI ON DYNAMICS There are multiple techniques employed for interpreting the behaviors of population dynamics models. However, not all c ontinuous models may be analyzed using the same toolset, and in many cases explicit solutions are impossible to achieve. We will be primarily considering two complimentary techniques of analysis: algebraic and geometric, which provide information regarding equilibria and their stability. All of the models we consider are based on ordinary differential equations, and thus, any person with a background in calculus should be capable of understanding the techniques covered. Solving Differential Equations Most continuous models of population dynamics are based on differential equations, which can be solved using a variety of techniques, which will in large part be omitted from this study Unfortunately, only the simplest of models are analytically solvable, leaving the necessity for other techniques of analysis for models with greater complexity. Examples of equations solved in a step by step fashion are in Chapter 2. In light of the fact that some models are too difficult to solve, or are simply unsolvable (e.g multispecies models discussed in Chapter 3), additional methods must be used in order to gain knowledge about the systems behaviors. Expressing in Dimensionless Form Several advantages are conferred by expressing a model in dimensionless or nondimensio nal terms. First, the units of measure are not important in calculations and in any case may be brought back into the model at the end of analysis. Recalling from Section 1.2, the criteria for simplicity and generality ; these features are upheld by express ing the PAGE 15 9 model in dimensionless terms without fear of any loss of generality. More importantly, reducing the number of relevant parameters into dimensionless groupings better illuminates the relationship s between parameters, while simultaneously allowing ca lculations to be made with greater ease. For example, c onsider Verhulsts logistic equation, which has a net growth rate parameter 0 r and a carrying capacity parameter 0 K 1, dNNrN dtK 0(0). NN ( 1 1 ) Here, we may introduce population and time in terms of dimensionless quantities, respectively, by N Q K and rt ( 1 2 ) Rewri ting the equation in the dimensionless terms, we obtain (1), dNdNdQd dtdQddt dQ Kr d rKQQ ( 1 3 ) Solving for Eq. (1.3) for dQ d and putting Eq. (1.1) back in, we get (1), dQ QQ d 0(0) QQ ( 1 4 ) where 00/ QNK is dimensionless. From this point, there are no parameters except 0 Q and the model can be solved using standard methods. PAGE 16 10 One Dimensional Models: Geometrical Analysis Qualitatively informed geometrical analysis of differential equations provides a visual representation of the systems dynamics, allowing one to gain a general insight into the behaviors of the system without the need to solve or compute. Determining a systems stability is an important yet easily achieved process for one dimensional systems Here we may recall the concepts of local versus global stability that were introduced in Section 1.1, but first let us define stability more precisely. Consider the following onedimensional differential equation: (), dN fN dt ( 1 5 ) where () fN is a continuously differentiable (typically nonlinear) function of N We say that *NN is an equilibrium point (also known as a fixed point steady state critical point or rest point ) where *0dN dt NN Equilibrium points can be calculated by solving *()0 fN That is to say, at equilibrium, there are no changes occurring in the system through time. It should be noted that there could be more than one value of *N that satisfies *()0 fN For instance, in addition to whatever equilibrium points a population N may reach (where *0N ), a trivial equilibrium generally found where*0, N indicating the biologically nontrivial fact that a population may not grow from a population of zero individuals. We can also note that if ()0, fN then N will increase, and if ()0, fN then N will decrease. By plotting the phase line of dN dt as a function of N it becomes an easy task to gain insight into the systems dynamics. Simply, the points of intersection at the N axis indicate that they are fixed points since *()0 fN at those points. PAGE 17 11 Equilibrium points are classified as either stable o r unstable In Figure 1.3, s table equilibrium points are represented graphically as filled in dots, and in stable equilibria perturbations dampen over time. By contrast, u nstable equilibrium points are represented as unfilled dots, and in unstable equilib r ia disturbances grow in time. Uns table equilibrium poi nts also may be referred to as sources or repellers and stable equilibrium points may be referred to as sinks or attractors An equilibrium point *N is asymptotically stable if all (sufficiently small) perturbations produce only small deviations that eventually return to the equilibrium. Suppose that *N is a fixed point and that () fN is a continuously differentiable function, and *'()0 fN Then the fixed point *N is considered asymptotically stable if *'()0 fN and asymptotically unstable if *'()0 fN Figure 1. 3 Phase line portrait of a population model () dN dt fN The trajectory has 3 non trivial equilibria 123,, NNN One Dimensional Models: Local Linearization The aforementioned geometrical techniques serve a utility by providing a means of intuitive analysis of equilibrium points that is qualitative in nature. A complimentary form PAGE 18 12 of steady state (equilibrium) analysis is achieved by linearizing the equation locally about the equilibrium points For this section, we will be following along the lines of (Kaplan & Glass, 1995) Reconsider Eq. (1.5) : (). dN fN dt We may recall that the equations equilibrium points are found by solving ()0 fN for N and such values of N are denoted *N Performing a Taylor series expansion of () fN for each equilibrium point *N in its neighborhood yields 2 *2 2 *1 ()()()(). 2NN NNdf df fNfNNN NN dN dN ( 1 6 ) In the neighborhood (i.e., within very close proximity) of *N all higher order terms such as *2() NN are insignificant compared to *() NN and therefore they are removed from the equation and *()0 fN yielding an approximated form of () fN given by the function *()().NNdf fNNN dN We may further simplify by defining two new variables df dN NNm and *xNN which gives *() dxd dN dtdt dtNN Therefore Eq. (1.6) becomes (). dx fNmx dt ( 1 7 ) This result is the linear equation for exponential growth or decay (described in further depth in Chapter 2). Thus, if 0 m then there is an exponential departure from the fixed PAGE 19 13 point, indicating that it is unstable By contrast, if 0 m then there will be an exponential convergence to the equilibrium point, indicating that it is stable Two Dimensional Models: Geometrical Analysis Prior qualitative analysis was limited to onedimensional models; here we will extend the case to include t wo dimensional systems. Consider the following two dimensional system of equations : 1 12(,), dN fNN dt ( 1 8 ) 2 12(,). dN gNN dt ( 1 9 ) The phase plane is the two dimensional phase space on which the systems trajectories are mapped, thus allowing certain behaviors to be visualized geometrically, and without the necessity for an analytic solution. The vector field or slope field of Eqs. (1.8) and (1.9) is plotted by choosing any arbitrary point 00 12(,) NN in the plane and substituting the point for 12(,) NN in the equations to obtain the slope at that point. Repeating this process at arbitrary but consistent intervals across the plane, while plotting each slope as a line segment achieves an approximate view of where the systems integral curves lie. These are unique parametric curves that lie tangent to the line segments. Finding the nullclines or zerogrowth isoclines of the system provides additional information on the systems dynamics. An isocline occurs where line segments in the vector field all have the same slope. Nullclines are a special case of isocline where the slope equals zero; thus, the 1 N nullcline is found when 12 (,)0 fNN and the 2N nullcline is found when 12 (,)0 gNN Their point of intersection marks the equilibrium point. PAGE 20 14 O ne type of equilibrium, illustrated by Figure 1.4 is called a center which behaves in a neutrally stable fas hion, much like the pathological frictionless pendulum, as May (2001) describes it. Here, a prey species and a predator species are represented by 1N and 2N respectively, while solid and dotted red lines represent their respective nullclines. We ca n observe that the equilibrium occurs at the point of intersection 12(,)(1,1) NN of both nullclines. ( T he nullclines in this case are straight lines; however they may take the shape of any curve.) Solutions travelling on the surrounding loops represent periodic oscillations, each of which remain stable on a closed trajectory. We will survey other classifications of equi libria in the following subsection. PAGE 21 15 Figure 1. 4 Phase portrait of Lotka Volterra predator prey system. Two Dimensional Models: Local Linearization The geometrical analysis considered previously indicates the existence and positioning of equilibria. A complimentary analysis involves linearizing about the systems equilibria to investigate the stability of each equilibrium point locally In the same vein as that which we performed on one dimensional systems, we will employ Taylors theorem to linearize the equations in the neighborhood of the equilibrium points, and determine the characteristics of the systems equilibria (Kaplan & Glass, 1995) Reconsider the system o f equations (1.8) and (1.9) : PAGE 22 16 1 12(,) dN fNN dt 2 12(,). dN gNN dt Phase trajectories are solutions of 112 212(,) (,) dNfNN dNgNN ( 1 10) So lutions to Eqs. (1.8) and (1.9) provide the parametric forms of the phase trajectories, where t is a parameter. A unique curve passes through any arbitrary point 00 12(,) NN except at equilibrium points ** 12(,) NN where **** 1212(,)(,)0. fNNgNN ( 1 11) Carrying out a Taylor series expansion of the nonlinear functions 12(,) fNN and 12(,) gNN in the neighborhood of the equilibrium point ** 12(,) NN gives us: ** ** 12 12 ** ** 12 12** 1212 11 22 12 ,, ** 1212 11 22 12 ,,(,)(,)()(), (,)(,)()(),NN NN NN NNff fNNfNN NN NN NNgg gNNgNN NN NN NN ( 1 12) where ellipses represent higher order (nonlinear) terms such as 2 *2 1 11 2 11 () 2 N NN N and 2 *2 2 22 2 21 (), 2 N NN N for f and g respectively. If we now let 11, XNN ( 1 13) 22,YNN ( 1 14) PAGE 23 17 t he n 1dN dX dtdt 2dN dY dtdt and ** **12 12(,)0(,) fNNgNN The original system defined by Eqs. (1.8) and (1.9) then b ecomes dX aXbY dt ( 1 15) dY cXdY dt ( 1 16) where ,,, abcd are given by the matrix ** 1212 12 ,. NNff NN ab cdgg NN A ( 1 17) This is the so called Jacobian matrix, which expresses via partial derivatives how each comp onent iN of the system changes with respect to itself and all other components. In the neighborhood of the equilibrium point, the higher order terms are negligible in comparison to the linear terms. Consequently, the nonlinear syste m can be ap proximated by a linear system dX aXbY dt ( 1 18) dY cXdY dt ( 1 19) The eigenvalues of the linearized equations describe the geometry of the vector field in the neighborhood of the equilibrium points. Let 1 and 2 be the eigenvalues of A : 2210 det () 0 01 ab adadbc cd 2 1,2()4 22 adbc ad ( 1 20) PAGE 24 18 We can observe how the eigenvalues of a matrix are also related to its determinant and trace: 12det() A ( 1 21) 12tr(). A ( 1 22) The eigenvalues of A can be found from the determinant and trace via: 2 1,2tr()tr()4det() 2 AAA ( 1 23) We can easily find the associated eigenvectors for each eigenvalue by solving : 11 22. NN ab NN cd ( 1 24) The eigenvalues reflect the rate of chang e of perturbations of size (0)i n that occur n ear the equilibrium point, and t hese rates of change occur along the eigenvectors passing through the equilibrium point. We can view the change in a perturbation as a function of the eigenval ues and eigenvectors at equilibrium: 1()(0) .iik tt ii intnecev ( 1 25) Thus, for di stinct eigenvalues, solutions are given by 121 1122 2,ttN cece N vv ( 1 26) where 1c and 2c are arbitrary constants, and 1v and 2v are the eigenvectors of A corresponding to their respective eigenvalues 1 and 2 If the eigenvalues are equal then 12()tccte describes the proportionality of the solutions. The eigenvectors are given by PAGE 25 19 21 (1),ii ip p v ,i ia p b 0, b 1,2. i ( 1 27) After the elimination of t p hase trajectories are mapped on the 12(,) NN plane. Several cases can be distinguished regarding the systems dynamics, depending on the state of the eigenvalues of matrix A in Eq. (1.17) Let us consider the resulting behaviors fo r various cases ( see Figure 1.6). Classification of Equilibria Focus. When the discriminant of Eq. (1.20) is negative, that is to say, when 2tr()4det()0 AA the eigenvalues are complex or imaginary This causes the trajectory to spiral around the equilib rium point. A focus can act as either a sink or source; i.e. its stability is classified as either stable or unstable. The imaginary part reveals how rapid the spiraling occurs, and the stability is reflected in the sign of the real part tr() 2 A If tr() 20 A then the focus is stable, and if tr() 20 A then the focus is unstable ( Figure 1.6). Center. A special limited case occurs when tr() 20 A namely that a neutrally stable trajectory remains on a closed path circling around the equilibrium point as in ( Figure 1.4). The resulting behavior is oscillations with a steady period. A perturbation of arbitrary magnitude would be required in order to move the trajectory onto a different closed path. PAGE 26 20 Node. When the discriminant 2tr()4det()0 AA and 2tr()tr()4det() AAA a node occurs. Here, both eigenvalues are real numbers with the same sign. If tr() 20 A then it is a stable node; if tr() 20 A then it is an unstable node. Additionally, w e may distinguish proper nodes from improper nodes For a stable proper node to occur, the eigenvalues must satisfy 120 and for an unstable proper node to occur, the eigenvalues must satisfy 120 A s table improper node has equal negative eigenvalues 120,, and an unstable improper node has equal positive eigenvalues 120. Saddle Point. When the discriminant 2tr()4det()0 AA and 2tr()tr()4det() AAA a saddle point occurs. Here, both eigenvalues are real, but have different signs ( Figure 1.6) The term saddle point originates from the fact that the trajectories behave in an analogous fashion as liquid poured onto a horses saddle; there is attraction towards the point center point, followed by a perpendicular repelling away from the point as the liquid repels off the sides of the s addle ( Figure 1.5). PAGE 27 21 Figure 1.5 A saddle point (blue dot) on the graph of 22zxy PAGE 28 22 Figure 1.6 Classification of equilibria and their associated eigenvalues. PAGE 29 23 1.4 AN HISTORICAL OVERVIEW OF POPULATION DYNAMI CS The study of dynamical systems has its origins in fifteenth century physics with Newton s invention of differential equations and solution to the tw o body problem; a two body problem is, for instance, to calculat e the motion of earth around the sun given the inversesquare law of gravi tational attraction between them (Strogatz, 1994) Newton and others in this time period (such as Euler, Leibniz, Gauss, and Laplace) worked to find analytic solutions to problems of planetary motion; yet, as it turned out, solutions to the threebody problem (e.g., sun, earth, and moon) were nearly impo ssible to achieve analytically in contrast to the two body problem ( 1994) A s a result, other approaches were developed. In the late 1800s Henri Poincar developed many of the graphical methods still used today for an alyzing the dynamics of systems that extend in complexity beyond the two body problem The geometrical approach pioneered by Poincar proved to be powerful approach to finding a global, qualitative understanding of a systems dynamics (Kaplan & Glass, 1995; Strogatz, 1994) While the geometrical approach to analyzing dynamical behaviors proved to be a powe rful method, there remained an additional source of analysis that could not be well harnessed until t he rise of computing in the 1950s With the tireless number crunching capabilities provided by the computer, numerical methods could finally realize a far greater potential. The computer allowed one to develop a more intuitive grasp of nonlinear equations by providing rapid numerical calculations This advancement in technology, coupled with the geometric methods of analysis, facilitated the surge of develop ments that PAGE 30 24 occurred in the field of nonlinear dynamics throughout the 19 60s and 1970s (Strogatz, 1994) For example, in 1963, Lorenz discovered the chaotic motion of a strange attractor. H e observed that the equations of his threedimensional system never settled down to an equilibrium state, but rather, they continued to oscillate in an aperiodic fashion (Lorenz, 1963) Additionally, running simulations from different, yet arbitrarily close, initial conditions led to unpredictably different behaviors. Plotted in 3 dimensions, the solutions to his equations fell onto a butterfly shaped set of points ( 1963) It was later shown that this set contained the properties of a fractal, and his example became a major influence in chaos theory (Strogatz, 1994) Today, the study of dynamics reaches far beyond applications in celestial mechanics, and it has achieved a truly interdisciplinary status. Significant roles have been es tablished for studying dynamical systems in biology, chemistry, physics, cognitive science, meteorology, the social sciences, finance, philosophy, and so forth Herein we will be considering dynamical systems solely from the standpoint of population ecology. The followi ng biographies are of key contributors to the study of population dynamics with a particular emphasis on individuals whose contributions and influences are most salient in the work outlined in the following chapters i.e., in continuous ordinary differential equation models of population dynamics. A key work used in outlining this section by Nicolas Baca r ( 2011) provides a compact yet thorough account of the historical figures associated with the development of population dynamics PAGE 31 25 Fibonacci Leonardo of Pisa, who posthumously became known as Fibonacci, finished writing Liber abaci in 1202, in which he explained various applications of the Arabic number system (decimal) in accounting, unit conversions interest rates, etc (Sigler, 2002) Appearing as a mere exercise in the midst of unrelated problems, Fibonacci outlined a problem that today wo uld be described as a problem in population dynamics {Document not in library: (Bacaer, 2011a)} He formulated his question with regard to a pair of mating rabbits and the number of offspring that could be expected after a given period of time. He wrote the following discrete difference equation: 11 ,nnnPPP ( 1 28) which states that the number of pairs of rabbits 1 nP after 1 n months is the sum of the number of pairs in month n and of the number of baby pairs in month 1 n ; however, baby rabbits cannot reproduce; therefore, they are considered to be the pairs that were present in month 1 n Fibonaccis rabbit problem was overlooked f or several centuries ; however, it is now recognized as one of the first models in population dynamics {Document not in library: (Bacaer, 2011a)} While the rab bit equation (1.28) turned out to be an unrealistic model (i.e., there are no limitations on growth, no mortality, etc.) the recurrence relation that bears Fibonaccis na me has an interesting relationship with naturally occurring geometries, and is found in numerous natural formations ranging from seashells to sunflowers {Document not in library: (Bacaer, 2011)} The ratio 1/nnPP approaches the so called PAGE 32 26 golden ratio 15 21.618 as n Despite the unrealistic nature of Fibonaccis model with regard to populations, it does share a common property with nearly all population models, namely geometrically increasing growth {Document not in library: (Bacaer, 2011)} Leonhard Euler Leonhard Euler was a Swiss mathematician born in 1701. He made numerous contributions in the fields of mechanics and mathematics and is considered to be one of the most prolific mathematicians of his time {Document not in library: (Bacaer, 2011a)} Given the breadth of his work and hi s display of interest in demograph y his work in population dynamics is only natural. Euler stated that a population nP in year n would satisfy the difference equation 1(1)nnPP ( 1 29) where n is a positive integer and the growth rate is a positive real number. With an initial condition 0P we find the population size in year n by the equation 0(1).n nPP ( 1 30) The form of growth assumed by this equation is called geometric growth (or exponential growth when dealing with continuous equations) As the son of a Protestant minister and having remained in strict religious faith Euler found this growth model to suit the biblical story in Genesis which held that the entire earths population descended from very few individuals, namely the three son s of Noah {Document not in library: (Bacaer, 2011)} Despite this ideological alignment, however, Euler recognized that the earth would never sustain such a high rate of growth g iven the fact that populations would have climbed upwards to 166 billion individuals in only 400 years Fifty years after Eulers PAGE 33 27 formulations Malthus considered the consequences of such growth with regard to human populations in his famous book titled An Essay on the Principle of Populatio n ( 1798) Daniel Bernoulli Daniel Bernoulli was born in 1700 into a family of already well established mathematicians: his father Johann Bernoulli and his uncle Jacob Bernoulli. His father did not want him to study mathematics, so Daniel began studying medicine, obtaining his doctorate in 1721 {Document not in library: (Bacaer, 2011a)} W ithin four years however, he published his first book on mathematics, titled Exercitationes quaedam mathematicae After his publication, he became involved in a series of professorships in botany, physiology, and physics, and around the year 1760, Bernoulli undertook studies analyzing the benefits of smallpox inoculation given the associated risk of death from inoculation. His model held the following assumptio ns: The number of susceptible individuals () St indicates those uninfected individuals at age t who remain susceptible to the smallpox virus The number of individuals () Rt indicates those whom are infected with the virus but who remain alive at age t The total number of individuals () Pt equals the sum of () St and () Rt The models parameters q and () mt respectively, represent each individuals probability of becoming infected with smallpox and each individuals probability of dying from other causes Given these assumptions, Bernoulli derived the following ODEs : (), dS qSmtS dt ( 1 31) PAGE 34 28 (1)(). dR qpSmtR dt ( 1 32) The sum of these equations yields (), dP pqSmtP dt ( 1 33) and using Eqs. (1.31) and (1.32) he yielded the fraction of susceptible individuals at age t by ()1 ()(1)qt St Ptpep ( 1 34) Bernoulli estimated the models parameters using Edmond Halleys life table, which provide d the distribution of living individuals for each age {Docume nt not in library: (Bacaer, 2011)} Choosing 1/8 q per year, and having eliminated () mt through mathematical trickery he computed the total number of susceptible people using Eq. (1.34) finding that approximately 1/13 of the populations deaths was exp ected to be due to smallpox. He further developed his model to examine the costs and benefits of inoculation, which he concluded were undoubtedly beneficial the life expectancy of an inoculated individual was raised by over three years Despite th ese findings the State never promoted smallpox inoculation and ironically the demise of King Louis XV in 1774 was a result of the smallpox virus {Document not in library: (Bacaer, 2011)} Thomas Robert Malthus Thomas Robert Malthus born 1766, was a British scholar who studied mathematics at Cambridge University, obtaining his diploma in 1791, and six years later becoming a priest of the Anglican Church {Document not in library: (Bacaer, 2011a)} PAGE 35 29 In 1798 Malthus anonymously published An Essay on the Principle of Population, as It Affects the Future Improvement of Society, With Remarks on the Speculations of Mr. Godwin, Mr. Condorcet and Other Writers ( 1798) In his book, he argued that the two named French authors optimistic views of an ever progressing society were flawed particularly they did not consider the rapid growth of human populations against the backdrop of limited resources (1798) For Malthus, the English Poor Laws, which favored population growth indirectly through subsidized feeding, did not actually help the poor, but to the contrary {Document not in library: (Bacaer, 2011a)} Given the growth of human populations proceeding at a far greater rate than the supply of food, Malthus predicted (albeit, incorrectly) a society plagued by misery and hunger The so called Malthusian growth model is described by the differential equation dN rN dt ( 1 35) where the growth of population N is governed by the net intrinsic growth rate parameter rbd which is the rate of fertility minus the mortality rate. Malthus emphasized that this equation holds true in capturing a growing populations dynamics only when growth goes unchecked (Malthus, 1798) However, the continued ex ponential growth of human populations against Earths limited resources Malthus argued, would ultimately lead to increased human suffering ( 1798) Malthus ideas proved to be influential in the work of numerous individuals, from Verhulsts density dependent growth model to ideas of natural selection pioneered by Charles Darwin and Alfred Russel Wallace {Document not in library: (Bacaer, 2011a)} Pierre Franois Verhulst PAGE 36 30 PierreFranois Verhulst was born in 1804 in Brussels, Belgium. A t the age of twenty one he obtained his doctorate in mathematics. While also bearing an interest in politics, Verhulst became a professor of mathematics in 1835 at the newly founded Free University of Brussels (Baca r, 2011) In 1835, Verhulsts mentor Adolphe Quetelet published A Treatise on Man and the Development of his Faculties ; he proposed that a populations long term growth is met with a resistance that is proportional to the square of the growth rate (Baca r, 2011) This idea encouraged Verhulsts developments found in Note on the law of population growth (1838 as cited in Baca r, 2011 ), in which he stated The v irtual increase of p opulation is limited by the size and the fertility of the country. As a result the population gets closer and closer to a steady state. Verhulst proposed the differential equation 1, dNNrN dtK ( 1 36) where the growth of the population N is governed by the Malthusian parameter r and the carrying capacity K (although at the time, these parameters were not named as such). The growth rate r decreases linearly against an increasing population densit y N However, w hen () Nt is small compared to K the equation can be approximated by the Malthusian growth equation dN rN dt ( 1 37) PAGE 37 31 which has the solution 0()rtNtNe where 0N is the initial number of individuals in the population However, we may also find the solution to Verhulsts logistic equation in Eq. (1.36) given by 0 0() 1(1)/rt rtNe Nt NeK ( 1 38) which describes the growth of population N increasing from an initial condition 0(0) NN to the limit, or carrying capacity, K which is reached as t Using available demographic data for various regions, Verhulst estimated parameters r and K using as few as three different but equally spaced years provided through census data. He showed that if the population is 0N at 0 t 1N at tT and 2N at 2 tT then both parameters can be estimated starting from Eq. (1.38) giving 011202 1 2 1022 NNNNNN KN NNN ( 1 39) 0 11/1/ 1 log 1/1/ NK r TNK ( 1 40) PAGE 38 32 Table 1.1 United States census data (1790 and 1840). Adapted from (Verhulst, 1845) Verhulsts work involving the logistic equation was overlooked for several decades; however, in 1922 biologist Raymond Pearl took notice of his work after rediscovering the same equation (Pearl & Reed, 1920) In following centuries the logistic equation proved to become highly influential; for instance, it is from the logistic models parameters r and K that r/K selection theory, pioneered by Robert MacArthur and E. O. Wilson ( 1967) took its name. Leland Ossian Howard and William Fuller Fiske Leland Ossian Howard was an American entomologist who served as Chief of Bureau of Entomology for the U nited S tates D epartment of A griculture (18941927) and W. F. Fiske headed T he Gypsy Moth Project in Massachusetts (19051911). Howar d had PAGE 39 33 been conducting research in Europe, and eventually arranged for parasites to be imported to the U.S. as agents of biological (pest) control. As experts in the rising field of biological control, collaboration between the two individuals ensued, resulting in a new set of concepts that had been ov erlooked prior, namely population regulatio n via functional relationships They proposed the terms facultative and catastrophic mortality, which respectively indicate different functional relationships between growth rate r and the population density (Howard & Fiske, 1911) Catastrophic mortality indicates a constant proportion of death in the population, regardless of density; the more familiar term for it now is den sity independence Facultative mortality indicates an increase of death in a population that is increasing in density, and is now more commonly referred to as density dependen ce Raymond Pearl Raymond Pearl was born in Farmington, New Hampshire in 1879. After obtaining his A.B. from Dartmouth in 1899, he studied at the University of Michigan, completing his doctorate in 1902 (Jennings, 1942; Pearl, 1999) During a brief stay in Europe, Pearl studied under Karl Pearson at University College, London, where he adopted a statistical view of biological systems, and eventually, after moving to Baltimore in 1918 to become professor of biometry at the Johns Hopkins University Pearl also became chief statistician at the Johns Hop kins Hospital (Jennings, 1942) While studying populations of Drosophila h e collected life expectancies, death rates, and so forth, and began discovering survivorship curves which turned out to be quite reminiscent of Verhulsts logistic curves (1 942) While Verhulsts logistic model ( Verhulst, 1845) would appear to be the precursor to PAGE 40 34 Pearls findings, it was in fact T. Brailsford Robertsons sigmoidal shaped chemical au tocatalytic curve that sparked Pearls insight ( Pearl, 1999) After s howing the consistency of the survivorship curves from organisms with varying life histories Pearl touted his finding as some law of population growth which in turn sparked a considerable controversy (Kingsland, 1985) In a paper co published with his associate L owell J. Reed, Pearl defended the logistic equation in its capacity to describe the growth of populations that should eventually reach a carrying capacity ( Lowell & Reed, 1920) : In a new and thinly populated country the population already existing there, being impressed with the apparently boundless opportunities, tends to rep roduce freely, to urge friends to come from older countries, and by the example of their well being, actual or potential, to induce strangers to immigrate. As the population becomes more dense and passes into a phase where the still unutilized potentialiti es of subsistence, measured in terms of population, are measurably smaller than those which have already been utilized, all of these forces tending to the increase of population will become reduce. While Robertsons sigmoidal autocatalytic curves were to o symmetrical to fit Pearl and Reeds data, they made adjustments to a ccommodate a more realistic fit, resulting in (), 1 atK Nt be ( 1 41) where population N has constant parameters ,, baK and, further, forming a generalized equation () 1tK Nt be ( 1 42) where 11 12 n naatat The number of terms and values of constants therefore determine the sigmoid curves precise shape. PAGE 41 35 The logistic equation, as opposed to the curve, is most commonly written in the ODE form shown in Eq. (1.36) ; however, it s curve can be written () 1 artK Nt e ( 1 43) where population N is marked by the maximum rate of growth m rr and parameters a and K represent the constant of integration and the carrying capacity, respectively. PAGE 42 36 Figure 1.7 Logistic growth of yeast population over time Data ( green dots ) collected from (Carlson, 1913 as cited in Raymond Pearl, Miner, & Parker, 1927) Logistic growth curve ( blue line ) fitted to data points with parameters 664.3 K 4.7 a and 0.536mr Figure 1.7 provides a visualization of the data collected from a yeast population with the corresponding fitted logistic growth curve (Carlson, 1913 as cited in Raymond Pearl, Miner, & Parker, 1927) Here, parameters were calibrated to 664.3 K 4.7 a and 0.536mr Pearl also fitted logistic curves to census data of several countries including France, Sweden, and the United States (1999) Alfred James Lotka and Vito Volterra Alfred James Lotka was born of American parents in the part of the Austro Hungarian Empire that is now Lviv, Ukraine. He studied physics and chemistry, receiving his bachelors from University of Birmingham in England, and eventually began work in PAGE 43 37 New Y ork for the General Chemical Company (Baca r, 2011) Despite his status as a physical chemist, his work has helped revolutionize the field of population ecology. Remaining unaware of Eulers work on the subject over a century prior, Lotka began studying the dynamics of agestructured populations, first marked by the publication of Relation between Birth Rates and Death Rates (Lotka, 1907 as cited in Baca r, 2011) His work follows a different approach than that of Euler, in that he uses continuous rather than discrete variables to represent age and time. Lotkas model is largely responsible for what has become known as stable population theory despite Euler having reached a similar result with his discrete model (Baca r, 2011) What is meant by stable population is that the populations age pyramid, that is to say, the distribution of ages within the population, remains stable regardless of the populations growth or decline ( 2011) Lotkas prior work involving oscillations in chemical dynam ics, along with his interest in the mathematics of ecological properties naturally led to his investigation of rhythms in ecological systems. In 1920, he published Analytical Note on Certain Rhythmic Relations in Organic Systems wherein he arrived at a system of equations used to describe the continuously oscillating dynamics of two populations: predators ( e.g., herbivores) and prey ( e.g., plants) where 1X and 2X represent the state of each species 1S and 2S respectively for all points in time 0 t (Lotka, 1920) He described the dynamics of the system verbally as: 1 1 11 2Other dead matter Mass of destroyed Rate of change ofMass of newly formed eliminated from per unit of time per unit of time by per unit of time per unit of time S S XS S and PAGE 44 38 2 2 2 1Mass of newly formed Rate of change of Mass of destroyed per unit time (derived per unit of time per unit of time from as injested food) S S X S Lotkas system of equations written in mathematical terms is thus : 1 1121 12(), dX XXXX dt XX ( 1 44) where and 2 122 21() dX XXX dt XX ( 1 45) where parameters ,,,, are functions of 1X and 2X (Lotka, 1920) After Lotkas publication ( 1920) was completed Raymond Pearl helped him obtain a scholarship from Johns Hopkins University, where Lotka was able to write his book in 1925, titled Elements of Physical Biology ( Lotka, 1925 as cited in Baca r, 2011) At the time, Lotkas book did not garner much attention, and it was not until Lotkas colleague Vito Volterra, who was a notable mathematical physicist, discovered the same equations that they earned their renowned status among ecologists (2011) Vito Volterra was born in a Jewish ghetto in Ancona, Italy, although at the time the city belonged to the Papal States. While remaining poor, Volterra performed well in school, completing his doctorate in physics in 1882 and subsequently obtaining a professorship at th e University of Pisa (Baca r, 2011) Volterra receiv ed considerable attention for his work in mathematical physics, and at the age of 65 he began investigating an ecolog ical problem proposed to him by his future son inlaw the zoologist Umberto dAcona. Volterra began PAGE 45 39 investigating the data collected between the years 1905 and 1923 on the varying proportions of sharks and rays landed in fishery catches in the Adriatic Sea {One or more documents not in library: (Bacaer, 2011b; Murray, 2002a)} DAcona had observed an increase in the populations of these predatory fishes during Wo rld War I, when harvesting activity was relatively reduced. Volterra unknowingly created the same mathematical model as Lotkas equations (1.44) and (1.45) to describe the dynamics of predator (shark) and prey (smaller fish) population The Lotka Volterra equations are standardly written as: 1 112, dN aNbNN dt ( 1 46) 2 212, dN cNdNN dt ( 1 47) where parameters ,,,0 abcd Coefficient a is the growth rate of prey in the absence of predators and c is the rate of decrease of the population of predators due to starvation (i.e., in the absence of prey ). The interaction terms 12bNN and 12dNN express the rates of mass transfer from prey to predators where db Lotka noticed that both populations satisfy the conditions for equilibrium in two scenarios. First, the so called trivial equilibrium occurs when ** 120, NN ( 1 48) where asterisks denote equilibrium. Here the prey population 1N is extinct, and likewise, the predator population 2N having no food source is extinct. Second, b oth populations coexist at nontrivial equilibrium when 1, c N d ( 1 49) PAGE 46 40 2. a N b ( 1 50) However, when both populations are not at eq uilibrium, then both functions 1(), Nt 2 () Nt be have in an oscillatory fashion with a period 0 T such that 11 ()() NtTNt and 22()() NtTNt for all 0 t For instance, if there is an abundant mass of prey in population 1 N then the predator population 2N will increase, in turn causing a decrease in 1 N When 1 N becomes too diminished to sustain feeding by 2N starvation occurs, causing 2N to decrease, and in turn, the mass of 1 N becomes rejuvenated. The process repeats itself indefinitely resulting in temporally offset oscillations of both populations (Baca r, 2011; E delsteinKeshet, 2005; Kot, 2001; Murray, 2002) These equations and their counterparts are elucidated in Chapter 3. Anderson Gray McKendrick and William Ogilvy Kermack Anderson Gray McKendrick was born in Edinburgh in 1876. He studied medicine at University of Glasgow before venturing abroad to fight diseases ( namely, malaria, dysentery, and r abies) in Sierra Leone and India (Baca r, 2011) He returned to E dinburgh in 1920 after contracting a tropical illness and began serving as superintendent of the Royal College of Physicians Laborator y. There, McKendrick met William Ogilvy Kermack, who served as head of t he chemistry division in the laboratory and with whom McKendrick would eventually begin collaborating (2011) In 1926, McKendrick published a paper titled Applications of mathematics to medical problems, in which he introduced continuous time models of epidemics with probabilistic effects determining infection and recovery (McKendrick, 1926 as cited in PAGE 47 41 Baca r, 2011) The paper served as the starting point for the famous S I R epidemic model, which was not fully developed until McKendrick and Kermack bega n collaborating ( Kermack & McKendrick, 1927) The S I R model derives its name from the progress ion of disease that individuals proceed through: susceptible ( S ), infected ( I), and recovered/resistant ( R) Figure 1.8 Kermack McKendricks S I R model Three possible states of progression: susceptible ( S ), infected ( I ), and recovered ( R ). A simplified form of the model demonstrating these disease dynamics follows as a threedimensional system of equations : Susceptible (): dS SSI dt ( 1 51) Infected (): dI ISII dt ( 1 52) Recovered/resistant (): dR RI dt ( 1 53) where parameters and respectively represent the rate of contact /infection and the rate of recovery ( which is proportional to the value of infected population I ) We can see that the quantity of new individuals belonging to the infected population I per unit time is proportional to the quantity of susceptible individuals and infected individuals while those PAGE 48 42 individuals in the susceptible population S are removed from S as they become infected () SI or later on, resistant () IR T he total population ()()()() NtStItRt must begin with a set of initial conditions (since the model assumes there is no birth, death, or migration) where (),(),() StItRt are 0 Thus, at the beginning of the epidemic ( 0 t ), the initial total population of size N contains a proper subset of infected individuals 0(0) II and susceptible individuals 00 (0) SSNI and we assume 0(0)0 RR since time must pass in order for individual s to pass from the infected state to the recovered state. There is no analytic solution t o this system; however, Kermack and McKendrick analyzed the properties of the system by other means. T hey observed that as t () St decreases to a limit 0 S while ()0 It and () Rt increases to a limit RN The equation log(), (0) S NS S ( 1 54) implicitly provides S and thus the final epidemic size may be obtained through RNS (Baca r, 2011; Kermack & McKendrick, 1927) T he S I R model thus provides the important biological indication that an epidemic ends before all susceptible individuals become infected (2011; 1927) Kermack and McKendrick continued developing disease models throughout the 1930s and their work has bec o me foundational in todays more complex epidemiological models (Baca r, 2011; Kingsland, 1985) PAGE 49 43 Georgy Frantsevich Gause Georgy F rantsevich Gause was a Russian biologist born in Moscow in 1910. He began his undergraduate studies at Moscow University under advisor W. W. Alpatov who was a student and friend of R aymond Pearl. Professor Alpatov may be credited for the direction of Gause s work, specifically in experimental population ecology (Kingsland, 1985) Foregoing field studies in favor of the controlled laboratory environment, Gause was able to control for potentially confounding variables in a series of ecological experiments performed in vitro In one experiment two competing species of Paramecium disp layed typical logistic growth when grown in isolation; h owever, when placed together in vitro one species always drove the other to extinction (Gause, 1934) By shifting environmental resource parameters (e.g., food and water), Gause found that the winner and loser species were not somehow predestined but rather dependent on the values of the resource parameters Similar results were yielded in experiments between two competing species of Saccharomyces yeast (Gause, 1932) These findings led to what has been called the principle of competitive exclusion or Gauses principle; stated briefly : complete competitors cannot coexist (Hardin, 1960) Restated if two sympatric nonbreeding populations (i.e., separate species occupying the same s pace) occupy the same ecological niche, and Species 1 has even an infinitesimally slightest advantage over Species 2, then Species 1 will eventually overtake Species 2, leading towards either extinction or towards an evolutionary shift to a different ecological niche for the less adapted species. As Hardin's First Law of Ecology states "You cannot do only one thing PAGE 50 44 Additionally, Gauses results were a confirmation of Lotka Volterras competition model, which by design, upholds the exclusion princ iple by assuming normal logistic growth for each species grown in isolation, and for species placed together, the mutual ly inhibitory effect of each species population (given the appropriate competition parameters) leads to the eventual demise of the dis advantaged population (Robert MacArthur & Levins, 1967) The phenomenon is als o called limiting similarity (1967) The Lotka Volterra competition model is described by the coupled system of equations: 1122 1 11 11, NN dN rN dt K ( 1 55) 2211 2 22 21, NN dN rN dt K ( 1 56) where populations of Species 1 and Species 2 are represented by 1N and 2N respectively This model and its counterparts are considered in further detail in Section 3.1 PAGE 51 45 CHAPTER 2 : SINGLE SPECIES POPULATION M ODELS The dynamics of single populations are g enerally described in terms of one dimensional differential equations. In this chapter, we consider onedimensional population models that were developed over the past several centuries to describe the growth and/or decay of single homogeneous populations. There is a pedagogical aim here in asserting the simple and fundamental principles at work in most continuous population models, and in this vein, we will begin from very simple foundations. Following a progression similar to that taken in EdelsteinKeshets text Mathematical Models in Biology ( 2005) we aim to elucidate the development of various models by augmenting in gradual increments The addition of each new parameter will be accompanied by an empirical and/or theoretical justification that will, in any case, provide the reader with a gradual (as opposed to what one might call saltationist) sense of the models evolution. It should be noted that, in most scenarios, the models outlined herein are inaccurate and oversimplified. They do not consider stochasticity (chance events), environmental effects spatial heterogeneity, or agestructure W ith regard to stochasticity, it is assumed that the deterministic model will produce result s that, on average, would be produced by the analogous stochastic model (Maynard Smith, 1974) The absence of c ertain realistic features does not negate the importance of these model s or the principles they convey Whereas the illustrative power of the models PAGE 52 46 outlined below is pedagogical in nature; therefore, the explanatory power of the underlying principles in large part, do minates the striving we might otherwise have for realism or accuracy. PAGE 53 47 2.1 MALTHUSIAN EXPONENTI AL GROWTH MODEL Recalling from Section 1.4 Thomas Robert Malt hus proposed that a populations growth will proceed exponentially if growth goes unchecked ( 1798) T he Malthus equation is denoted dN rN dt ( 2 1 ) where N represents the number o f individuals in the population (or more precisely, the biomass of the population) and r is a constant representing the intrinsic rate of growth. The growth rate r is also called the Malthusian parameter or the net intrinsic growth rate (i.e. rbd where b and d are intrinsic birth and death rates, respectively) Units of time t vary depending on the organism of inquiry. For instance, for rapidly multiplying organisms (e.g., bacteria), t may be measured in minutes, whereas for slowly multiplying organisms (e.g., elephants ), t may be measured in years. Analytic Solution The solution to Eq. (2.1) is easily achieved by separation of variables and integrating both sides of the equation, assuming that (0)(0) NNt to yield () (0) 0 NT tT NtdN rdt N ( 2 2 ) () 0 (0)ln ,NT T NNrtc ( 2 3 ) and evaluating th e upper and lower limits yields PAGE 54 48 ln()ln(0)()(*0), ln()ln(0) () (0)rTNTNrTcrc NTNrT NT e N ( 2 4 ) where () NT and (0) N are both positive. Rearranging the equation, we get the exact solution 0(),rtNtNe ( 2 5 ) w here the initial condition 0(0) NN Geometrical Analysis Malthusian growth described by Eq. (2.1) can manifest as both exponential growth and exponential decay (Figure 2.1) ; for instance, exponential growth occurs for all 0 r ( Figure 2.2); however, reversing the sign of r the model becomes one in which a population decays exponentially in time as the fraction r of individuals is removed per unit time (F igure 2. 3 ). Viewing the phase line ( Figure 2.1), the linear rate of change in population size, or density, is portrayed for both exponential growth and exponential decay. The only equilibrium solution 0dN dt occurs when *0N Qualitatively, we can judge the equilibrium points stability based on whether trajectories approach for all 0 r or zero for all 0 r PAGE 55 49 Figure 2.1 Phase line portrait of exponentia l model given by Eq. (2.1) : phase trajectory reveals the linear rate of change in growth as a function of N Figure 2.2 Dynamics of exponential growth given by Eq. (2.1) : exponential growth for a set of arbitrary positive growth rates r PAGE 56 50 Figure 2 .3. Dynamics of exponential growth described by Eq. (2.1) : exponential decay for a set of arbitrary negative growth rates r Assumptions of the Model The Malthus model is one of the simplest models of growth for any reproducing population; however, it is too simple to be useful in most circumstances As such, it makes the following assumptions: the population is homogeneous ( i.e. all members are identical); the population inhabits a uniform and unvarying environment; an infinite supply of nutrients is available; there are no spatial limit ations and growth is density independent. R ealistically, popula tion growth is limited by various factors from resource availability to predation. Additional limitations arise from the systems internal dynamics, such as overcrowding. The Malthus model may accurately describe the growth of a population for a limited period of time; however, unrestrained growth is never sustainable, and thus additional components are necessary to obtain a more realistic model. PAGE 57 51 2.2 CLASSICAL LOGISTIC GROWTH MODEL The logistic equation, developed by Verhulst (Sectio n 1.4) anticipates a l imit, or carrying capacity on population growth. This carrying capacity is symbolically represented K P lotting the populations growth as a function of time shows N approaching K along a sigmoid ( S shaped ) curve when the population s initial state 0N is below 2 K ; above 2 K solutions exponentially converge towards K (Figure 2.5) The addition of the new term K to our model is an intuitive advance from the Malthusian model since we know realistically that individuals cannot propagate infinitely in a finite space, and that the growth rate should decline as population density increases. The classic logistic model assumes t hat the individual growth rate ( ar ) is a linearly decreasing function of N such that ().arfN ( 2 6 ) We define mr as the maximum growth rate, which should decrease linearly as N increases. When NK the rate of growth will be zero, and growth rate will be come negative in the case of NK This new linearly decreasing growth rate is developed starting with an equation for a straight line yaxb where a represents the slope and b is the y intercept (here mr ). We calculate the slope by 21 210 0mmrr yy a xxKK ( 2 7 ) a nd the relationship between ar and N is PAGE 58 52 1.amN rr K ( 2 8 ) Now we may substitute ar for r in the original equation given in Eq. (2.1) to yield the logistic equation: 1.mdNN rN dtK ( 2 9 ) Analytic Solution The solution to Eq. (2.9) may be achieved via separation of variables and integrating both sides assuming (0)(0) NtN to yield 1N KdN rdt N ( 2 10) () (0) 0. 1NT tT N Nt KdN rdt N ( 2 11) Integrating requires the use of partial fractions: 1 (1)(1) 11, 1,NN KKAB NN N ABN K AN ABN K ( 2 12) where 1 A and 0A KB and further, we get () () (0) (0) 01 (1)NT NT T N NN KA B KK dNdN rdtNK ( 2 13) Integration and exponentiation on both sides yields PAGE 59 53 () (0)()(0) ln()ln(0)ln1ln1 expln()ln(0)ln1ln1 ,NT N rTT KKNTN NTN rT KK NTN e ( 2 14) and f urther simplifying, we get (0) () (0) ()expln()expln1 expln(0)expln1 ()(1) (0)(1)N K rt NT K N rt K NT KNT e N NT e N ( 2 15) F inally, solving for () NT and simplifying further provides us with the solution: 0 00() ,()rtNK NT NKNe ( 2 16) where the initial condition 0(0). NN Obtaining Equilibrium P oint s W e obtain the systems equilibrium points *N by finding all values of N that satisfy 0dN dt : *0 0mdN dt rN 10, or 10,mNrN K N K ( 2 17) and we get *0, N ( 2 18) *.NK ( 2 19) Thus, the logistic equation has exactly two equilibrium points. PAGE 60 54 N d N /d t NK 0 dN dt 0 NK 0 dN dt NK 0 dN dt 0 N 0 dN dt Table 2.1 Behavior logistic growth, described by Eq. (2.9) for different cases of N Geometrical Analysis Viewing the phase line in Figure 2.1 we can discern a number of facts co ncerning the systems dynamics; in fact, we see that the equilibria are already obtained graphically. We can also observe that any point 0N on the trajectory will approach K as t with the exception of the case 0 N in which there is no population. Table 2.1 illustrates the sign of dN dt for values of N The trivial equilibrium point *0 N is unstable, and the second equilibrium point *NK represents the stable equilibrium, where N asymptotically approaches the carrying capacity K In terms of the limit, we can say lim(), (0)0TNTKN A point of inflection occurs at 2 KN for all solutions that cross it, and we can see graphically that growth of N is rapid until it passes the inflection point 2 KN From there, subsequent growth slows as N asymptotically approaches K As shown in Table 2.1 if NK then 0dN dt and N decreases exponentially towards K This case should only occur when the initial condition 0(0) NNK In the PAGE 61 55 following section we will confirm the stability of equilibria by linearization about each equilibrium solution. Figure 2.4 Phase line portrait of logistic growth as described by Eq. (2.9) Figure 2. 5 Dynamics of the logistic model given by Eq. (2.9) PAGE 62 56 Local Linearization For a more quantitative measure of the systems stability, we may linearize the equation in the neighborhood of its equilibrium points. Let *()() ntNtN where () nt is a small perturbation in the neighborhood of an equilibrium point denoted *N We are interested in whether the perturbation grows or decays, so consider **()()(). dnd dN NNfNfNn dtdt dt ( 2 20) Performing a Taylor series expansion on Eq. (2.20) yields ** *()() ,NNdf fNnfNn dN ( 2 21) where ellipsis denotes quadratically small nonlinear terms in n that we will henceforth ignore. We may also eliminate the term *() fN since it is equal to zero, and we are provided with the approximated equation (). NNdf fNnn dN ( 2 22) Thus, 2()1 ,m mmrN N fNrNrN KK ( 2 23) 2 () .m mrN dfN r dNK ( 2 24) Hence, near the equilibrium points *0 N and *NK we obtain 02 ,m mm NrN dnnr rn dtK ( 2 25) PAGE 63 57 2 (2).m m mmm NKrN dn nr nrrrn dtK ( 2 26) Recalling the Malthus equation from previously, we see that dn m dt rn takes its form. Since 0mr t hese results indicate that the equilibrium point *0 N is unstable since the perturbation () nt grows exponentially if *'()0 fN On the other hand, *NK is stable since () nt decays exponentially if *'()0 fN Additionally, the magnitude of *'() fN tells how rapidly exponential growth or decay will occur, and its reciprocal 1 *'() fN is called the characteristic time scale which gives the amount of time it takes for () Nt to vary significantly in the neighborhood of *N (Strogatz, 1994) In this case, the characteristic time scale is 1 *1 '()mfNr for both equilibrium points. Assumptions of the Model The assumptions of the logistic model are the same as those of the Malthusian model, except that the reproduction rate is positively proporti onal to the size of the population when the population size is small and negatively proportional when the population is large. The point K towards which the population converges is the carrying capacity, and the parameter K is d etermined phenomenologically. As Sewall Wright cautioned, any flexible mathematical formula resulting in a sigmoid shape could be made to fit the data (Kingsland, 1985) Thus it would not be difficult to produce a curve fitted to the data by producing an algebraic expression and simply deriving a differential equation from which it is the solution (Murray, 2002) PAGE 64 58 2.3 THETA LOGISTIC GROWTH MODEL A simple variation on the classic logistic model incorporates a new term which provides additional generality and flexibility in term s of the change in per capita growth rate ar with respect to population density N With this augmentation, the model may yield more fitting result s under circumstances where density dependence is of importance. As su ch, the model provides additional complexity over the classic logistic model in terms of the shape of its growth curve (Figures 2. 6 and 2 .7 ) The updated per capita growth rate parameter becomes 1,amN rr K ( 2 27) where 0 Note that z ero growth would be given by 0 and in cases where 0 growth under K would decay to 0 while growth above K would proceed unbounded. The respective convexity or concavity of the curve is determined by whether 1 or 1 (Figure 2.6). S ubstituting the updated per capita growth rate ar into the original logistic equation yield s the thetalogistic growth model : 1.mdN N rN dt K ( 2 28) Varying the parameter is intended to reflect relation between intraspecific competition and population density. As such, the linear density dependence held by the classical logistic model can be altered to become curvilinear ( Figure 2.6). This is clear because PAGE 65 59 (), 1 (), 1 (), 1.aa aa aarr NN N rr N N r N r ( 2 29) That is to say, if we set 1 the carrying capacity term will be given more weight, in turn, weakening the density dependence for low values of N and thus reflecting scenarios in which crowding holds a lesser prominence ( since crowding has a lesser e ffect at low densities). When 1 the model is i dentical to the classic logistic model so its density dependence is a linearly decreasing function of N When 1 density dependence is strengthened for low values of N leading to slowed population growth (Figures 2.6 and 2. 7 ) Mechanistically, the value of should depend on the functional relationships between individuals at varying densities; however the parameter is phenomenological and therefore, does not possess a mechanistic significance per se Analyzing timeseries data from ~3200 different populations of insects, birds, mammals, and fish using the Global Population Dynamics Database, Sibly, et al. ( 2005) fou nd respective theta values for each population using a leastsquares approach In the majority of cases (~75%), populations displayed a concaveup density functional relationship (i.e., 1 ), indicating that many animals spend the majority of their time at or above carrying capacity ( 2005) PAGE 66 60 Figure 2.6 Plot of per capita growth rate ar described by Eq. (2.27) as a function of population density N for arbitrary values of 0.5, 1, 3, 10 Figure 2.7 Phase line portrait of theta logistic model described by Eq. (2.28) for arbitrary v alues of 0.5, 1, 3, 10 PAGE 67 61 2.4 LOGISTIC MODEL WITH ALLEE EFFECT We now consider an elaborated derivation of the logistic model intended to describe the situations in which a sparsity of individuals leads in turn to the reduced survival of offspring Restated in biological terms, the Allee e ffect is said to occur when dwindling population levels lead, in turn, to increasingly diminished reproduction despite the lack of intraspecific competition (Figure 2.8) In fisheries science literature, the effect is often called depensation (Courchamp, Clutton Brock, & Grenfell, 1999) Figure 2.8 Schematic plot of (a) negat ive (classical logistic type) and (b) positive (Allee effect) relationships between individual fitness and population density. Note Adapted from (Courchamp, Lud k & Gascoigne, 2008) Allee effects can be explained by several mechanisms, including limited mate availability and impaired cooperative behavior ( for instance, if too few individuals are available for cooperative foraging, hunting, and defens e ) Courchamp, Lud k, & Gascoign e ( 2008) devote a whole chapter in Allee Effects in Ecology and Conservation to elucid ating such mechanisms PAGE 68 62 The Allee effect is evidenced to occur i n numerous species ranging from the colonial Damaraland molerats to African Wild Dogs (Courchamp, Clutton Brock, & Grenfell, 1999) although it is not believed to affect the populations of most taxa (Sibly, Barker, Denham, Hone, & Pagel, 2005) A good review of the Allee effect is provided by Courchamp, et al. ( 1999) To begin incorporating the Allee effect into our m odel w e set out to f ind a critical threshold value (also called the Allee threshold ) above which the population will continue by ordinary logistic growth, and below which the population will decay Consider a population following normal logistic growth We begin by rever sing the sign of the righ t hand side of Eq. (2.9) to yield 1mdN N rN dt K ( 2 30) where 0mr As before, there are still two fixed points, only their stabilities h ave reverse d; now *0 N is stable and *NK is unstable. Finally, incorporating the Allee threshold parameter T yields 11,mdN NN rN dt TK ( 2 31) where 0mr and 0 TK (Gruntfest, et al. 1997 as cited by Courchamp, Lud k, & Gascoign e, 2008). As expressed previously, the stability of equilibrium points can be assessed qualitatively by analyzing the phase line portrait and additionally by linearizing the equation in the neighborhood of its equilibrium points. PAGE 69 63 Geometrical Analysis The phase line portrait ( Figure 2.9) depicts all three equilibrium solutio ns of Eq. (2.31) ( We will rely solely on this graphical method of determination while keeping in mind that equilibria can be found just as easily algebraically. ) Stable equilibria are found at *0 N and *NK and a single unstable equilibrium lies at *NT Solutions starting above the unstable equilibrium *NT converge to K as t and those below *NT converge to zero as t The concavity or conv exity of the solution curves ( Figure 2.10) is determined in the usual manner by finding the slope of the line tangent to the curve of the phase line (Figure 2.9). Points of inflection occur where the slope of the line tangent to the phase line equals zero Here, the first point of inflection is fo und between zero and T and the other is found between T and K R elative degrees of local stability may be distinguished by determining the steepness of the slopes around each equilibrium point, respectively We observe that the behavior of Eq. (2.31) reflects that of an Allee effect insofar as individuals below the Allee threshold T become extinct and those above the threshold progress towards their environmental carrying capacity K The specific type is called a strong Allee effect because populations below the threshold are driven to extinction, whereas, in cases of a weak Allee effect, populations below the threshold are merely hampered in their rates of growth. PAGE 70 64 Figure 2.9 Phase line portrait of the Allee effect, as described by Eq. (2.31) where T represents the critical Allee threshold value and K represents carrying capacity. Figure 2.10 Dynamics of the Allee effect, as described by Eq. (2.31) where T represents the critical Allee threshold value and K represents carrying capacity. PAGE 71 65 2.5 GROWTH MODEL WITH MULTIPLE EQUILIBRIA In the models outlined thus far populations always return to the same state of equilibrium following a perturbation ( unless they are pushed to extinction) Multistability is the possibility of alternative nontrivial equilibria such that the timing and magnitude of a disturbance may push the population into an alternate equilibrium state. Recalling the heuristic rolling ball analogy from Section 1.1 we can imagine a landscape upon which the ball is placed, where dips in the landscape represent basins of attraction, towards which the ball will roll, and mounds in the landscape represent unstable domains that repel the ball. In the current model, there are two basins of attraction towards which the ball might roll, depending on its positioning in the landscape; i.e., there are two stable nontrivial equilibria. If the population is at equilibrium, then a perturbation of sufficient magnitude is required in order for the population to converge towards the alternate equilibrium. Eq. (2.31) of the prior model possesses two nontrivial equilibria; however, the threshold equilibrium T is unstable, so there remains only one stable nonzero equilibrium solution. A ugment ing Eq. (2.31) such that a new term 1N L is included, and reversing the sign of the righthand side of the equation, we obtain 111,mdNNNN rN dtTKL ( 2 32) where 0 TKL and 0mr In ge neral, we may consider an equation of the form PAGE 72 66 11,n m i idN N rN dt K ( 2 33) having n many equilibrium points, where 12 nKKK Equilibri um points occur at iNK with alternating stability such that 2 iNK are stable and 21 iNK are unstable. Geometrical Analysis Figures 2.11 and 2.12 reveal the existence of two u nstable equilibria ( *0 N and *NK ) and two stable equilibria ( *NT and *NL ) of Eq. (2.32) The population may persist at either of the two alternative stable equilibria. If the population contains any number of individuals initially, then it is guaranteed to converge towards one of the equilibria, such that if the population s size N is below the ecological threshold K then it will converge towards the lower equilibrium T and if NK then it will converge towards the higher equilibr ium L PAGE 73 67 Figure 2.11 Phase line portrait of multistable growth model given by Eq. (2.32) with potential for coexistence at two levels, T and L where TL Figure 2.12 Dynamics of multistable growth model, given by Eq. (2.32) with potential for coexistence at two levels, T and L where TL PAGE 74 68 CHAPTER 3: MULTISPECIES POPULATION MODELS Until this point, our models have assumed that the population of a single species is constitutive of the entire system. We previously considered one dimensional models wherein single homogeneous populations fluctuate in the absence of i nterspecific relationships. Naturally ecosystems are constituted by populations belonging to multiple species and thus, the effects of these species on one anoth er should come under consideration. The narrowest case in which community dynamics can be modeled involves two interacting species These two species are modeled under the assumption that everything apart f rom that pair (that is to say, the environment, other species, etc.) is held constant. Therefore the two species are said to exist in isolation Maynard Smith ( 1974) presents an important inquiry Does the extent to which actual ecosystems show properties of persistence or stability depend on the fact that the pairwise interactions between species would likewise, in isolation, lead to stability and persistence? Henceforth w e turn our considerations to two dim ensional models, which will allow us to account for the effects of two species on one another. We will explore the ecological ramifications of competing populations, mutualistic populations, and predator prey interactions all of which exhibit characteristic nonlinear behaviors. Using qualitative approaches outlined prior, we will aim to elucidate these models and their ecological PAGE 75 69 significance. A dditionally, we will make use of the Jacobian (community) matrix of partial derivatives and its eigenvalues to as sess the systems and their stability. Ecological relationships are typically categorized by virtue of their interspecific interactions Table 3.1 describes these categories, most important of which are competition, mutualism, commensalism, and predation. The term symbiosis describes the interactions of species acting within the limits of any of these criteria; however, the term is restrictive insofar as species must live together in order to be called symbiotic Interaction type Effect on Species 1 Effect on Species 2 Competition ( ) Negative ( ) Negative Mutualism (+) Positive (+) Positive Predation (+) Positive ( ) Negative Commensalism (+) Positive (0) Neutral Amensalism ( ) Negative (0) Neutral Indifference (0) Neutral (0) Neutral Table 3.1. Basic categories of interspecific relationships T he (+) positive, or accelerating effect on a species S should be read as an increase in the birth rate of S or otherwise a decrease in the death rate of S Along the same lines, the ( ) negative, or inhibitory effect on a species S should be read as a decrease in the birth rate of S or an increase in the death rate of S The interaction types in Table 3.1 may be placed into three broad categories: cooperation, competition, and predation W e may broadly categorize competitive behaviors as those which occur when multiple species compete for the sam e commodity or resource. This may take the form of competitive exclusion ( ) or more rarely, amensalism ( 0). By contrast, cooperative behaviors are those having a positive net effect on the species PAGE 76 70 involved, namely mutualism (+ +) and commensalism (+ 0) Examples of two forms of mutualism are provided in Sections 3.2 and 3.3. Lastly, it should be noted that accor ding to these classifications, pred ation subsumes both predator prey interactions and hostparasitoid interactions. This assumption is favorabl e in terms of its simplicity, but it has the unfortunate conseq uence of emphasizing predation over analogous interactions such as parasitism or herbivore plant interactions (Maynard Smith, 1974) We will distinguish these terms on a caseby case basis. PAGE 77 71 3.1 INTERSPECIFIC C OMPETITION MODEL A major ecological concern involves competition between species sharing a habitat. How does a populations rate of change depend on its own population density and on the densities of competitor populations? We will begin with a classical model of competition based on the wor k of Lotka and Volterra wherein two species are assumed to have an inhibitory effect on each other. We start with the assumption that two species with respective populations 1N and 2N each grow logistically in the absence of the other as described by the following uncoupled logistic equations: 11 11 11, dNN rN dtK ( 3 1 ) 22 22 21. dN N rN dt K ( 3 2 ) The term s /iiNK for 1,2 i can be understood to represent intra specific competition, as we recall from the logistic model The carrying capacity parameter K therefore, is not explicitly determined by the environment, and thus its connection with the environment (including other species) is determined phenomenologically (Pastor, 2008) I f 1N and 2N are two species competing for a shared resource, however, then we may assume that the carrying capacity becomes a shared resource. As a result, each species inhibits the other: each individual of the first species causes a decrease in per capita growth of the second species, and vice versa The result is a symmetrical system of equations; that is, symmetrical with respect to the identity of each species (Pastor, 2008) Th is symmetry is PAGE 78 72 sensitive to parameters 12, KK as well as the p air of competition coefficients which describes the degree of competition each species has on the other W e arrive at the following coupled system of equations : 12 1 11 11, NN dN rN dt K ( 3 3 ) 21 2 22 21, NN dN rN dt K ( 3 4 ) where 1212,,, rrKK, and are positive. Obtaining Equilibrium P oint s Next, we may obtain the systems equilibrium points by finding values of 1N and 2N for which 120dNdN dtdt is satisfied. We may begin by finding values for which 10dN dt : ** 1 12 11 1 *** 1 1120 1 0, 0 or dN NN rN dt K NNKN ( 3 5 ) F inding values for which 20dN dt is satisfied: ** 2 21 22 20 1 0, dN NN rN dt K ( 3 6 ) and plugging in the first value 10N in Eq. (3. 5) into Eq. (3.6) it is apparent that ** 2 220 or .NNK ( 3 7 ) Plugging the second value from Eq. (3.5) ** 112NKN into Eq. (3.6) we get *** 212 22 210 NKN rN K ( 3 8 ) PAGE 79 73 and values of 2N satisfying Eq. (3.8) are ** 21 22 0 and 1 KK NN We can identify the following equilibrium points: ** 12,0,0 NN ( 3 9 ) ** 121,,0 NNK ( 3 10) ** 122,0, NNK ( 3 11) ** 1221 12,, 11 KKKK NN ( 3 12) We should note that the equilibrium in Eq. (3.12) may not be positive for all possible values; however, a biologically relevant equilibrium must be positive. With this information at hand, we may continue our analysis by determining the stability of each equilibrium point, and viewing the qualitative beh aviors for each case. Geometrical Analysis The nullclines or zero growth isoclines of these equations allow us to better characterize the systems dynamics Nullclin es for 1N and 2N occur, respectively, whe re 10 dN dt and 20. dN dt ( 3 13) T he 1N nullclines 10dN dt are given by the equations 10, N ( 3 14) 112, NKN ( 3 15) PAGE 80 74 and the 2N nullclines 20dN dt are given by 20, N ( 3 16) 221. NKN ( 3 17) We may already see how nullclines could be come helpful in finding equilibrium solutions. In Figures 3.53.6, n ullcli nes are represented as dashed and solid red lines for 1N and 2N respectively The general behavior of the vector field depends on whether the nullclines intersect, their degrees of orientation, and their relative positioning In addition, we may assess the stability of each equilibrium point to determine the flow of trajectories in the phase plane. From the phase portraits illustrated in Figures 3. 2 3. 6 w e can observe the outcomes delineated in Table 3.2 and the parameter space diagram of Figure 3.1. PAGE 81 75 Figure 3.1 Parameter space of four competition scenarios, wh ere, 2 1K KA and 1 2K KB Value of Qualitative Observation Corresponding Figure 12 21, KK KK Competitive advantage : 1N prevails, 2N goes extinct. 3.2 3. 3 12 21, KK KK Competitive advantage : 2N prevails, 1 N goes extinct. 3. 4 12 21, KK KK Strong competition : Bistability : winners success depends on initial conditions. 3. 5 12 21, KK KK Weak competition : Coexistence: both populations remain in stable equilibrium. 3. 6 Table 3.2 Four possible scenarios of the Lotka Volterra competition model PAGE 82 76 Figure 3.2 Phase plane portrait of Lotka Volterra competition model described by Eqs (3.1) and (3.2) for 12 21,KK KK indicating a competitive advantage of 1N over 2N Figure 3 .3 Dynamics of Lotka Volterra competition model described by Eqs. (3.1) and (3.2) for 12 21,KK KK The corresponding solution curve is denoted by in Figure 3.2 Here, 110 NN and 220 NN as t PAGE 83 77 Figure 3.4 Phase plane portrait of L otka Volterra competition model, described by Eqs. (3.1) and (3.2) for 12 21,KK KK indicating a competitive advantage of 2N over 1N In the phase plane portraits of the two monoculture scenarios (Figures 3.2 and 3.4), there are no critical points in the first quadrant because the nullclines for 1N and 2N (not shown) do not intersect. The equilibria for both cases lie on the axes/boundaries, so they are called boundary equilibria (Pastor, 2008) The solid dots represent nodal sinks (stable), and the hollow dots re present nodal sources (unstable). Under the conditions 12 21,KK KK which we may rewrite as 12KK and 21KK we observe that the maximum carrying capacity of 1 N namely 1 K exceeds that of 2N when the maximum competitive effect of 2N is less than the maximum carrying capacity of 1N; thus, the1N nullcline is positioned above that of 2N in the phase plane, and the trajectories converge toward stable equilibrium at ** 121(,)(,0) NNK (Figure 3.2). PAGE 84 78 The reverse case is also true; under the conditions 12 21,KK KK the 2N nullcline is positioned above that of 1N, and trajectories converge towards a stable boundary equilibrium at ** 122(,)(0,) NNK (Figure 3.4). These results uphold Gauses principle of competit ive exclusion ; namely that when the competitive effect of species 1N does not overcome the carrying capacity 2K of species 2N then 1N will be driven toward extinction, resulting in a monoculture of 2N or vice versa This type of competition is referred to as interference competition (Hardin, 1960; Meszna et al. 2006) Figure 3.5 Phase plane portrait of Lotka Volterra competition model described by Eqs. (3.1) and (3.2) for 12 21,KK KK indica ting bistability: either 1N or 2N will prevail. A strong degree of competition is suggested by the case where 12 21,KK KK (Figure 3.5). The nullclines of 1N and 2N are represented by dashed pink and orange lines, PAGE 85 79 respectively. Their intersection marks the formation of a half stable saddle point (recall Figure 1.6) at 1221** 12 11(,),,KKKKNN which is graphically represented as a half filled dot. The resulting behaviors are described as bistable because the trajectories may converge to two possible equilibrium points given the same parameter values. Viewing the phase portrait, we observe that either species may prevail as the winning monocu lture, while the loser species is driven to extinction. The initial conditions determine which equilibrium point will be approached by population trajectories. The line dividing the two locally stable regions is called the separatrix (Pastor, 2008) The trajectories divided by this line converge locally to the nearest stable boundary equilibrium, namely either 1(,0) K or 2(0,) K Figure 3.6 Phase plane portrait of Lotka Volterra competition model described by Eqs. (3.1) and (3.2) for 12 21,KK KK indicating coexistence. PAGE 86 80 U nder conditions of weak competition, as depicted in Figure 3.6, both competition coefficients are low, such that 12 21,KK KK and as a result, both populations exist at a stable equilibrium Thus, for stable coexistence to occur, inter specific competition coefficients must remain below the intraspecific competition thresho lds (i.e., carrying capacities), which are imposed regardless of the presence or absence of the other species. O ne might interpret this result, prima facie in conflict with Gauses pr inciple of competitive exclusion. B ecause, however, interspecific competition levels are weak (in fact, weak er than those f or intraspecific competition) it can be concluded that the two species do not compete to a high enough degree for them to be considered true competitors. Thus, they are said to inhabit independent ecological niches and Gauses principle is contested (Hardin, 1960) Local Linearization While our geometrical analys i s appear s sensible, let us verify the results of our model by linearizing it about the systems equilibria. Near equilibrium points, the dynamics of the system may be approximated by du AuBv dt ( 3 18) dv CuDv dt ( 3 19) where 11uNN and 22vNN At equilibrium t he associated Jacobian matrix, or community matrix is given by PAGE 87 81 ** 12** 1112 11 12 1 1 ** 2221 22 12 222 2NNrKNN rN ff NN K K gg rKNN rN NN KK J ( 3 20) We can see that in the specific case of coexistence t he Jacobian is written 112112 11 coexistence 221221 22()() (1)(1) ()() (1)(1) rKKrKK KK rKKrKK KK J ( 3 21) The condition for stability of coexistence requires that *tr()0 J and *det()0 J We know that the trace is negative if 1,2 1,2,K K and that the determinant is positive if 1 Therefore the product of the intra specific density dependence coefficients is greater than those of inter specificity. The biological relevance of this reflects Gauses observations discussed prior, n amely that complete competitors cannot coexist (Britton, 2003) PAGE 88 82 3.2 FACULTATIVE MUTU ALISM MODEL Facultative mutualism is a condition in which both species benefit from their mutual association. It is distinguished from obligate mutualism to the extent that facultative species may survive in the absence of each other, and obligate species will perish in the absence of one another An interesting example of facultative mutualism is found in the case of the Boran people of Kenya, who, in search of honey, use the guidance of a bird called the greater honeyguide ( Indicator indicator ) to locate colonies of honeybees In doing so, a mutual benefit is granted to both species: the people receive food subsistence from the honey, and the bird receives the increased ability to feed on honeybee larvae and hive wax. Isack & Reyer s ( 1989) statistical analysis reveals significant correlations between each species interspecific communications as well as the increased mutual success in locating the honeybee hives I n reality mutualisms are not necessarily or often, symmetrical For instance, the fitness of species 1S may depend wholly on 2S while the fitness of 2S may depend only slightly on 1S This relation would be called an obligatefacultative mutualism. For the sake of brevity, we will o nly be considering symmetrical associations (i.e., facultativefacultative and obligateobligate) We will approach the problem in a similar fashion as the LotkaVolterra competition model by first assuming two species with populations 1N and 2N that grow logistically in each others absence. C hanging the sign of the interaction coefficients and from PAGE 89 83 negative to positive we obtain interaction terms that enhance growth rates rather than inhibit them. The adjusted equations governing mutual benefaction are denoted 1 12 11 1() 1, dN NN rN dt K ( 3 22) 2 21 22 2() 1, dN NN rN dt K ( 3 23) w h ere and clump many mechanisms together into phenomenological entities, which determine the strength of mutualistic benefaction from each species In the case of facultative mutualism, we set parameters 1r 2r 1K, 2K to positive qu antities, so that for a species in absence of its mutualist, the equilibrium population density will be equivalent to its carrying capacity iK, 1,2 i That is to say 1N will grow towards its carrying capacity 1K even in the absence of 2N and vi c e versa Obtaining Equilibrium P oint s We obtain the systems equilibrium points by finding values of 1 N and 2N for which 120dNdN dtdt is satisfied. The following equilibrium points are obtained: ** 120, 0NN ( 3 24) ** 112,0 NKN ( 3 25) ** 1220, NNK ( 3 26) ** 1221 12, 11 KKKK NN ( 3 27) We note that these equilibrium solutions bear a strong resemblance to those of the competition model, except for the change of signs of the critical point for coexistence. PAGE 90 84 Geometrical Analysis The 1N nullclines 1(0)dN dt are given by the equations 10, N ( 3 28) 112, NKN ( 3 29) and the 2N nullclines 2(0)dN dt are given by 20, N ( 3 30) 221. NKN ( 3 31) The dynamics of facultative mutualism differ between a weak case where 1 and strong case where 1 Nullclines intersect if 1, and they diverge if 1 Therefore, the only case in which a critical point occurs is the weak case. A stable node form s in the first quadr ant where bo th nullclines cross denoting stable coexistence ( Figure 3.7); however, if nullclines diverge, then they do not cross at any point and the populations undergo unbounded growth (Figure 3.8) in what Robert May has called an orgy of mutual benefaction ( 1981) The case of unbounded growth driven by positive feedback between b oth mutualists is not a realistic scenario. S ome researchers have modified the problem such that limits are imposed on the mechanism of mutual positive feedback. For instance, Wolin and Lawlor ( 1984) consider the impact of mutualism with respect to recipient density via six different models: one model with per capita benefits of mutualism independent of recipient density, three models with mutualism effects most pronounced at a high density of recipients, and two models with mutualism effects most pronounced at a low density of recipients The latter two low density models were unique in th e sense that they always produced a PAGE 91 85 stable equilibrium, even in cases of strong facultative mutualism ; that is to say, where interaction coefficients were 12 21,KK KK ( 1984) Figure 3.7 Phase plane portrait of Lotka Volterra weak facultative mutualism described by Eqs. (3.1) and (3.2) with 1. PAGE 92 86 Figure 3.8 Phase plane portrait of the Lotka Volterra strong facultative mutualism described by Eqs. (3.1) and (3.2) with 1. Local Linearization Given that coexistence occurs only in the weak case of facultative mutualism we may begin our analysis evaluating the stable critical point that occurs when 1 The Jacobian matrix mirrors that of the competition model, except for the reversal of signs within the numerators: PAGE 93 87 112112 11 coexistence 221221 22()() (1)(1) ()() (1)(1) rKKrKK KK rKKrKK KK J ( 3 32) We find the trace and determinant of *J by 22 12121221 12() tr() (1) rrKKrKrK KK J ( 3 33) 121221 12()() det() (1) rrKKKK KK J ( 3 34) I f *tr()0 J and *det()0 J then the system is stable. Examining the case for weak facultative mutualism where 1 we find that, indeed, *tr()0 J and *det()0 J indicating that the critical point is a stable node given the real, negative eigenvalues ( Figure 3.7). Therefore, coexistence is guaranteed for the case of weak facultative mutualism but not for strong facultative mutualism PAGE 94 88 3.3 OBLIGATE MUTUALISM MODEL In the case of obligate mutualism, we may use the same equations that were used for facultative mutualism : (3.22) and (3.23) except that the sign of parameters 1r 2r 1K, 2K is reversed from positive to negative. Changing the carrying capacities may appear counterintuitive; however, it simply requires that neither species can survive in the absence of the other. Therefore b oth species are said to be obligate mutualists. Ob ligate mutualism is exemplified by numerous species that rely on intracellular bacterial symbionts T hese endosymbionts in turn, rely on their hosts for survival and fecundity. T he results of Wernegreens ( 2002) genomic analysis of two obligate mutualists ( B. aphidicola of aphids and W. glossinidia of tsetse flies) reveal marked gene loss and an integration of metabolic function between endo symbiont and host. In a similar vein, the integration of functional biological machinery arises in e ndosymbiotic theory pioneered by Lynn Margulis, wherein the endosymbiotic union of bacteria is held to be responsible for the origins of mitochondria and chloroplasts in eukaryotic cells (Kozo Polyansky, 2010) Using the same equations as in the prior case of facultative mutualism to define the 1N and 2N nullclines, namely Eqs. (3.28) through (3.31) we should note the changes that occur due to 1K and 2K becoming negative quantities, namely that both species become reliant on each other for survival. Geometrical Analysis In the case where 1, the nullclines of 1N and 2N do not intersect in the first quadrant, and the stable node at (0,0) is the only equilibrium so both populations decay PAGE 95 89 to wards extinction (Figure 3.9) In this case, the two species interdependent relations are too weak for either species to benefit the other to the point of survival. The probability of a weak obligate relationship occurring in nature would be rare since the two species are wholly interdependent. Conversely, when 1 a saddle point fo rms at the point of intersection between the two nullclines (Figure 3.10) Here, if densities of mutualists are below the saddle point threshold then both populations decay towards extinction despite the strong nature of the ir interaction. If mutualist densities are sufficiently high once more, both populations engage in an orgy of mutual benefaction, where orbits diverge to infinity These results indicate that coexistence between mutualistic species in the LotkaVolterra models, whether facultativefacultative or obligateobligate, is possibly only if interspecific interactions are sufficiently weak. Pastor (2008) speculates that strong interspecies interactions ap pear to destabilize food webs. Additionally, the absence of complex eigenvalues of the Jacobian matrices prohibits the possibility of stable limit cycles from occurring. PAGE 96 90 Figure 3.9 Phase plane portrait of the Lotka Volterra weak obligate mutualism model described by Eqs. (3.22) and (3.23) with 1. PAGE 97 91 Figure 3.10 Phase plane portrait of the Lotka Volterra strong obligate mutualism model, described by Eqs. (3.22) and (3.23) with 1. PAGE 98 92 3.4 PREDATORPREY MODEL Recalling the phenomenological nature of the carrying capacity term K of the logistic model, we found that K therefore is independently derived, and has little to do with the actual surrounding environment. Here, we will begin by eliminating the phenomenological term K for both species such that a new carrying capacity is determined instead, from the interactions between both species (Pastor, 2008) Along these lines, it should be easy to determine whether the exponential growth of prey is stabilized by predation, and likewise whether the growth of predators is stabilized by the decline of their food source, prey. Consider two populations 1N and 2N which represent preys and a predator s, respectively. The coupled system of equations of the LotkaVolterra predator prey model follows as 1 1 12, dN rNN d N t ( 3 35) 2 122, dN NNN dt ( 3 36) where, i n the equation of the prey population (3.35) r is the Malthusian growth parameter and is an interaction coefficient determining the rate at which predation (i.e., prey death) may occur The second term of (3.35) assume s that the product of both species densities 12NN accounts for the fact that both species must meet in order for predation to occur, and the coefficient determines the probability of a predation event successfully occurring. In Eq. (3.36) of the predator population, is the interaction coefficient determining the amount of biomass transferred from prey to predator for each successful predation event. PAGE 99 93 The constant serves as the predator death rate parameter, which mediates the constant exponential decay of predators Th e Lotka Volterra predator prey model can be further elucidated by outlining its assumptions Prey growth proceeds exponentially without limit in the absence of predators, causing dynamics to proceed in a Malthusian fashion such that 11 dN dtr N (0) r Predator growth is dependent on prey abundance, and therefore, if prey are absent then predator growth decays exponentially such that 22 dN dtN (0) The predation rate is dependent o n the likelihood of a predator individual meeting a prey individual in a spatially homo geneous population distribution, providing the terms governing prey death and predator birth. That is to say, the predator growth rate is proportional to the number of p rey present ( 12NN where is a positive constant), and likewise, prey death is proportional to the number of predators present ( 12NN where is a positive constant). Geometrical Analysis Initially, we may view the general behavior of 1 N and 2N by plotting the 12(,) NN vector field (Figure 3.11) and by determining the systems nu llclines, which occur where 120dNdN dtdt Following along the lines of Pastor ( 2008) we factor out 1N and 2N on the right hand side of each equation, yielding 1 12(), dN NrN dt ( 3 37) 2 21(). dN NN dt ( 3 38) PAGE 100 94 The intersection of 1 N and 2N reveals th e systems point of equilibrium. T rivial nullclines are found at 10 N and 20 N with each axis representing the nullcline for the other species (Figure 3.11) Setting the terms in parentheses of Eqs. (3.37) and (3.38) equal to zero such that 20, rN ( 3 39) 10, N ( 3 40) we achieve nullclines for 1N and 2N For 1N, we get the nullcline 2 rN and for 2N we get 1N (Figure 3.11) Notice that the nullcline for each species is dependent on the values of the other species as o pposed to its own population density. Alternatively, the net dynamics of the model can be viewed in Figure 3.12, which confirms the steady oscillatory behavior, suggested initially by the vector field (Figure 3.11). When 2 rN the dec line of the prey population is less than its growth rate, and thus the prey pop ulation increases: 10dN dt L ikewise, when 2 rN the predation upon prey by predators dominates the prey growth factor, and therefore the prey population decreases : 10dN dt The same mode of analysis can be performed for viewing predator dynamics When 1N the conversion of prey biomass into predator biomass by means of predation is outweighed by the constant mortality of predators; therefore, the predator population declines: 20dN dt and when 1N the growth of predators via predation outweighs the mortality of predators, and therefore, the predator population grows: 20dN dt PAGE 101 95 Figure 3.11 Phase plane portrait of Lotk a Volterra predator prey system described by Eqs. (3.35) and (3.36) PAGE 102 96 Figure 3.12 Dynamics of Lotka Volterra predator prey system, given by Eqs. (3.35) and (3.36 ) Local Linearization Setting both equations (3.35) and (3.36) equal to zero, and solving for 1 N and 2N we yield the following equilibrium points: ** 120, 0 NN ( 3 41) and ** 12, r NN ( 3 42) Near equilibrium point s, the dynamics of the system may be approximated by du AuBv dt ( 3 43) dv CuDv dt ( 3 44) PAGE 103 97 where 11uNN and 22vNN Constants are given by the community matrix ** 12** 12 21 ** 21 12 (,) NNff NN rNN gg NN NN J ( 3 45) Evaluating the first equilibrium at the origin (0,0) the Jacobian becomes (0,0)0 0 r J ( 3 46) The eigenvalues determine the stability of the point, and are found by 2 1,2() 22 r r ( 3 47) 1, r 2. ( 3 48) In the LotkaVolterra equations both ,0 r ; therefore, 10 and 20 T herefore, the equilibrium point at (0,0) is a saddle point (recall Figures 1.5 and 1.6) Evaluating the second equilibrium at coexistence (,)r the Jacobian becomes coexistence0 0 r J ( 3 49) The e igenvalues therefore, of the equilibrium point at 1N and 2 rN are given by 1,2. r ( 3 50) These are purely complex eigenvalues which indicates that the equ ilibrium is a neutrally stable center (recall Figure 1.6). PAGE 104 98 CHAPTER 4: CONCLUDING REMARKS In Chapter 1 we provided an outline of general dynamical systems, modeling methodology, and the associated history of differential equations in modeling population dynamic s. We discovered the multitude of uses that differential equations have served in modeling ecological dynamics, not to mention their ease of analysis and parsimonious assumptions. Beginning with the simplest one dimensional model s in Chapter 2, we c arried out a gradual modification of each subsequent model and analyzed its behaviors and ecological consequences. In Chapter 3, we followed the same basic process for two dimensional models. The basis of this progression has served to demonstrate how one might go about formulating a complex model from simple foundations in a piecemeal fashion. The current state of affairs in ecological population modeling continues to find numerous uses for differential equations based modeling. The bulk of e cological rese arch however, has become drastically refocused as the need to address global issues, from climate change to naturally driven catastrophes, has become ever more prevalent. The classical approach of ecology, as Miao, et al. (2009) put it, h as shifted to a new era in which ecological science must play a greatly expanded role in improving the human condition by addressing the sustainability and resilienc e of socioecological sy stems Temporal and spatial scales are expanded in the se increasingly sophisticated global models, in turn availing long er term and long er range predictions (Miao, Carstenn, & Nungesser, 2009) PAGE 105 99 However, much like the unpredictability of weather patterns, most ecological systems are chaotic and in turn, difficult to forecast. Both the historical development of population dynamics (outlined in Section 1.4) and schematization of ecological modeling (Section 1.2) illuminate the progression towards globally informed ecological modeling practices. The mathematical models considered throughout this work have focused primarily on small sc ale systems containing relatively few processes, yet their principles are, to a large degree, generalizable to higher dimensions Stochasticity and spatiality were ignored throughout this work; however, stochastic (i.e., probabilistic ) techniques have beco me well established, and there are numerous methods that account for spatial dynamics, including agentbased models such as cellular automata, or partial differential equations. The most intense issues of global ecology include anthropogenic impacts on ear ths ecosystems and climate (Yue, Jorgensen, & Larocque, 2010) and global health issues. Humanitys response to mitigating such impacts relies to a large d egree, on modeling. The rapid changes sustained by human populations and their environment in recent times present an unprecedented demand for ecologists to forge new connections with planning, policy making, and risk and decisionanalysis For example, the MIDAS project funded by the National Institute of Health, develops models intended for policymakers, public health workers and other researchers interested in emerging infectious diseases. E nvironmental conservation policy and decisionmaking also relies heavily on well formulated models The CREAM project funded by the European Union serves as one such example, having initiated over a dozen ecological modeling projects involving chemical risk assessment (Schmolke, Thorbek, DeAngelis, & Grimm, 2010) Additionally, a number of PAGE 106 100 intergovernmental panels share data and enact decisionmaking based on applications of ecological modeling; these include Man and the Biosphere (MAB), World Climate Research Program (WCRP), the Intergovernmental Panel on Climate Change (IPCC), International Geosphere Biosphere Program (IGBP), International Human Dimensions Program on Global Environmental Change (IHDP), International Program of Biodiversity Scienc e (DIVERSI TAS) Millennium Ecosystem Assessment (MA), Earth System Science Partnership (ESSP) and Global Earth Observation System of Systems (GEOSS) among others (Yue, Jorgensen, & Larocque, 2010) C lassi cal infectious disease dynamics, including Kermack and McKendricks SIR model (Section 1.4), rely on differential equations and have been praised for their maximal parsimony (Epstein, 2009) Such models, however, turn out to be poorly suited to capturing the complex behaviors of global spread, and many investigators have opted instead for individual or a gentbased models (ABMs ). These models are formulated from the bottom up by individuals whose interactions and behaviors are governed by a set of rules. A s Epstein writes, a gents can be made to behave something like real people: prone to error, bias, fear and other foibles (2009) Thus, ABMs embrace the complexity of social networks and the interactions between individuals (2009) albeit at a heightened computational expense Their utility is exemplified by virtue of the fact that they provide a base case ( as opposed to crystal ball forecast ) of material scenarios W ith the base case in hand, the model can be run repeatedly under varying circumstances to examine the effects of different environmental conditions and agent based configurations (2009) thereby informing p olicy and decisionmaking. PAGE 107 101 REFERENCES Bacar, N. (2011). A Short History of Mathematical Population Dynamics (p. 160). London: Springer. Berryman, A. A., & Kindlmann, P. (2008). Population systems: a general introduction Population (English Edition) (p. 222). London: SpringerVerlag. Britton, N. F. (2003). Essential Mathematical Biology Mathematical Medicine and Biology (Vol. 20, p. 370). London: Springer. Courchamp, F., Clutton Brock, T., & Grenfell, B. (1999). Inverse density dependence and the Allee effect. Trends in ecology & evolution (Personal edition) 14(10), 405410. Allee Effects in Ecology and Conservation (p. 256). Oxford: Oxford University Press. EdelsteinKeshet, L. (2005). Mathematical Models in Biology (p. 586). Philadelphia, PA: SIAM. Epstein, J. M. (2009). Modelling to contain pandemics. Nature 460(7256), 687. Gause, G. F. (1932). Experimental studies on the struggle for existence. Journal of Experimental Biology 9 389 402. Gause, G. F. (1934). The Struggle for Existence. Baltimore, MD: Williams & Wilkins. Gilpin, M. E., & Ayala, F. J. (1973). Global models of growth and competition. Proceedings of the National Academy of Sciences of the United States of America 70(12), 35903. Hardin, G. (1960). The Competitive Exclusion Principle. Science 131(3409), 12921297. Howard, L. O., & Fiske, W. F. (1911). Importation into the U.S. of the Parasites of the Gipsy Moth and the Browntail Moth (p. 311). Isack, H. A., & Reyer, H. U. (1989). Honeyguides and honey gatherers: interspecific communication in a symbiotic relationship. Science 243 (4896), 1343 1343. American PAGE 108 102 Association for the Advancement of Science, 1333 H St, NW, 8 th Floor, Washington, DC, 20005, USA. Jennings, H. S. (1942). Biographical Memoir of Raymond Pearl (1879 1940). National Academy of Sciences of the United States of America Biographical Memoirs (Vol. 31, pp. 293348). Kaplan, D., & Glass, L. (1995). Understanding Nonlinear Dynamics (p. 440). New York, NY: Springer. Kermack, W. O., & McKendrick, A. G. (1927). A Contribution to the Mathematical Theory of Epidemics. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 115(772), 700721. The Royal Soc iety. Kingsland, S. E. (1985). Modeling nature: Episodes in the history of population ecology (p. 267). Chicago: University of Chicago Press. Kot, M. (2001). Elements of Mathematical Ecology (p. 447). Cambridge: Cambridge University Press. Kozo Polyansky, B. M. (2010). Symbiogenesis: a new principle of evolution (V. Fet & L. Margulis, Eds.) Theory in biosciences = Theorie in den Biowissenschaften (Vol. 128, p. 240). Cambridge, MA: Harvard University Press. Lorenz, E. N. (1963). Deterministic Nonperiodic Flo w1. Journal of the Atmospheric Sciences 20(2), 130141. Lotka, A. J. (1907). Relation between birth rates and death rates. Science 26(653), 2122. Lotka, A. J. (1920). Analytical note on certain rhythmic relations in organic systems. Proceedings of the N ational Academy of Sciences of the United States of America 6 (7), 410415. National Academy of Sciences. MacArthur, R., & Wilson, E. O. (1967). The Theory of Island Biogeography (p. 203). Princeton, NJ: Princeton University Press. MacArthur, Robert, & Lev ins, R. (1967). The Limiting Similarity, Convergence, and Divergence of Coexisting Species. The American Naturalist 101(921), 377385. The University of Chicago Press for The American Society of Naturalists. Malthus, T. (1798). An Essay on the Principle of Population, as it Affects the Future Improvement of Society with Remarks on the Speculations of Mr. Godwin, M. Condorcet, and Other Writers (Vol. 6, p. 340). London. May, R. M. (1981). Theoretical Ecology: principles and applications (2nd ed., p. 489). O xford: Wiley Blackwell. PAGE 109 103 May, R. M. (2001). Stability and complexity in model ecosystems (p. 235). Princeton University Press. Maynard Smith, J. (1974). Models in Ecology (p. 146). Cambridge: Cambridge University Press. Meszna, G., Gyllenberg, M., Psztor, L., & Metz, J. a J. (2006). Competitive exclusion and limiting similarity: a unified theory. Theoretical population biology 69 (1), 6887. Miao, S., Carstenn, S., & Nungesser, M. (Eds.). (2009). Real World Ecology: LargeScale and Long Term Case Studies and Methods (p. 303). New York, NY: Springer. Murray, J. (2002). Mathematical Biology: An Introduction, Volume 1 (3rd ed., p. 551). New York, NY: Springer. Pastor, J. (2008). Mathematical Ecology of Populations and Ecosystems Ecosystems (p. 319). Oxford: Wiley Blackwell. Pearl, R. (1999). The Raymond Pearl Collection. The Alan Mason Chesney Medical Archives of The Johns Hopkins Medical Institutions Retrieved January 5, 2011, from http://www.medicalarchives.jhmi.edu/sgml/pearl.html. Pear l, R., Miner, J. R., & Parker, S. L. (1927). Experimental Studies on the Duration of Life. XI. Density of Population and Life Duration in Drosophila. The American Naturalist 61(675), 289318. The University of Chicago Press for The American Society of Nat uralists. Pearl, R., & Reed, L. J. (1920). On the Rate of Growth of the Population of the United States Since 1790 and its Mathematical Representation. Proceedings of the National Academy of Sciences 6 (6), 275288. Popper, K. R. (1992). The logic of scien tific discovery New Yorker, The (p. 513). New York, NY: Routledge. Schmolke, A., Thorbek, P., DeAngelis, D. L., & Grimm, V. (2010). Ecological models supporting environmental decision making: a strategy for the future. Trends in Ecology & Evolution 25 4 79486. Retrieved from http://linkinghub.elsevier.com/retrieve/pii/S016953471000100X. Schoener, T. W. (1986). Mechanistic Approaches to Community Ecology: A New Reductionism? American Zoologist 26(1), 81106. Oxford University Press. Sibly, R. M., Barker, D., Denham, M. C., Hone, J., & Pagel, M. (2005). On the regulation of populations of mammals, birds, fish, and insects. Science 309 (5734), 60710. Sigler, L. E. (2002). Fibonacci s Liber Abaci. London: Springer Verlag. PAGE 110 104 Strogatz, S. H. (1994). Nonlinear D ynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering (p. 498). Reading, MA: Perseus Books. Verhulst, P. F. (1845). Recherches mathmatiques sur la loi d accroissement de la population. Nouv. mm. de l Academie Royale des Sci et Belles Lettres de Bruxelles 18 1 41. Wernegreen, J. J. (2002). Genome evolution in bacterial endosymbionts of insects. Nature reviews. Genetics 3 (11), 85061. Wolin, C. L., & Lawlor, L. R. (1984). Models of facultative mutualism: density effects. A merican Naturalist 124 (6), 843 862. JSTOR. Yue, T. X., Jorgensen, S. E., & Larocque, G. R. (2010). Progress in global ecological modelling. Ecological Modelling Elsevier B.V. 