• Lotka Volterra R Code
• Lotka Volterra In R

Curator:Frank Hoppensteadt

Aim: I am trying to numerically solve a Lotka-Volterra ODE in R, using de sde.sim function in the sde package. I would like to use the sde.sim function in order to eventually transform this system into an SDE. So initially, I started with an simple ODE system (Lotka Volterra model) without a noise term. The Lotka-Volterra ODE system. Simecol (simulation of ecological systems) is a lightweight R package that employs an object oriented paradigm for the implementation of dynamic simulation models. The package is intended to give users (students and scientists) an interactive environment to implement, distribute, simulate and document basic and advanced ecological models without the need to write long simulation. Regular Expression Syntax in R; Regular Expressions (regex) Reproducible R; Reshape using tidyr; Reshaping data between long and wide forms; RESTful R Services; RMarkdown and knitr presentation; RODBC; roxygen2; Run-length encoding; Scope of variables; Set operations; Shiny; Solving ODEs in R; The Lorenz model; Lotka-Volterra or: Prey vs.

Predator-prey models are arguably the building blocks of the bio- and ecosystems as biomasses are grown out of their resource masses. Species compete, evolve and disperse simply for the purpose of seeking resources to sustain their struggle for their very existence. Depending on their specific settings of applications, they can take the forms of resource-consumer, plant-herbivore, parasite-host, tumor cells (virus)-immune system, susceptible-infectious interactions, etc. They deal with the general loss-win interactions and hence may have applications outside of ecosystems. When seemingly competitive interactions are carefully examined, they are often in fact some forms of predator-prey interaction in disguise.

Lotka Volterra model competition model and predator prey model with equation - Duration: 30:13. Shomu's Biology 90,510 views.

A General Predator-Prey Model

Consider two populations whose sizes at a reference time (t) aredenoted by (x(t) ,) (y(t) ,) respectively. The functions (x) and (y) might denote population numbers or concentrations (number per area) orsome other scaled measure of the populations sizes, but are taken tobe continuous functions. Changes in population size with time aredescribed by the time derivatives (dot x equiv dx/dt) and (dot yequiv dy/dt ,) respectively, and a general model of interactingpopulations is written in terms of two autonomous differential equations[ dot x = x f(x,y) ][ dot y = y g(x,y)](i.e., the time (t) does not appear explicitly in thefunctions (x f(x,y)) and (y g(x,y))). The functions (f) and (g) denote the respectiveper capita growth rates of the two species. It is assumed that ( df(x, y)/dy<0 ) and ( dg(x, y)/dx>0. )This general model is often called Kolmogorov's predator-prey model (Freedman 1980, Brauer and Castillo-Chavez 2000).

Lotka-Volterra Model

In 1926, the famous Italian mathematician Vito Volterra proposed a differential equationmodel to explain the observed increase in predator fish (and correspondingdecrease in prey fish) in the Adriatic Sea during World War I. At the same time in the United States, the equations studied by Volterra werederived independently by Alfred Lotka (1925) to describe a hypothetical chemical reaction in which thechemical concentrations oscillate.The Lotka-Volterra model is the simplest model of predator-prey interactions. It is based on linear per capita growthrates, which are written as [f= b-p y] and (g=r x-d .)

• The parameter (b) is the growth rate of species (x) (the prey) in the absence of interaction with species (y) (the predators). Prey numbers are diminished by these interactions: The per capita growth rate decreases (here linearly) with increasing (y ,) possibly becoming negative.
• The parameter (p) measures the impact of predation on (dot x/x .)
• The parameter (d) is the death (or emigration) rate of species (y) in the absence of interaction with species (x .)
• The term (r x) denotes the net rate of growth (or immigration) of the predator population in response to the size of the prey population.

The Prey-Predator model with linear per capita growth rates is [dot x = (b - p y) x] (Prey)[dot y = (r x - d) y] (Predators)This system is referred to as the Lotka-Volterra model: it represents one of the earliest models in mathematical ecology.

Lotka Volterra R Code

The system can be integrated directly. In particular, any solution((x(t),y(t))) of the system satisfies the identity[ C = b ln y(t) - p y(t) - r x(t) + d ln x(t) ]for all (t ,) where the constant (C= b ln y(0) - p y(0) - r x(0) + dln x(0)) is determined by initial conditions and system parameters.(Here (ln x) denotes the natural logarithm of (x ,) etc.) Thequantity on the right hand side of the identity above is referredto as a conservation law, since it is constant along anysolution. Having a conserved quantity facilitates visualizingsolutions. In the figure we draw level-set contours of thesurface[z = b ln y - p y - r x + d ln x]in the first quadrant of the (xy)-plane. The contours describesolutions of the system determined by their initial data, and since they are closed curves, the solutions are periodic oscillations.

