Mathematical Biology
The course is planned for the M.Sc. in Engineering Mathematics at ICT in their third semester. The course is also open for PhD candidates enrolled at ICT. The following topics will be covered. I mainly use the black board as the best medium for communication with the students. However, I shall add many materials time to time. The materials will be uploaded here. Each module will be supported by R Programming also. Total hours = 60. The course feedback will be appreciated that would help me to enhance the effective learning and teaching in this course.
Acknowledgment: I sincerely acknowledge the contribution of Kamlesh Suthar for his help in writing the R codes. Without his contribution, it would be extremely difficult for me to update the page timely throughout the semester continuously. He sat for many of my classes and immediately put those concepts into R codes.
The assignments will be posted online (in this page only).
Evaluation Plan:
Continuous assessment - I (Module - I) : 10
Continuous Assessment - II(Module - II & III) : 10
Continuous Assessment - III(module-IV & V) : 10
(Best two out of the three will be considered)
Midsem Examination : 30
Endsem Examination : 50
Click the link for tutorials. The tutorials mainly contain selected problems from the books mentioned in the reference list.
Tutorial - I Tutorial - II Tutorial - III Tutorial - IV Tutorial - V Tutorial - VI Tutorial - VII
Continuous Assessment - I, Continuous Assessment - II, Mid Semester, Continuous Assessment - III, End Semester
The solution of some selected problems in Midsem
Notes on Linear Birth-Death Process (Also contains solutions of Tutorial -VII)
The tutorial sessions will be arranged with the students' consultation.
Introduction and basic concepts
Concept of population, species, metapopulation, community
Relative growth rate (RGR) or per capita growth rate (PGR), Absolute growth rate (AGR)
Module - I: Single Population dynamics
The concept of overlapping generations (Differential equation models)
Population growth models: Mathematical formulation with biological interpretations
Unbounded growth aka exponential model,
The logistic growth model,
The theta-logistic growth model,
Allee growth model (Model with mate limitation)
Harvesting models
Growth models with varying coefficients: Gompertz model, sub-exponential model, Power law model, the concept of variable carrying capacity
Computation of growth characteristics: Concept of density-dependent and time-dependent growth model and computation of points of inflections for both the cases.
Checking for the stability of equilibrium points for single population dynamics
R companion: Solving differential equations in R (one dimensional); Basic plots of the above growth models (RGR or AGR curves). Almost all of these materials are covered in Bioinformatics and Statistics Lab-I.
Module - II: Single Population Dynamics
The concept of nonoverlapping generations ( Difference equation models)
Population growth models: Mathematical formulation with biological interpretations
Geometric growth
Discrete logistic growth model
Ricker growth model
Discrete Allee growth model (model with mate limitation and presence of generalist predator)
Stability of fixed points for one-dimensional map, Lyapunov exponents, Period doubling, Chaos
R companion: Bifurcation diagrams, computation of Lyapunov exponent, Time series plots, Understanding basic characteristics of the above models through plots.
Bonus: Concept of immigration and its importance in controlling chaotic dynamics
Delay models: lagged density dependence. Stability analysis
Module-III: Interacting populations
Classical Lotka-Volterra Predator-Prey system (historical importance)
Modeling predator-prey dynamics under different biological assumptions
Functional Response
Stability of equilibrium points for predator-prey system
Routh-Hurwitz criteria
Existence of limit cycle
Bendixson's negative criterion (with proof)
Bendixson-Dulac negative criterion
Checking for the stability of limit cycle
Discrete time predator-prey models
R Companion: Solving differential equations describing the predator-prey system in R. Time series plots, phase portrait. Visualizing stable-unstable equilibria, limit cycles graphically
Module-IV: Stochastic Population Dynamics
Demographic and environmental stochasticity
Birth and death process
Stationary and quasi-stationary distribution
Computation of extinction time in a stochastic environment
R Companion: Simulation of stochastic population processes
Module-V: Some additional topics (Only one of the following will be covered/ May be used for self-study or group-study and presentation by the students)
Concept of age-structured models
Concept of spatially structured models
Introduction to Disease models (Basic reproduction number, SI, SIR, SIRS, Next generation matrix)
Introduction to Food chain models
References:
Mark Kot, 2001. Elements of Mathematical Ecology, Cambridge University Press, Cambridge.**
Murray, J. D. 1989. Mathematical Biology, Springer-Verlag, Berlin **
Horst R. Thieme, 2003, Mathematics in Population Biology, Princeton University Press, Princeton.
Eric Renshaw. 1991. Modelling Biological Populations in Space and Time. Cambridge University Press
Stevens, M. Henry, A Primer in Ecology with R, Springer, 2009 **
Brauer, F. and Castilllo-Chavez, Mathematical Models in Population Biology and Epidemiology (2nd Edition, Springer, 2012)**
Edelstein-Keshet, Mathematical Models in Biology, Classics in Applied Mathematics
Module-II: R Companion
Population growth models governed by discrete difference equations:
Discrete time logistic difference equation. The following function is used to simulate the population time series data with fixed parameter values. The parameters are r and K, which represents the intrinsic growth rate and carrying capacity respectively. The discrete logistic equation is given by N(t+1) = N(t) + rN(t)[1-N(t)/K], N(0) = N0. No is the initial population size. We shall investigate the nature of the solutions for different values of the parameter r.
discrete.logistic.sim = function(K = 100, rd = 2.8, N0 = 5, t = 30) {
N = c(N0, numeric(t))
for(i in 1:t) {
N[i + 1] = N[i] + rd * N[i] * (1 - N[i]/K)
}
return(N)
}
d = discrete.logistic.sim(K=100, rd = 2.8, N0=50, t=50)
plot(d, type = "b", xlab = expression(t), ylab = expression(N[t]))
Let us create the plots using different values of r and try to understand the nature of the solutions of the discrete logistic equation. The carrying capacity K is asymptotically stable only for 0<r<2. For 0<r<1, solutions tend monotonically towards the carrying capacity. For larger r, 1<r<2, small perturbations about the carrying capacity still return to the carrying capacity, but in an oscillatory manner. Finally, when r>2, small perturbations from the carrying capacity diverge, with oscillations.
par(mfrow= c(2,2))
d = discrete.logistic.sim(K=100, rd = 1, N0=50, t=50)
plot(d, type = "b", xlab = expression(t), ylab = expression(N[t]), ylim = c(0, 140), main = "stable fixed point")
d = discrete.logistic.sim(K=100, rd = 1.9, N0=50, t=50)
plot(d, type = "b", xlab = expression(t), ylab = expression(N[t]), main = "asymptotic to a fixed point")
d = discrete.logistic.sim(K=100, rd = 2.2, N0=50, t=50)
plot(d, type = "b", xlab = expression(t), ylab = expression(N[t]), main = "stable 2-point limit cycle")
d = discrete.logistic.sim(K=100, rd = sqrt(6)+0.1, N0=50, t=50)
plot(d, type = "b", xlab = expression(t), ylab = expression(N[t]), main = "stable 4-point limit cycle")
Plot the return rate of the population. This is obtained by plotting N(t) along the x-axis and N(t+1) along the y-axis. In the following, I plot the time series and return map together.
par(mfrow=c(1,2))
Nts = discrete.logistic.sim()
t = 30; a = 0.01
plot(0:t, Nts, xlab="time (0:t)", ylab= "Population Size", type = "b")
abline(h = 1/a, lty = 3)
plot(Nts[1:t], Nts[1:t+1], type = "p", xlab=quote("N"["t"]), ylab=quote("N"["t+1"]))
Bifurcation diagram with respect to the parameter r. In the following, we implement the algorithm: We subdivide the interval of r in [0,3] into 401 points. For each value of r we generate the time series data governed by the logistic difference equation. We retain the last 50% observation (the choice is subjective but the diagram will not change for large t values) and plot them against the r values.
num.rd = 401; t = 200
rd.s = seq(1, 3, length.out = num.rd)
tmp = sapply( rd.s, function(r) discrete.logistic.sim( rd = r, N0 = 25, t = t ) )
tmp.s = stack(as.data.frame(tmp))
names(tmp.s) = c("N", "Old.Column.ID")
tmp.s$rd = rep(rd.s, each = t + 1)
tmp.s$time = rep(0:t, num.rd )
N.bif = subset(tmp.s, time > 0.5 * t)
plot(N ~ rd, data = N.bif, pch = ".", xlab = quote("r"["d"]),ylab="Population Size (N)")
Cobweb analysis for the logistic difference equation. Modify the following function for Cobweb analysis for discrete logistic model.
discrete.logistic_cobweb = function(K = 100, rd = 2.3, N0 = 5, t = 30) {
N = c(N0, numeric(t))
for(i in 1:t) {
N[i + 1] = N[i] + rd * N[i] * (1 - N[i]/K)
}
curve((1+rd)*x-rd*x^2/K , xlim = c(0,160), ylim = c(0,140) ,ylab = expression(N[t+1]),xlab = expression(N[t]), col = "red", lwd = 2, main = substitute(paste(r[d],"=",rd)))
curve(return(x), add = T,lty =4, lwd =2)
arrows(N[1], 0, N[1], N[2],length = 0)
for(i in 1:(length(d)-1)){
arrows(N[i], N[i+1], N[i+1], N[i+1], length = 0.00, col="blue" )
arrows(N[i+1], N[i+1], N[i+1], N[i+2], length = 0.00, col = "blue")
}
return(N)
}
d = discrete.logistic_cobweb(K=100, rd = 2.3, N0=40, t=50)
plot(d, type = "b", xlab = expression(t), ylab = expression(N[t]))
We shall understand the cobweb analysis for different values of r.
par(mfrow=c(2,2))
d = discrete.logistic_cobweb(K=100, rd = 1, N0=40, t=50)
d = discrete.logistic_cobweb(K=100, rd = 1.9, N0=40, t=50)
d = discrete.logistic_cobweb(K=100, rd = 2.2, N0=40, t=50)
d = discrete.logistic_cobweb(K=100, rd = sqrt(6)+0.1, N0=40, t=50)
Simulation of discrete time logistic difference equation with immigration
discrete.logistic.im = function(K = 100, rd = 2.8, N0 = 5, t = 30, c = 10) {
N = c(N0, numeric(t))
for(i in 1:t) {
N[i + 1] = N[i] + rd * N[i] * (1 - N[i]/K) + c
}
return(N)
}
Nts = discrete.logistic.im()
t = 30;K=100
plot(0:t, Nts, xlab= "time (0:t)", ylab="Population Size (N)", type = "b")
abline(h = K, lty = 3)
plot(Nts[1:t], Nts[1:t+1], type = "p", xlab=quote("N"["t"]),ylab=quote("N"["t+1"]), main="Return Map")
Bifurcation diagram with respect to the intrinsic growth rate r by varying the immigration parameter c. Check the figures carefully and note that increasing the value of the immigration parameter does not stabilize the population (recall the lecture). So, if we want to stabilize the population following a discrete logistic growth equation by means of constant immigration at each generation, will not be helpful.
par(mfrow = c(2,3))
num.rd = 201; t = 200
rd.s = seq(1.3, 3, length.out = num.rd)
c = seq(3, 8, 1)
for (i in 1:length(c)){
tmp = sapply( rd.s, function(r) discrete.logistic.im( rd = r, N0 = 25, t = t, c = c[i] ) )
tmp.s = stack(as.data.frame(tmp))
names(tmp.s) = c("N", "Old.Column.ID")
tmp.s$rd = rep(rd.s, each = t + 1)
tmp.s$time = rep(0:t, num.rd )
N.bif = subset(tmp.s, time > 0.5 * t)
plot(N ~ rd, data = N.bif, pch = ".", xlab = quote("r"["d"]), ylab="Population Size", main=paste("c=",as.character(c[i])))
}
Let us check the stabilization of the population following a Ricker growth equation with immigration. The following function is used to generate the population size following a discrete Ricker growth equation N(t+1) = N(t)exp[r(1-N(t)/K)], where r and K have the usual meaning. Note that the Ricker growth equation can be obtained from the continuous logistic growth equation by approximating 1/xdx/dt as d/dt(log x) and approximating it by taking dt = 1. The 1/dt (log x) is approximately equal to log(x(t+1)) - log(X(t)) and subsequently get the expression for X(t+1) as a function of X(t).
discrete.ricker.sim = function(K = 100, rd = 2.5, N0 = 80, t = 30){
N = c(N0, numeric(t))
for(i in 1:t){
N[i + 1] = N[i] * exp(rd * (1 - N[i]/K))
}
return(N)
}
par(mfrow=c(1,2))
Nts = discrete.ricker.sim()
t = 30; K = 100
plot(0:t, Nts, type = "l", xlab="Time(0:t)", ylab="Population Size")
abline(h = K, lty = 3)
plot(Nts[1:t], Nts[1:t+2], type = "p", xlab=quote("N"["t"]),ylab=quote("N"["t+1"]))
Bifurcation diagram with respect to the parameter r for the Ricker growth equation.
num.rd = 201; t = 200
rd.s = seq(1.8, 4, length.out = num.rd)
tmp = sapply( rd.s, function(r) discrete.ricker.sim( rd = r, N0 = 25, t = t ) )
tmp.s = stack(as.data.frame(tmp))
names(tmp.s) = c("N", "Old.Column.ID")
tmp.s$rd = rep(rd.s, each = t + 1)
tmp.s$time = rep(0:t, num.rd )
N.bif = subset(tmp.s, time > 0.5 * t)
plot(N ~ rd, data = N.bif, pch = ".", xlab = quote("r"["d"]))
Ricker growth equation with immigration
discrete.ricker.im = function(K = 100, rd = 2.8, N0 = 5, t = 30, c = 10) {
N = c(N0, numeric(t))
for(i in 1:t) {
N[i + 1] = N[i] * exp(rd * (1 - N[i]/K)) + c
}
return(N)
}
par(mfrow=c(1,2))
Nts = discrete.ricker.im()
t = 30; K = 100
plot(0:t, Nts, xlab="Time(0:t", ylab="Population Size")
abline(h = K, lty = 3)
plot(Nts[1:t], Nts[1:t+1], type = "p", xlab=quote("N"["t"]), ylab=quote("N"["t+1"]))
In the following, the bifurcation diagrams are drawn with respect to the growth parameter r at different values of the immigration rate c. Note that increasing c helps stabilizing the population.
par(mfrow = c(2,3))
num.rd = 201; t = 200
rd.s = seq(1.5, 4, length.out = num.rd)
c = c(2, 4, 5, 6, 7, 14)
for (i in 1:length(c))
{
tmp = sapply( rd.s, function(r) discrete.ricker.im( rd = r, N0 = 25, t = t, c = c[i] ) )
tmp.s = stack(as.data.frame(tmp))
names(tmp.s) = c("N", "Old.Column.ID")
tmp.s$rd = rep(rd.s, each = t + 1)
tmp.s$time = rep(0:t, num.rd )
N.bif = subset(tmp.s, time > 0.9 * t)
plot(N ~ rd, data = N.bif, pch = ".", xlab = quote("r"["d"]), main = paste("c = ", as.character(c[i])))
}
Observe that for the Ricker growth model the chaotic region is diminished by periodic halving. Why is it happening?
To answer this question, using the following code, we investigate the slope of the tangent to the return map at equilibria. Run the following piece of code and analyze the graphs generated.
r = 3.5; K = 80; c = c(1, 3, 5, 8, 10, 30)
par(mfrow=c(2,3)) for(i in 1:length(c)){
curve(x + r*x*(1-x/K)+c[i], 0 ,110, lwd = 2, ylim=c(0,150),
ylab = expression(N[t+1]), col = "red", xlab = expression(N[t])) abline(h=0, lwd=2)
curve(x*exp(r*(1-x/K)) +c[i], 0 ,120, lwd = 2, ylim=c(0,290),
curve(return(x), add = TRUE, col = "red", lwd = 2) } par(mfrow=c(2,3)) for(i in 1:length(c)){
}
ylab = expression(N[t+1]), col = "red", xlab = expression(N[t])) abline(h=0, lwd=2)
curve(return(x), add = TRUE, col = "red", lwd = 2)
Generating population size from Hassel model
dHassel = function(a = 0.5, ld = 2.2, b = 5.088, N0 = 2, t = 30){
N = c(N0, numeric(t))
for (i in 1:t) {
N[i + 1] = ld * N[i] / (1 + (a * N[i])^b )
}
return(N)
}
Nts = dHassel()
t = 30; a = 0.5
plot(0:t, Nts, type = "l", xlab="Time (0:t)", ylab="Population Size")
plot(1:t, Nts[1:t], type = "o")
par(mfrow=c(1,2))
Nts = dHassel()
t = 30; a = 0.5
#First order return map
plot(Nts[-length(Nts)], Nts[-1], type = "p", xlab=quote("N"["t"]),ylab=quote("N"["t+1"]))
#Second order return map
plot(Nts[-c(length(Nts),length(Nts)-1)], Nts[-c(1,2)], type = "p", xlab=quote("N"["t"]), ylab=quote("N"["t+2"]))
Bifuraction diagram with respect to lambda
num.l = 201; t = 200
l.s = seq(1, 4, length.out = num.l)
tmp = sapply(l.s, function(l) dHassel(ld = l, N0 = 10, t = t))
tmp.s = stack(as.data.frame(tmp))
names(tmp.s) = c("N", "Old.Column.ID")
tmp.s$l = rep(l.s, each = t + 1)
tmp.s$time = rep(0:t, num.l )
N.bif = subset(tmp.s, time > 0.5 * t)
plot(N ~ l, data = N.bif, pch = ".", xlab = expression(paste(lambda)))
Discrete Hassell model with immigration
dHasselIM = function(a = 0.6, ld = 1.2, b = 10, N0 = 5, t = 50, c = 0.5){
N = c(N0, numeric(t))
for (i in 1:t) {
N[i + 1] = ld * N[i] / (1 + (a * N[i])^b ) + c
}
return(N)
}
par(mfrow=c(2,2))
Nts = dHasselIM()
t = 50; a = 0.5
plot(0:t, Nts, type = "l", xlab="Time (0:t)", ylab="Population Size (N)")
Plot for one dimensional return map
plot(Nts[1:t], Nts[1:t+1], type = "p", xlab=quote("N"["t"]), ylab=quote("N"["t+1"])) #First order return map
plot(Nts[1:t], Nts[1:t+2], type = "p", xlab=quote("N"["t"]), ylab=quote("N"["t+2"])) #Second order return map
num.c = 201; t = 200
c.s = seq(0, 1, length.out = num.l)
tmp = sapply(c.s, function(c) dHasselIM(c = c, N0 = 10, t = t))
tmp.s = stack(as.data.frame(tmp))
names(tmp.s) = c("N", "Old.Column.ID")
tmp.s$c = rep(c.s, each = t + 1)
tmp.s$time = rep(0:t, num.l )
N.bif = subset(tmp.s, time > 0.5 * t)
plot(N ~ c, data = N.bif, pch = ".", xlab="Immigration c", ylab="Population Size (N)")
Hassel Model Special cases from the paper "Effect of Immigration on the Dynamics of Simple Ecological Models” lambda = 2.2, a = 0.5, b = 5.088, c= 0.5 lambda = 1.2, a = 0.6, b = 10, c= 1
Lyapunov exponent: Lyapunov exponent is a measure of sensitive dependence and an indicator of chaos. If it is positive, small differences are, on average, magnified and the resulting orbit is chaotic.
The following code is developed to compute the Lyapunov exponent for the logistic difference equation for different values of r. Run the code and understand the figures.
#Define function that returns N[t] for each t
F = function(K = 100, rd = 2.8, N0 = 5, t = 30) {
N = c(N0, numeric(t-1))
for(i in 1:t) {
N[i + 1] = N[i] + rd * N[i] * (1 - N[i]/K)
}
return(N)
}
#Define derivative of F which return df/dt at single point
dF = function(N, K = 100, rd = 2.8, t = 30) {
return( 1 + rd - 2 * rd * N/K)
}
#Ploting the Lyapunov Exponent
r_range = c(0,3)
plot(0, bty='n', pch='', xlim = r_range, ylab = expression(lambda), xlab = "r")
#plot(0,0, xlim = c(start_point, end_point), ylim = c(-5,1))
for(r in seq(r_range[1], r_range[2], by = 0.001)){
N = F(K=100, rd = r, N0=50, t=500)
lambda = 1/length(N)*(sum(log(abs(dF(N, rd=r, t=500, K=100)))))
points(r,lambda, pch = '.')
}
abline(h=0)
Module-III: R Companion
Solving the Predator-Prey Differential equations using R:
#########################################################################################################
# The following code is used to solve the predator-prey system,
# (a) The prey population N(t) grows exponentially with intrinsic growth rate r
# (b) The predator population P(t) dies exponentially in absence of preys at a rate m (mortality)
# (c) The loss of prey due to predator is proportional to both the number of preys and the number
# predators (mass action)
# (d) c is the per-prey capture rate by an individual predator
# (e) b is the nutritional gain of the predator.
#
# Under the above assumptions, the predator-prey system is given by
# dN/dt = rN - cNP; dP/dt = bNP - mP
# After nondimnesionalizing the system x=bN/m, y = cP/r,
# The following code solves the systems differential equation after
# after nondimensionalization
# dx/dt = r(1-y)x; dy/dt = m(x-1)y
#########################################################################################################
r = 0.8; m = 0.3
plot(1, 1, xaxt = "n", yaxt = "n", xlim = c(0,2), ylim = c(0,2), xlab ="x", ylab = "y")
# Prey zero isoclines
abline(h = 0, lwd = 2, col= "red")
abline(v = 1, lwd =2, col = "red")
# Predator zero isoclines
abline(v = 0, col = "blue", lwd= 2)
abline(h = 1, col = "blue", lwd = 2)
box()
# Equilibrium points
points(0,0, lwd= 2, col="black", cex=2)
points(1,1, lwd= 2, col="black", cex=2)
# Adding legends
legend("topright", legend = c("Prey isoclines", "Predator isoclines"), col = c("red", "blue"), bty = "n", lwd = c(2,2)); box()
library(deSolve)
Predator_Prey_Eqn = function(t, state, parameters){
with(as.list(c(state, parameters)),{
dx = r*(1-y)*x # equation for prey
dy = m*(x-1)*y # equation for predator
list(c(dx, dy)) # return the rate of change
})
}
state = c(x=0.7, y=0.8) # Initial condition (x_0, y_0)
parameters = c(r = r, m = m) # Value of the model parameters
times = seq(0,100, by = 0.01) # Time interval in which the solution to be computed
out = ode(y = state, times = times, func = Predator_Prey_Eqn, parms = parameters)
head(out) # Check first few entries of the solution vector
# Phase-portrait plot of the solution trajectories
lines(out[,"x"], out[,"y"], xlab = "Prey", ylab = "Predator", type = "l", col = "green", lwd = 2)
s = c(10, 500, 1000, 1200,1950, 2500, 5500)
arrows(out[s,"x"], out[s,"y"], out[s+1,"x"], out[s+1,"y"], col = 1, lwd =2, length = 0.15)
# Time series plot
plot(out[,"x"], xlab="time", ylab="density", type="l", lwd=2, col="red", ylim = c(0, max(out[,"y"]+0.5)))
lines(out[,"y"], xlab= "time", ylab = "density", type="l", lwd= 2, col ="blue")
legend("topright", legend = c("Prey", "Predator"), lwd = c(2,2), col = c("red", "blue"), bty = "n")
#########################################################################################################
# The following code is used to solve the predator-prey system
# (a) The prey population N(t) grows logistically with a fixed carrying capacity K and intrinsic growth rate r
# (b) The predator population P(t) dies exponentially in absence of preys at a rate m (natural mortality)
# (c) The loss of prey due to predator is proportional to both the number of preys and the number of predators (mass action)
# (d) c is the per-prey capture rate by an individual predator
# (e) b is the nutritional gain of the predator.
#
# Under the above assumptions, the predator-prey system is given by
# dN/dT = rN(1-N/K) - cNP; dP/dT = bNP - mP
# After nondimnesionalizing the system
# x=N/K, y = cP/r, t=rT, alpha = m/bK and beta = bK/r
# The following code solves the systems differential equation after
# after nondimensionalization
# dx/dt = x(1-x-y); dy/dt = beta(x-alpha)y
#########################################################################################################
beta = 1.2;
alpha = 0.2 # Check for alpha < 1 and alpha > 1
x_star = alpha; y_star = 1-alpha # Interior equilibria
discrim = alpha^2 - 4*alpha*beta*(1-alpha)# >0 implies stable node and <0 implies stable focus
# Plot of null clines
par(mfrow=c(1,1))
# Prey zero isoclines
curve(1-x, 0, alpha+1, type="l", xaxt = "n", yaxt = "n", bty="n", xlab = "x", ylab = "y", col = "red", lwd=2, xlim= c(0,alpha+1), ylim = c(0,2), cex.lab = 1.5)
abline(v = 0, lwd =2, col = "red")
# Predator zero isoclines
abline(v = alpha, col = "blue", lwd= 2)
abline(h = 0, col = "blue", lwd = 2)
# Equilibrium Points. Intersection of red and blue colour lines.
points(0,0, lwd= 2, col="black", cex=2)
points(1,0, lwd= 2, col="black", cex=2)
points(x_star, y_star, lwd = 2, col = "black", cex=2)
# Adding legends
legend("topright",legend=c("Prey isoclines","Predator isoclines"),col=c("red","blue"), bty="n", lwd=c(2,2))
box()
# Defining the systems of equations to be solved using deSolve package
install.packages("deSolve", dependencies = TRUE)
library(deSolve)
Predator_Prey_Eqn = function(t, state, parameters){
with(as.list(c(state, parameters)),{
dx = x*(1-x - y) # equation for prey
dy = beta*(x-alpha)*y # equation for predator
list(c(dx, dy)) # return the rate of change
})
}
state = c(x=1.0, y=1.5) # Intital condition (x_0, y_0)
parameters = c(beta = beta, alpha = alpha)
times = seq(0,100, by = 0.01)
# Solving the equations numerically
out = ode(y = state, times = times, func = Predator_Prey_Eqn, parms = parameters)
head(out) # Check first few entries of the solution vector
# Phase-portrait plot of the solution trajectories
lines(out[,"x"], out[,"y"], xlab = "Prey", ylab = "Predator", type = "l", col = "green", lwd = 2)
points(x_star, y_star, lwd = 2, pch=20, col = "black", cex = 1.5)
s = c(10, 500, 1000, 1200,1950, 2500, 5500)
arrows(out[s,"x"], out[s,"y"], out[s+1,"x"], out[s+1,"y"], col = 1, lwd =2, length = 0.15)
# Time series plot
plot(out[,"x"], xlab= "time", ylab="density", type="l", lwd= 2, col ="red", ylim = c(0, max(out[,"y"]+0.1)))
lines(out[,"y"], xlab= "time", ylab = "density", type="l", lwd= 2, col ="blue")
legend("topright", legend = c("Prey", "Predator"), lwd = c(2,2), col = c("red", "blue"), bty = "n")
# Create the Jacobian matrix symbolically for equilibrium analysis
F = expression(x - x^2 - x*y)
G = expression(beta*(x-alpha)*y)
# Computation of partial derivatives symbolically
F_x = D(F, "x"); F_x
F_y = D(F, "y"); F_y
G_x = D(G, "x"); G_x
G_y = D(G, "y"); G_y
# Create the Jacobian matrix that can be evaluated at equilibria
J = expression(matrix(c(eval(F_x), eval(F_y), eval(G_x), eval(G_y)), nrow = 2, byrow = TRUE))
print(J)
# Stability analysis of (0,0), check J evaluated at (0,0)
x = 0; y = 0
J_00 = eval(J); eigen(J_00); sum(diag(J_00)); det(J_00)
# Equilibrium analysis of (1,0)
x = 1; y = 0
J_10 = eval(J); eigen(J_10); sum(diag(J_10)); det(J_10)
# Equilibrium analysis of (alpha,1 - alpha)
x = alpha; y = 1-alpha
J_int = eval(J); eigen(J_int); sum(diag(J_int)); det(J_int)
# Plotting multiple trajectories in phase plane, with different initial conditions
curve(1-x, 0, alpha+1, type="l", xaxt = "n", yaxt = "n", bty="n", xlab = "x", ylab = "y", col = "red", lwd=2, xlim= c(0,alpha+1),
ylim = c(0,2), cex.lab = 1.5)
abline(v = 0, lwd =2, col = "red")
# Predator zero isoclines
abline(v = alpha, col = "blue", lwd= 2)
abline(h = 0, col = "blue", lwd = 2)
require(deSolve)
for(i in 1:7){
state = c(x=runif(1), y=runif(1)) # Intital condition (x_0, y_0)
parameters = c(beta = beta, alpha = alpha)
times = seq(0,100, by = 0.01)
# Solving the equations numerically
out = ode(y = state, times = times, func = Predator_Prey_Eqn, parms = parameters)
lines(out[,"x"], out[,"y"], xlab = "Prey", ylab = "Predator", type = "l", col = "green", lwd = 2)
points(x_star, y_star, lwd = 2, pch=20, col = "black", cex = 1.5)
s = c(10, 500, 1000, 1200,1950, 2500, 5500)
arrows(out[s,"x"], out[s,"y"], out[s+1,"x"], out[s+1,"y"], col = 1, lwd =2, length = 0.1)
}
box()
#########################################################################################################
# The following code is used to solve the predator-prey system,
# (a) The prey population N(t) grows logistically with a fixed carrying capacity K and intrinsic growth
rate r
# (b) The predator population P(t) dies exponentially in absence of preys at a rate m (natural mortality)
# (c) The loss of prey due to predator is modelled by the type-II functional response with a is the half saturation constant
# (d) c is the saturation point
# (e) b is the nutritional gain of the predator.
#
# Under the above assumptions, the predator-prey system is given by
# dN/dT = rN(1-N/K) - cNP/(a+N); dP/dT = bNP/(a+N) - mP
# After nondimnesionalizing the system
# x=N/a, y = cP/ar, t=rT, alpha = m/b and beta = b/r, gamma = K/a
# The following code solves the systems differential equation after after nondimensionalization
# dx/dt = f(x)(g(x)-y); dy/dt = beta(f(x)-alpha)y with f(x) = x/(1+x) and g(x) = (1-x/gamma)(1+x)
#########################################################################################################
gamma = 4.0; beta = 1.2; alpha = 0.8; x_star = alpha/(1-alpha)
f = function(x){
x/(1+x)
}
g = function(x){
(1-(x/gamma))*(1+x)
}
# Plot of predator and prey isoclines for different values of gamma
par(mfrow=c(2,2))
gamma_values = c(3, 5, 8, 11)
for(gamma in gamma_values){
curve(g(x), 0, gamma, type="l", bty="n", xlab = "x", ylab = "y",
col = "red", lwd=2, xlim= c(0,11), ylim = c(0,4), cex.lab = 1.5, lty=1)
abline(v=0, lwd=2, col="red")
abline(h= 0, lwd=2, col = "blue")
abline(v = alpha/(1-alpha), col = "blue", lwd= 2)
points(gamma,0, lwd=2, col="black", cex= 2)
points(0,0, lwd= 2, col="black", cex=2)
points(x_star, g(x_star), lwd = 2, col = "black", cex=2)
legend("topright", legend = c("Prey", "Predator"),
col = c("red", "blue"), bty = "n", lwd = c(2,2))
box()
}
require(deSolve)
Predator_Prey_Eqn = function(t, state, parameters){
with(as.list(c(state, parameters)),{
dx = x*(1-(x/gamma)) - (x*y)/(1+x) # Prey growth equation
dy = beta*(x/(1+x)-alpha)*y # Predator growth equation
list(c(dx, dy)) # return the rate of change
})
}
# Phase-plane analysis
gamma = 9.5 # Change the values of gamma and see the solution
par(mfrow=c(1,1))
curve(g(x), 0, gamma, type="l", xaxt='n', yaxt="n", bty="n", lty=1, xlab = "x", ylab = "y", col = "red", lwd=2,
xlim= c(0,max(alpha/(1-alpha), gamma)+1), ylim = c(0,5), cex.lab = 1.5)
points(x_star, g(x_star), lwd = 2, col = "black", cex=2)
abline(v=0, lwd=2, col="red")
abline(h= 0, lwd=2, col = "blue")
abline(v = alpha/(1-alpha), col = "blue", lwd= 2)
points(gamma,0, lwd=2, col="black", cex= 2)
points(0,0, lwd= 2, col="black", cex=2)
box()
times = seq(0,1000, by = 0.01)
state = c(x=8, y=4);
parameters = c(gamma = gamma, beta = beta, alpha = alpha)
out = ode(y = state, times = times, func = Predator_Prey_Eqn, parms = parameters)
lines(out[,"x"], out[,"y"], type = "l", col = "green", lwd = 2)
s = c(10, 500, 1000, 1200, 1700, 1950, 2500, 3500, 5500)
arrows(out[s,"x"], out[s,"y"], out[s+1,"x"], out[s+1,"y"], col = 1, lwd =2, length = 0.15)
if(gamma < x_star){
points(gamma,0, lwd=2, col="black", cex= 2, pch=20)
}
#Time series plot
plot(out[,"x"], xlab= "time", ylab = "density", type="l", lwd= 2, col ="red", ylim = c(0, 10))
lines(out[,"y"], xlab= "time", ylab = "density", type="l", lwd= 2, col ="blue")
legend("topright", legend = c("Prey", "Predator"), lwd = c(2,2), col = c("red", "blue"), bty = "n")
# Hopf bifurcation: As discussed in the class, for a Hopf bifurcation to occur, a stable or unstable focus must change its
stability and purely imaginary roots of the Jacobian matrix appears at equilibria. We study the behavior of the roots of
the Jacobian matrix for the above system
gamma = 5.0; beta = 1.2; alpha = 0.8; x_star = alpha/(1-alpha)
f = function(x){
x/(1+x)
}
g = function(x){
(1-(x/gamma))*(1+x)
}
F = expression(x*(1-(x/gamma)) - (x*y)/(1+x))
G = expression(beta*(x/(1+x)-alpha)*y)
F_x = D(F, "x"); F_x
F_y = D(F, "y"); F_y
G_x = D(G, "x"); G_x
G_y = D(G, "y"); G_y
# Create the Jacobian matrix that can be evaluated at equilibria
J = expression(matrix(c(eval(F_x), eval(F_y), eval(G_x), eval(G_y)), nrow = 2, byrow = TRUE))
print(J)
x=0; y=0
J_00 = eval(J); eigen(J_00); sum(diag(J_00)); det(J_00)
x =gamma; y =0
J_10 = eval(J); eigen(J_10); sum(diag(J_10)); det(J_10)
# Graphically demonstrating Hopf bifurcation; alpha is treated as the bifurcation parameters
par(mgp=c(3,1,0), mar=c(5,4,4,2)+0.1)
par(mfrow=c(2,2))
alpha_values = c(0.57, 0.65, 0.7)
for(alpha in rev(alpha_values)){
x_star=alpha/(1-alpha)
curve(g(x),- 0, gamma, type="l", bty="n", lty=1, xlab = "x", ylab = "y", col = "red", lwd=2,
xlim= c(0,max(alpha/(1-alpha), gamma)+1), ylim = c(0,6), cex.lab = 1.2)
points(x_star, g(x_star), lwd = 2, col = "black", cex=2)
abline(v=0, lwd=2, col="red")
abline(h= 0, lwd=2, col = "blue")
abline(v = alpha/(1-alpha), col = "blue", lwd= 2)
points(gamma,0, lwd=2, col="black", cex= 2)
points(0,0, lwd= 2, col="black", cex=2)
box()
times = seq(0,1000, by = 0.01)
state = c(x=8, y=4);
parameters = c(gamma = gamma, beta = beta, alpha = alpha)
out = ode(y = state, times = times, func = Predator_Prey_Eqn, parms = parameters)
lines(out[,"x"], out[,"y"], type = "l", col = "green", lwd = 2)
s = c(10, 500, 1000, 1200, 1700, 1950, 2500, 3500, 5500)
arrows(out[s,"x"], out[s,"y"], out[s+1,"x"], out[s+1,"y"], col = 1, lwd =2, length = 0.15)
}
alpha_values = seq(0.01, 0.9, by=0.001)
plot(1,1,xlim =c(0, gamma/(1+gamma)),ylim=c(-0.5,0.5), xlab = expression(alpha), ylab = expression(Re(lambda)), cex.lab = 1.2 )
for(alpha in alpha_values){
x = alpha/(1-alpha); y = g(x); x_star=alpha/(1-alpha)
J_int = eval(J)
lambda = eigen(J_int)$values
if(!is.complex(lambda)){
points(alpha, max(lambda), col = "blue", pch=20, lwd = 2)
}
if(is.complex(lambda)){
if(Re(eigen(J_int)$values[1]) < 0)
points(alpha,Re(eigen(J_int)$values[1]), col = "red", lwd=2, pch=20)
if(Re(eigen(J_int)$values[1]) > 0)
points(alpha,Re(eigen(J_int)$values[1]), col = "green", lwd=2, pch=20)
}
}
points(0.666, 0, cex=2, col = "black", pch = 20)
abline(h=0)
# Discrete time predator-prey system with linear functional response
Predator_Prey_Eqn = function(r = 0.5, b = 0.01, K = 100, c=0.01, m =1, t = 50, N0 = 80, P0 = 30){
NPs = data.frame(Ns = numeric(t+1), Ps = numeric(t+1))
NPs[1, ] = c(N0, P0)
for (t in 1:t) NPs[t + 1, ] = {
N = NPs[t, 1] + r * NPs[t, 1] * (1 - NPs[t, 1]/K) - c * NPs[t, 1] * NPs[t, 2]
P = NPs[t, 2] + b * NPs[t, 2] * NPs[t, 1] - m*NPs[t,2]
c(N, P)
}
return(NPs)
}
out = Predator_Prey_Eqn()
par(mfrow= c(1,2))
plot(out[,1], xlab = "time", ylab = "prey", pch = "*")
plot(out[,2], xlab = "time", ylab = "predator", pch = "*")
# Host Parasitoid system
HostParasitoid = function(lambda = 2, a = 1, c = 1, t = 10, N0 = 20, P0 = 10){
NP = data.frame(N = numeric(t+1), P = numeric(t+1))
NP[1, ] = c(N0, P0)
for (t in 1:t) NP[t + 1, ] = {
N = lambda* NP[t, 1] * exp(-a* NP[t,2])
P = c* NP[t, 1]*(1-exp(-a*NP[t,2]))
c(N, P)
}
return(NP)
}
out = HostParasitoid(lambda = 1.5, a=0.8, t=50)
par(mfrow=c(1,2))
plot(out[,1], xlab = "time", ylab = expression(N[t]), pch= "*", col = "red")
plot(out[,2], xlab = "time", ylab = expression(P[t]), pch = "*", col = "blue")