7 Multi-species model

7.1 Model runs

We build post-breeding/post-reproduction population models.

# TT: number of years the model will run
# nsim: number of trajectories 

wolfN0 <- matrix(nrow=TT, ncol=nsim)  # number of wolf puppies
wolfN1 <- matrix(nrow=TT, ncol=nsim)  # number of first year wolfs
wolfN2 <- matrix(nrow=TT, ncol=nsim)
wolfN3plus <- matrix(nrow=TT, ncol=nsim)
wolfN <- matrix(nrow=TT, ncol=nsim) # total wolf population
wolfP <- matrix(nrow=TT, ncol=nsim)

deerInN0 <- matrix(nrow=TT, ncol=nsim)
deerInN1 <- matrix(nrow=TT, ncol=nsim)
deerInN2 <- matrix(nrow=TT, ncol=nsim)
deerInN3plus <- matrix(nrow=TT, ncol=nsim)
deerOutN0 <- matrix(nrow=TT, ncol=nsim)
deerOutN1 <- matrix(nrow=TT, ncol=nsim)
deerOutN2 <- matrix(nrow=TT, ncol=nsim)
deerOutN3plus <- matrix(nrow=TT, ncol=nsim)
deerInN <- matrix(nrow=TT, ncol=nsim) # total deer population in 
deerOutN <- matrix(nrow=TT, ncol=nsim) # total deer population out 
deerN <- matrix(nrow=TT, ncol=nsim) # total deer
deerpropIn <- matrix(nrow=TT, ncol=nsim) # proportion of deer from outside that move into the SNP during hunting

chamoisInN0 <- matrix(nrow=TT, ncol=nsim)
chamoisInN1 <- matrix(nrow=TT, ncol=nsim)
chamoisInN2 <- matrix(nrow=TT, ncol=nsim)
chamoisInN3plus <- matrix(nrow=TT, ncol=nsim)
chamoisInN <- matrix(nrow=TT, ncol=nsim)
chamoisOutN0 <- matrix(nrow=TT, ncol=nsim)
chamoisOutN1 <- matrix(nrow=TT, ncol=nsim)
chamoisOutN2 <- matrix(nrow=TT, ncol=nsim)
chamoisOutN3plus <- matrix(nrow=TT, ncol=nsim)
chamoisOutN <- matrix(nrow=TT, ncol=nsim)
chamoisN <- matrix(nrow=TT, ncol=nsim)