If (b>0 ,) there are two equilibria, (x=0, y=0) (extinction), and(x=d/r, y=b/p) (co-existence), and the surface [z = b ln y - p y - r x + d ln x]has a single peak at the latter equilibrium. Thecontour lines in the figure describe the classicprey-predator cycles observed in ecological systems.

The model above has been derived independently in the following fields:

• epidemics (Kermak and McKendrick 1927, 1932, 1933) ((b=0))
  • (x) are susceptible individuals and
  • (y) are infective individuals,
• ecology (Lotka 1925, Volterra 1926)
  • (x) are prey and
  • (y) are predators,
• combustion theory (Semenov 1935)
  • (x) and (y) are chemical radicals formed during H₂ O₂ combustion,
• economics (Galbraith 2006)
  • (x) is the populace and
  • (y) is a predatory institution,

and numerous other studies from diverse disciplines.

Alfred James Lotka (1880 - 1949, USA, chemist, demographer, ecologist and mathematician) was born in Lviv (Lemberg), at that time situated in Austria, now in Ukraine. He came to the United States in 1902 and wrote a number of theoretical articles on chemical oscillations during the early decades of the twentieth century, and authored a book on theoretical biology (1925). He then left (academic) science and spent the majority of his working life at an insurance company (Metropolitan Life). In that capacity he became president of the PAA (the Population Association of America).

Vito Volterra (May 3, 1860 - October 11, 1940, Italian, mathematician and physicist, best known for his contributions to mathematical biology) was born in Ancona, into a very poor family. Volterra showed early promise in mathematics before attending the University of Pisa. His work is summarised in his book Theory of Functionals and of Integral and Integro-Differential Equations (1930).After World War I, Volterra turned his attention to the application of his mathematical ideas to biology, principally reiterating and developing the work of Pierre François Verhulst. The most famous outcome of this period is the Lotka-Volterra model. In 1922, he joined the opposition to the Fascist regime of Benito Mussolini and in 1931 refused to take a mandatory oath of loyalty. He was compelled to resign his university post and membership of scientific academies, and, during the following years, he lived largely abroad, returning to Rome just before his death.

Kermack-McKendrick Model

There is herd immunity in predation and in epidemics. It isconvenient to frame this in terms of epidemiology where now we referto the prey as being susceptibles and the predators as beinginfectives. The infection dynamics is depicted by the graph (xto yto ,) indicating that susceptibles can become infectives and thatinfectives can be removed from the process (e.g., through death,quarantine, or inoculation). Consider a time interval that is shortcompared to reproduction of the susceptible population, i.e., let(b=0 .) If the initial susceptible population is so large that (x(0)> d/r ,)then we see from the second equation in the predator-prey model that initially(dot y>0 ,) which indicates that the infectives will initially morethan replace themselves by passing on the infection. However, ifthis condition is not satisfied, the infective population willdecrease. The critical value (Requiv rx(0)/d=1) is referred to asbeing an epidemic threshold or a tipping point for the process.

A short calculation shows that (x(t)) converges to a constant, say (x(t)to x^*,) where (x^*) can be found by solving the equation (C=r x^* - d ln x^* ,) as shown in the figure. Surprisingly, thisnumber is always greater than zero, which shows that somesusceptibles will always survive!

This phenomenon, which is referredto as herd immunity is observed in practice; in fact, thenumber (R) is published regularly for various diseases and localesas an epidemic control measure. It reflects the fact that thesusceptible population can be reduced to a level below whichinfectives will not increase. The model in this case is referred toas being the Kermack-McKendrick model of susceptible-infective interactions in epidemiology.

Jacob-Monod Model

Another approach to modeling the interaction between prey andpredators was developed to account as well for organisms (such asbacteria) taking up nutrients. There is a limited uptake rate thatsuch organisms are capable of, and the next model accounts forlimited uptake rates. Suppose now that (x) denotes the size of apopulation of feeders and that they are feeding on a chemicalspecies of concentration (y .) These are usually taken to representconcentrations of feeders and nutrients in solution rather than'head counts'. The Jacob-Monod model is[dot x = frac{V y}{K+y}x][dot y = -frac 1Y frac{V y}{K+y}x]where

• (V) is the uptake velocity,
• (K) is the saturation constant, and
• (Y) is the yield of (x) per unit (y) taken up.

