Moose and wolf population dynamics on Isle Royale

If the virus affects only the wolves (canine parvovirus), why might the moose have crashed so badly?

Phase plane diagram of predator prey dynamics

## # A tibble: 53 × 3
##     Year Moose  Wolf
##    <int> <int> <int>
##  1  1959   560    20
##  2  1960   600    22
##  3  1961   625    22
##  4  1962   640    23
##  5  1963   660    20
##  6  1964   700    26
##  7  1965   720    28
##  8  1966   760    26
##  9  1967   900    22
## 10  1968  1050    22
## # … with 43 more rows

trying to estimate model parameters

Raw data on N, P, and annual per capita growth rates.

n <- nrow(dfw)
Pt <- dfw$Wolf[-n]
Nt <- dfw$Moose[-n]
Pt1 <- dfw$Wolf[-1]
Nt1 <- dfw$Moose[-1]

pcgrN <- log( Nt1/Nt )

qplot(Pt, pcgrN, geom=c("point")) + geom_smooth(method="lm", se=FALSE)
## `geom_smooth()` using formula 'y ~ x'

Maximum prey per capita growth rate, \(r\), should be the y-intercept.

Using slopes and intercepts to estimate parameters

\[\frac{1}{N}\frac{dN}{dt} = r-aP\]

r as the y-intercept, a as the slope.

cfsN <- coef( mN <- lm(pcgrN ~ Pt) )
(r <- as.numeric(cfsN[1]));
## [1] 0.08554478
( a <- abs( as.numeric(cfsN[2] ) ) )
## [1] 0.003709956

\[\frac{1}{P}\frac{dP}{dt} = baN - m\]

m as the intercept, and b as the slope of annual per capita growth rate vs. aN.

pcgrP <- log( Pt1/Pt)
aNt <- a*Nt
qplot(aNt, pcgrP, geom=c("point")) + geom_smooth(method="lm", se=FALSE)
## `geom_smooth()` using formula 'y ~ x'

cfsP <- coef( mP <- lm(pcgrP ~ aNt)  )
(m <- abs(as.numeric( cfsP[1]) ) ) ; 
## [1] 0.1684191
( b <- as.numeric( cfsP[2]) )  
## [1] 0.04345842

Solve the ODE

predpreyLV <- function(t, y, p) {
  N <- y[1]; P <- y[2]
    with( as.list( c(y,p) ), {
      dN <- r*N - a*N*P
      dP <- b*a*N*P - m*P
      return( list( c(dN, dP) ))
    })
  }
predpreyPRD <- function(t, y, p) {
    N <- y[1]; P <- y[2]
    with( as.list( c(y,p) ), {
      dN <- r*N - a*N*P/P^z
      dP <- b*a*N*P/P^z - m*P
      return( list( c(dN, dP) ))
    })
  }
  parms <- c(r = r, a=a, b=b, m=m, z=0)
  y0 <- c(N=m/(b*a) + 5, P=r/a + 100)
  y0 <- c(N=560, P=20)
    times <- seq(0, 52, by=1)
    out <- as.data.frame( ode(y=y0, t=times, fun=predpreyPRD, parms=parms) )
    op <- pivot_longer(out, cols=-time, names_to = "Populations", 
                       values_to="Number")
    ggplot(op, aes(x=time, y=Number, colour=Populations)) + geom_line() +
      labs(x="Years") + scale_color_discrete(breaks=c("N","P")) + 
      facet_grid(Populations~., scales="free")

Using mixed methods to estimate parameters

  • r is the maximum annual per capita growth rate of the prey.
  • a is the solution to \(P^* = r/a\).
  • m is the absolute value of the minimum annual per capita growth rate of the predator.
  • b is the solution to \(N^* = m/(ba)\)

We estimate the equilibria with mean N and P.

r <- max(pcgrN)
a <- r/mean(Pt)
m <- abs( min(pcgrP) )
b <- m/(a*mean(Nt))

# What is the period for the resulting cycles
J <- matrix(c(0, b*r, -m/b, 0), ncol=2)
ev <- eigen(J)$values
2*pi/Im(ev[1])
## [1] 9.793018
parms3 <- c(r = r, a=a, b=b, m=m, z=0)
y0 <- c(N=560, P=20)
times <- seq(0, 52, by=1)
out3 <- as.data.frame( ode(y=y0, t=times, fun=predpreyPRD, parms=parms3) )
op <- pivot_longer(out3, cols=-time, names_to = "Populations", 
                       values_to="Number")