for(r in 1:nsim){
# starting populations (TODO: stable age composition)
  # TO DISCUSS: do we need to separate males and females in which species from which age on?
  wolfN0[1,r] <- 9
  wolfN1[1,r] <- 8
  wolfN2[1,r] <- 2
  wolfN3plus[1,r] <- 2
  wolfN[1, r] <- wolfN0[1,r]  + wolfN1[1,r] + wolfN2[1,r] +wolfN3plus[1,r]
  wolfP[1,r] <- 1  # number of packs
  
  deerInN0[1,r] <- 300
  deerInN1[1,r] <- 300
  deerInN1[1,r] <- 300
  deerInN2[1,r] <- 150
  deerInN3plus[1,r] <- 1000
  deerInN[1, r] <- deerInN0[1,r]  + deerInN1[1,r] + deerInN2[1,r] +deerInN3plus[1,r]
  
  deerOutN0[1,r] <- 500
  deerOutN1[1,r] <- 500
  deerOutN2[1,r] <- 250
  deerOutN3plus[1,r] <- 1500
  deerOutN[1, r] <- deerOutN0[1,r] + deerOutN1[1,r] + deerOutN2[1,r] +deerOutN3plus[1,r]
  deerN[1,r] <- deerInN[1,r] + deerOutN[1,r]
  
  chamoisInN0[1,r] <- 100
  chamoisInN1[1,r] <- 100
  chamoisInN2[1,r] <- 50
  chamoisInN3plus[1,r] <- 1000
  chamoisInN[1, r] <- chamoisInN0[1,r] + chamoisInN1[1,r] + chamoisInN2[1,r] +chamoisInN3plus[1,r]
  chamoisOutN0[1,r] <- 400
  chamoisOutN1[1,r] <- 400
  chamoisOutN2[1,r] <- 200
  chamoisOutN3plus[1,r] <- 2400
  chamoisOutN[1, r] <-  chamoisOutN0[1,r] +chamoisOutN1[1,r] + chamoisOutN2[1,r] +chamoisOutN3plus[1,r]
  chamoisN[1, r] <- chamoisInN[1,r] + chamoisOutN[1,r]
  
  
  for(t in 2:TT){
    
  # Wolf 
  wolfS0 <- FUNwolfS0((deerN[t-1,r]+chamoisN[t-1,r])/ungulatesmaxN)
  wolfS1 <- FUNwolfS1((deerN[t-1,r]+chamoisN[t-1,r])/ungulatesmaxN)
  wolfS2plus <- FUNwolfS2plus((deerN[t-1,r]+chamoisN[t-1,r])/ungulatesmaxN)
  
  wolfN1[t,r] <- rbinom(1, size=wolfN0[t-1,r] , prob=wolfS0)
  wolfN2[t,r] <- rbinom(1, size=wolfN1[t-1,r], prob=wolfS1)
  wolfN3plus[t,r] <- rbinom(1, size=wolfN2[t-1,r]+wolfN3plus[t-1,r], prob=wolfS2plus)
  wolfP[t,r] <- 2 # pack building
  
  wolfnrepfemales <- min(c(wolfN3plus[t,r]/2, wolfP[t,r]*2)) # number of reproducing females
  wolfN0[t,r] <- FUNwolfnpups(wolfnrepfemales) 
  wolfN[t, r] <-  wolfN0[t,r] + wolfN1[t,r] + wolfN2[t,r] +wolfN3plus[t,r]

  
  # Red deer
  propdeerInindiet <- deerInN[t,r]/(deerInN[t-1,r]+chamoisInN[t-1,r])
  if(is.na(propdeerInindiet)) propdeerInindiet <- 0
  propdeerOutindiet <- deerOutN[t-1,r]/(deerOutN[t-1,r]+chamoisOutN[t-1,r])
  if(is.na(propdeerOutindiet)) propdeerOutindiet <- 0
  
  deerInS0 <- FUNdeerInS0(wolfN[t-1,r]/wolfmaxN, propindiet=propdeerInindiet)
  deerInS1 <- FUNdeerInS1(wolfN[t-1,r]/wolfmaxN, propindiet=propdeerInindiet)
  deerInS2plus <- FUNdeerInS2plus(wolfN[t-1,r]/wolfmaxN, propindiet=propdeerInindiet)
  
  deerOutS0 <- FUNdeerInS0(wolfN[t-1,r]/wolfmaxN, propindiet=propdeerOutindiet)
  deerOutS1 <- FUNdeerInS1(wolfN[t-1,r]/wolfmaxN, propindiet=propdeerOutindiet)
  deerOutS2plus <- FUNdeerInS2plus(wolfN[t-1,r]/wolfmaxN, propindiet=propdeerOutindiet)
  
  
  deerInN1[t,r] <- rbinom(1, size=deerInN0[t-1,r], prob=deerInS0)
  deerInN2[t,r] <- rbinom(1, size=deerInN1[t-1,r], prob=deerInS1)
  deerInN3plus[t,r] <- rbinom(1, size=deerInN2[t-1,r]+deerInN3plus[t-1,r], prob=deerInS2plus)
  deerOutN1[t,r] <- rbinom(1, size=deerOutN0[t-1,r], prob=deerOutS0)
  deerOutN2[t,r] <- rbinom(1, size=deerOutN1[t-1,r], prob=deerOutS1)
  deerOutN3plus[t,r] <- rbinom(1, size=deerOutN2[t-1,r]+deerOutN3plus[t-1,r], prob=deerOutS2plus)

   deerInnrepfemales <- round((deerInN2[t-1,r]+deerInN3plus[t-1,r])/2) # number of reproducing females
  deerOutnrepfemales <- round((deerOutN2[t-1,r]+deerOutN3plus[t-1,r])/2) # number of reproducing females
   deerInN0[t,r] <- FUNdeerncalves(deerInnrepfemales, npop=deerInN[t-1,r], capacity=2000)
   deerOutN0[t,r] <- FUNdeerncalves(deerOutnrepfemales, npop=deerOutN[t-1,r], capacity=5000)
  
   deerInN[t, r] <- deerInN0[t,r] + deerInN1[t,r] + deerInN2[t,r] +deerInN3plus[t,r]
   deerOutN[t, r] <- deerOutN0[t,r] + deerOutN1[t,r] + deerOutN2[t,r] +deerOutN3plus[t,r]
   deerN[t, r] <- deerInN[t,r] + deerOutN[t,r]
    
  # Chamois
  
  propchamoisInindiet <- chamoisInN[t-1,r]/(deerInN[t-1,r]+chamoisInN[t-1,r])
  if(is.na(propchamoisInindiet)) propchamoisInindiet <- 0
  propchamoisOutindiet <- chamoisOutN[t-1,r]/(deerOutN[t-1,r]+chamoisOutN[t-1,r])
  if(is.na(propchamoisOutindiet)) propchamoisOutindiet <- 0

  chamoisInS0 <- FUNchamoisS0(wolfN[t-1,r]/wolfmaxN, propindiet=propchamoisInindiet)
  chamoisInS1 <- FUNchamoisS1(wolfN[t-1,r]/wolfmaxN, propindiet=propchamoisInindiet)
  chamoisInS2plus <- FUNchamoisS2plus(wolfN[t-1,r]/wolfmaxN, propindiet=propchamoisInindiet)

  chamoisOutS0 <- FUNchamoisS0(wolfN[t-1,r]/wolfmaxN, propindiet=propchamoisOutindiet)
  chamoisOutS1 <- FUNchamoisS1(wolfN[t-1,r]/wolfmaxN, propindiet=propchamoisOutindiet)
  chamoisOutS2plus <- FUNchamoisS2plus(wolfN[t-1,r]/wolfmaxN, propindiet=propchamoisOutindiet)
  
  chamoisInN1[t,r] <- rbinom(1, size=chamoisInN0[t-1,r], prob=chamoisInS0)
  chamoisInN2[t,r] <- rbinom(1, size=chamoisInN1[t-1,r], prob=chamoisInS1)
  chamoisInN3plus[t,r] <- rbinom(1, size=chamoisInN2[t-1,r]+chamoisInN3plus[t-1,r], prob=chamoisInS2plus)
  
  chamoisOutN1[t,r] <- rbinom(1, size= chamoisOutN0[t-1,r], prob=chamoisOutS0)
  chamoisOutN2[t,r] <- rbinom(1, size=chamoisOutN1[t-1,r], prob=chamoisOutS1)
  chamoisOutN3plus[t,r] <- rbinom(1, size=chamoisOutN2[t-1,r]+chamoisOutN3plus[t-1,r], prob=chamoisOutS2plus)
  
   chamoisInnrepfemales <- round((chamoisInN2[t-1,r]+chamoisInN3plus[t-1,r])/2) # number of reproducing females
  chamoisInN0[t,r] <- FUNchamoisncalves(chamoisInnrepfemales, npop=chamoisInN[t-1,r], capacity=2000)
  chamoisOutnrepfemales <- round((chamoisOutN2[t-1,r]+chamoisOutN3plus[t-1,r])/2) # number of reproducing females
  chamoisOutN0[t,r] <- FUNchamoisncalves(chamoisOutnrepfemales, npop=chamoisOutN[t-1,r], capacity=8000)
 
  chamoisInN[t, r] <- chamoisInN0[t,r] + chamoisInN1[t,r] + chamoisInN2[t,r] +chamoisInN3plus[t,r]
  chamoisOutN[t, r] <- chamoisInN0[t,r] + chamoisInN1[t,r] + chamoisOutN2[t,r] +chamoisOutN3plus[t,r]
  chamoisN[t, r] <- chamoisInN[t,r] + chamoisOutN[t,r]
  } # close t
} # close r