Note that when (y=K ,) the uptake velocity is (V/2 ,) half the maximum; in practice, (y=K) is taken as a tipping point: If (y<K ,) then uptake is ignored. This is the canonical model of nutrient uptake (nutrient (y) is taken up by species (x)), and it underlies many calculations in biology,microbiology, and food engineering (digestion, beer, etc.). Thismodel was discovered independently in several diverse applications.It is akin to the Haldane-Briggs model and Michaelis-Menten model inbiochemistry, the Jacob-Monod model in microbial ecology, and theBeverton-Holt model in fisheries. It serves as one of the importantbuilding blocks in studies of complex biochemical reactions and inecology (Smith and Waltman 1997).

The quantity (C=x +Y y) is conserved (since the derivative withrespect to time of the right hand side is zero) where(C=x(0)+Yy(0) .) Substituting this into the first equation gives[dot x =frac{V (C-x)}{Y K+(C-x)}x .]The solutions can be found by quadratures, and these show that if(x(0)>0 ,) then (x(t)to C) as (ttoinfty ,) at which point the nutrient has been depleted.

Lotka Volterra In R

A typical use of this model is to describe a continuous-flow growthdevice, such as a chemostat, where there is continuous removal of nutrient and feedersand a continuous supply of fresh nutrient. The Jacob-Monod model is used to describe such a bacterial growth device, for example to determine conditions for a sustained dynamic equilibrium to exist by balancing growth due to uptake of nutrient with wash out of feeders.

Some predator-prey models use terms similar to those appearing in the Jacob-Monod model to describe the rate at which predators consume prey. More generally, any of the data in the Lotka-Volterra model can be taken to depend on prey density as appropriate for the system being studied. This is referred to as a functional response, an idea that is introduced and discussed by C. S. Holling (1959). Several different forms of functional response have been used in population models, but the Jacob-Monod form, also called the Holling type 2 form by ecologists, is one of the more common ones. Many other investigations of predator-prey models have involved functional responses. For example, (r) is replaced by (r ln(K/x)) (Gompertz, 1825), (r(K-x)/(K+epsilon x),) (Smith, 1963), and (rleft((K/x)^g-1right)) with (0< g le 1 ,) (Rosenzweig, 1971). These and other functional responses are also discussed in May (1974). Such mechanisms in the Lotka-Volterra model can stabilize or destabilize the system, for example resulting in predator extinction or in co-existence of prey and predators. This is in contrast to the plurality of cycles predicted by the original Lotka-Volterra model.

Logistic Equation

An interesting case for [dot x =frac{V (C-x)}{Y K+(C-x)}x]is when (V) and (YK)are very large compared to the other data in the model, but withtheir ratio being of moderate size, say (V/(YK)approx r .) Then wecan ignore the second term in the denominator and get[dot x = r(C-x)x .]This is referred to as the logistic equation. It, too, has arisen invarious disciplines, but one of its first appearance was due toVerhulst in the mid 19th century who used it to correct certaindeviations of Malthus's model ((dot x = r x)) from certain humanpopulation data. The number (C) is now referred to as being thecarrying capacity for the population; this corresponds tothere being no remaining nutrient in the Jacob-Monod model.

The logistic equation can be solved inclosed form by quadratures. This shows that (x(t)to C) as(ttoinfty ,) if (x(0)> 0 .)

Predation with Time Delays: Chaos in Ricker's Reproduction Equation

Time delays occur in biological systems, and they can producecomplicated dynamics. To model age structure (and other time delays)in a system, we take the approach that was introduced by Euler inthe 18th century. Let (x(t)) denote the number born at time(t .) Then (x(t-a)) denotes the number who were born (a) units ago.Suppose that there is a nonlinear effect such that the number ofnewborns (a) units ago who can participate in reproduction at time(t) is (h(x(t-a))) where the function (h) is called thereproduction function. This nonlinearity might be due to predation or environmental factors, as discussed earlier. Suppose those of age (a) havefertility (m(a) .) Then the equation below shows how manywill be born at time (t .) It is obtained by adding up across allages those who can participate weighted by their fertility:[ x(t) = x_0(t) + int_0^t m(a) h(x(t-a)), da ,]where (x_0(t)) denotes the lingering contributions to later birthsof the initial population. In the case where (h(x)= x) is a linearfunction, this equation is referred to as being therenewal equation, which is widely used in demography. Itssolution can be found using Laplace transforms (Keyfitz and Flieger 1971).However, if (h) is a nonlinear function, then that method ofsolution is not available.