ggplot(op, aes(x=time, y=Number, colour=Populations)) + 
  geom_line() + labs(x="Years") + 
  scale_color_discrete(breaks=c("N","P")) + 
      facet_grid(Populations~., scales="free")

ggplot(out3, aes(N, P, colour=time)) + 
  geom_path(arrow=arrow(length=unit(3, "mm"), angle = 15, type = "open")) + 
  labs(subtitle="Each arrow is one annual time step.")

ggsave("MooseWolfPP_model3.png", height=5, width=5.5)

optimization

https://www.r-bloggers.com/learning-r-parameter-fitting-for-models-involving-differential-equations/

dfw$time <- dfw$Year-min(dfw$Year)
ssq <- function(parms){
 
 # inital concentration
 #parms <- c(r = r, a=a, b=b, m=m, z=0)
y0 <- c(N=560, P=20)
    t <- seq(0, 52, by=1)
 # time points for which conc is reported
 # include the points where data is available
 # parameters from the parameter estimation routine

 # solve ODE for a given set of parameters
 # out=ode(y=cinit,times=t,func=rxnrate,parms=list(k1=k1,k2=k2))
 out <- as.data.frame( ode(y=y0, t=times, fun=predpreyLV, parms=parms) )
 
 # Evaluate predicted vs experimental residual
 preddf=pivot_longer(out,-time,names_to="species",values_to="N")
 expdf=pivot_longer(dfw,-c(time, Year), names_to="species",values_to="N")
 ssqres=preddf$N-expdf$N
 str(ssqres)
 # return predicted vs experimental residual
 return(ssqres)
 
}
ssq2=function(parms){
 
 # inital concentration
 #parms <- c(r = r, a=a, b=b, m=m, z=0)
y0 <- c(N=560, P=20)
    t <- seq(0, 52, by=1)
 # time points for which conc is reported
 # include the points where data is available
 # parameters from the parameter estimation routine

 # solve ODE for a given set of parameters
 # out=ode(y=cinit,times=t,func=rxnrate,parms=list(k1=k1,k2=k2))
 out <- as.data.frame( ode(y=y0, t=times, fun=predpreyLV, parms=parms) )
 
 # Evaluate predicted vs experimental residual
 preddf=dt_pivot_longer(out,-time, names_to="species",values_to="N")
 preddf <- preddf[order(preddf$time, preddf$species),]
 expdf=dt_pivot_longer(dfw,-c(time, Year), names_to="species",values_to="N")
 expdf <- expdf[order(expdf$time, expdf$species),]
 ssqres=preddf$N-expdf$N
 #str(ssqres)
 # return predicted vs experimental residual
 return(ssqres)
 
}
parms <- parms3[1:4]
system.time(fitval<-nls.lm(par=parms,fn=ssq2))
##    user  system elapsed 
##   0.745   0.012   0.761
summary(fitval)
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## r 3.263e-01  4.106e-04   794.9   <2e-16 ***
## a 1.730e-02  2.023e-05   855.0   <2e-16 ***
## b 8.884e-02  7.205e-05  1233.0   <2e-16 ***
## m 1.379e+00  9.718e-04  1418.9   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 315.3 on 102 degrees of freedom
## Number of iterations to termination: 10 
## Reason for termination: Relative error between `par' and the solution is at most `ptol'.
parms <- coef(fitval)
y0 <- c(N=560, P=20)
times <- seq(0, 52, by=.01)
out4 <- as.data.frame( ode(y=y0, t=times, fun=predpreyLV, parms=parms3) )
op <- pivot_longer(out4, cols=-time, names_to = "Populations", 
                       values_to="Number")
df$time <- df$Year - min(df$Year)
df$Populations <- ifelse(df$Species == "Moose", "N", "P")
ggplot(op, aes(x=time, y=Number, colour=Populations)) + geom_line() + 
  geom_line(data=df, aes(time, N, colour=Populations)) +
  labs(x="Years") + scale_color_discrete(breaks=c("N","P")) + 
  facet_grid(Populations~., scales="free") + theme(legend.position="none")

ggsave("mooseWolffitted.png", width=4, height=3)

Discretized version for use in spreadsheets

dt <- 0.01
time <- seq(dt, 52, by=dt)
N <- numeric(52/dt)
P <- numeric(52/dt)
N[1] <- 560
P[1] <- 20

for(t in 1:(length(N)-1) ){
  N[t+1] <- N[t] + dt*(r*N[t] - a*N[t]*P[t])
  P[t+1] <- P[t] + dt*(b*a*N[t]*P[t] - m*P[t])
}
summary(P)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   10.30   12.81   19.45   23.29   33.50   45.00
qplot(time,N, geom="line")