7.2 Results

ATTENTION: This model is under construction, thus the results here show how a possible outcome could look like. However, the parameters used to simulate the model are not yet worked out properly. Therefore, be aware: unreliable input gives unreliable output.

par(mfrow=c(7,1), mar=c(1,4,0.1,1), oma=c(3,2,0,0))
plot(1:TT, seq(0, max(wolfN), length=TT), type="n", xlab=NA, ylab=NA, xaxt="n")
for(i in 1:nsim) lines(1:TT, wolfN[,i], col=rgb(0,0,0,0.2))
text(1,max(wolfN), adj=c(0,1) , labels="Wolf region")
mtext("Population size",side=2, line=0, outer=TRUE)


plot(1:TT, seq(0,max(deerN), length=TT), type="n", xlab="year", ylab="", xaxt="n")
for(i in 1:nsim) lines(1:TT, deerN[,i], col=rgb(0,0,0,0.2), xaxt="n")
text(1,max(deerN), adj=c(0,1) , labels="Red deer region")

plot(1:TT, seq(0,max(deerInN), length=TT), type="n", xlab="year", ylab="", xaxt="n")
for(i in 1:nsim) lines(1:TT, deerInN[,i], col=rgb(0,1,0,0.2), xaxt="n")
text(1,max(deerInN), adj=c(0,1) , labels="Red deer in park")


plot(1:TT, seq(0,max(deerOutN), length=TT), type="n", xlab="year", ylab="", xaxt="n")
for(i in 1:nsim) lines(1:TT, deerOutN[,i], col=rgb(0,0,1,0.2), xaxt="n")
text(1,max(deerOutN), adj=c(0,1) , labels="Red deer outside park")


plot(1:TT, seq(0,max(chamoisN), length=TT), type="n", xlab="year", ylab="", xaxt="n")
for(i in 1:nsim) lines(1:TT, chamoisN[,i], col=rgb(0,0,0,0.2))
text(1,max(chamoisN), adj=c(0,1) , labels="Chamois region")


plot(1:TT, seq(0,max(chamoisInN), length=TT), type="n", xlab="year", ylab="", xaxt="n")
for(i in 1:nsim) lines(1:TT, chamoisInN[,i], col=rgb(0,1,0,0.2))
text(1,max(chamoisInN), adj=c(0,1) , labels="Chamois in park")

plot(1:TT, seq(0,max(chamoisOutN), length=TT), type="n", xlab="year", ylab="")
for(i in 1:nsim) lines(1:TT, chamoisOutN[,i], col=rgb(0,0,1,0.2))
text(1,max(chamoisOutN), adj=c(0,1) , labels="Chamois outside park")