The fertility function (m(a)) describes the fertility of those ofage (a .) If (m) is constant, then this equation is equivalent to adifferential equation. The other extreme occurs when fertility isfocused all at one age, as for salmon or cicadas. In this case,suppose that all fertility is due to those of age (a^* .) We write(m(a)=r,delta(a-a^*)) where (delta) is Dirac's delta function and(r) is the reproduction rate. The equation becomes[x(t) = x_0(t)+r h(x(t-a^*)) .]Let us focus only on the births at the times (a^*, 2a^*, 3a^*,ldots ;)we define (x_n=x(na^*)) with (x_1=x_0(a^*) .) Substituting this inthe equation above gives the recursion[x_{n+1}= r h(x_n) ,] for (n=1,2,3,4,ldots .)In the case of Malthus's model, (h(x) = x ,) and the solution is asimple geometric progression (x_n = r^n x_0 .) However, if (h) is asin Verhulst's case, where (h(x) = x ( C - x)_+ equiv max(C-x,0),) the recursionbecomes[x_{n+1}=r x_n(C-x_n)_+ ,]whose iterates can exhibit quite wild behavior. This was firstnoticed by the ecologist W.E. Ricker in the 1950's who used thefunction (h(x)=xexp{(-x)}) in studies of the dynamics of fisheries,although his work was largely ignored at the time. This reproductionfunction accounts for cannibalism (self predation) in that if thepopulation is small, the model looks like Malthus's, but for largepopulations, reproduction is strongly suppressed. The work was laterrediscovered by Robert May, who stimulated the nowprominent area of chaos (Gleick 1987).

Chaotic dynamics can be illustrated by a simple computer simulationusing Ricker's model. Since the maximum of (h) is 1, we need onlyconsider values of (xin [0, r] .) We divide this interval up into300 equal bins. The parameter (r) is fixed at one of 300 equallyspaced numbers in (rin [0,20] .) The recursion is[x_{n+1}=r x_ne^{-x_n} .]Beginning with (x_1=1 ,) we determine the sequence({x_{50},x_{51},x_{52},x_{53},dots,x_{300}}) and record in anarray, say (D(r,x) ,) the number of sequence elements that are ineach bin. We only start counting after 50 initialiterates to avoid transients. The result is a surface (D(r,x)) thatdescribes the dynamics. We plot this record looking from above, asshown in the figure (Hoppensteadt and Hyman 1977). The solution always settles into some structure, which might be highly complex.At first, there is a unique stable state. Near (r=8) this statebifurcates into an oscillation of period 2. Near (r=12) the period 2solution bifurcates into a period 4 solution, etc. By (r=18) thesolution set is quite complicated. This is in the range of chaoticdynamics. When the fertility function (m) has broader support, the solutionset is more complicated still.

References

• F. Brauer and C. Castillo-Chavez, Mathematical Models in Population Biology and Epidemiology, Springer-Verlag, Heidelberg, 2000.
• H. I. Freedman, Deterministic Mathematical Models in Population Ecology. New York: Marcel Dekker, 1980.
• J.K. Galbraith, The Predator State, Mother Jones, May/June 2006.
• J. Gleick, Chaos: The Making of a New Science, Viking Press, New York, 1987.
• C. S. Holling, The characteristics of simple type of predation and parasitism, Canadian Entomologist 91 (1959), 385-398.
• F.C. Hoppensteadt, J.M. Hyman, Periodic solutions of a logistics difference equation, SIAM J. Appl. Math. 58(1977)73-81.
• O. Kermack, A.G. McKendrick, Proc. Roy. Soc. A, 115(1927)700-721, 138(1932)55-83, 141(1933)94-122.
• N. Keyfitz, W. Flieger, Populations: Fact and Methods of Demography, W.H.Freeman, San Francisco, 1971.
• M. Kot, Elements of Mathematical Ecology, Cambridge University Press, 2001
• A.J. Lotka, Elements of physical biology. Williams and Wilkins, Baltimore, 1925.
• R. M. May, Stability and Complexity in Model Ecosystems, Princeton U. Press, NJ, 1974.
• N.N. Semenov, Chemical Kinetics and Chain Reactions, Clarendon, Oxford, 1935.
• H. Smith, P. Waltman, The Mathematical Theory of Chemostats, Cambridge U. Press, 1997.
• V. Volterra, Variazioni e fluttuazioni del numero d'individui in specie animali conviventi. Mem. R. Accad. Naz. dei Lincei 2(1926)31�113

Internal references

• Jeff Moehlis, Kresimir Josic, Eric T. Shea-Brown (2006) Periodic orbit. Scholarpedia, 1(7):1358.
• Philip Holmes and Eric T. Shea-Brown (2006) Stability. Scholarpedia, 1(10):1838.

External Links

See also

Dynamical Systems,Periodic Orbit,Population Dynamics,Phase Space,Stability