In noamross/age-infects: Package for manuscript and code of age-infection interaction paper

% The null model for age effects with overdispersed infectio and a reservoir host % Noam Ross % 13-06-11 17:32:30

Now I add a reservoir host $B$ to the multi-infection model. In Sudden Oak Death systems, Bay Laurel plays this role: it is infected and broadcasts new spores but suffers no mortality effects. For completeness, I include non-hosts $R$, which compete for space but are not infected, too, but for now I'll keep their population at zero. The resevoir and inert hosts have no age structure for now, though there may be reasons to include it in the future.

$$\begin{aligned} \frac{dJ}{dt} &= A f_A \left(1 - \frac{J+A+B+R}{K} \right) + J \left(f_J \left(1 - \frac{J+A}{K} \right) - d_J - g\right) - \alpha P_J \ \frac{dA}{dt} &= J g - A d_A - \alpha P_A \ \frac{dB}{dt} &= B f_B \left(1 - \frac{J+A+B+R}{K}\right) - B d_B \ \frac{dR}{dt} &= R f_R \left(1 - \frac{J+A+B+R}{K}\right) - R d_R \ \frac{dP_J}{dt} &= \lambda \frac{J}{K} (P_J + P_A + P_R) - P_J \left(d_J + \mu + g + \alpha \left(1 + \frac{(k+1)P_J}{kJ} \right) \right) \ \frac{dP_A}{dt} &= \lambda \frac{A}{K} (P_J + P_A + P_R) + P_J g - P_A \left(d_A + \mu + \alpha \left(1 + \frac{(k+1)P_A}{kA} \right) \right) \ \frac{dP_B}{dt} &= \lambda \frac{B}{K} (P_J + P_A + P_R) - P_A \left(d_B + \mu\right) \end{aligned}$$

r opts_chunk$set(tidy=FALSE, warning=FALSE, message=FALSE, fig.path=paste(gsub(".Rmd", "", basename(knitr:::knit_concord$get('infile')) ),"-", sep=""))

od.model <- function(t, y, parms) {
  list2env(as.list(y), environment())
  list2env(as.list(parms), environment())
  dJ <- A*f_a*(1 - (J+A+B+R)/K) + J * (f_j*(1 - (J+A)/K) - d_j - g) - a_j * PJ
  dA <- J*g - (A * d_a) - a_a * PA
  dB <- B*f_b*(1 - (J+A+B+R)/K) - B * d_b
  dR <- R*f_r*(1 - (J+A+B+R)/K) - R * d_r
  dPJ <- lambda * (PJ + PA + PB) * J/K - 
         PJ*(d_j + mu + g + a_j*(1 + (k+1)*PJ/(k*J)))
  dPA <- lambda * (PJ + PA + PB) * A/K + PJ * g - 
         PA*(d_a + mu + a_a*(1 + (k+1)*PA/(k*A)))
  dPB <- lambda *(PJ + PA + PB) * B/K - PB*(d_b + mu)
  dPR <- 0
  return(list(c(dJ=dJ, dA=dA, dB=dB, dR=dR,
                dPJ=dPJ, dPA=dPA, dPB=dPB, dPR=dPR)))
}

fec <- c(f_j=0.01,
         f_a=0.01,
         f_b=0.01,
         f_r=0.01)

mort <- c(d_j=0.005,
          d_a=0.005,
          d_b=0.005,
          d_r=0.005)

alpha <- c(a_j=0.05,
          a_a=0.05,
          a_b=0,
          a_r=0)

parms <- c( 
  fec,
  mort,
  alpha,
  g=0.1,
  lambda=0.3,
  K=50,
  mu=0.00,
  k=0.000
  )

#A_ss = with(as.list(parms), K/(d_a/g + 1))
#J_ss = with(as.list(parms), K - A_ss)
init <- c(J=0.93792, A=18.75839, B=5.50823, R=0, PJ=0.094, PA=1.876, PB=0.55, PR=0)

require(deSolve)
require(reshape2)
require(plyr)
ks <- c(10000)
df <- adply(ks, 1, function(x) {
  parms["k"] <- x
  df <- data.frame(k.val=as.factor(x), 
                   as.data.frame(lsoda(y=init, times=1:100, func=od.model, 
                                       parms=parms)))
  })

names(df)[names(df)=="time"] <- "Time"
df$X1 <- NULL

df <- cbind(melt(df[,1:6], id.vars=c("k.val", "Time")),
      melt(df[,c(1,2,7,8,9,10)], id.vars=c("k.val", "Time"))[4])

names(df)[3:5]<- c("Class", "Population", "Infections")


#mor <- mort[-which(df$Class==zeroclass)]

df <- ddply(df, .(k.val, Time), summarize,
            Population=Population,
            Infections=Infections,
            Class=Class,
            FracPop = Population/sum(Population),
            InfectPerIndiv = ifelse(is.nan(Infections/Population), 0, Infections/Population),
            FracInfected = 1 - exp(-InfectPerIndiv),
            TotalInfected = rep(sum(Population*FracInfected),length(Population)),
            MortalityRate = mort + alpha * Infections / (Population * FracInfected),
            YearsToMortality = 1/MortalityRate,
            InfectionRate = 1 - exp(-parms["lambda"] * sum(Infections) * 
                                    Population * exp(-InfectPerIndiv) / 
                                    parms["K"]),
            YearsToInfection = 1/InfectionRate
            )

check.df <- ddply(df, .(k.val, Class), summarize, SUM=sum(Population))
zeroclass <- check.df$Class[which(check.df$SUM==0)]
df <- droplevels(df[-which(df$Class==zeroclass),])

require(ggplot2)
require(gridExtra)
theme_nr <- theme_nr + theme(legend.title=element_text(size=22),
                             legend.text=element_text(size=16),
                             legend.key.size=unit(.75, "cm"))
JAlab <- scale_color_discrete(labels=c("Small Tanoak","Big Tanoak", "Bay"), name="") 
p1 <- ggplot(df, aes(x=Time, y=Population, col=Class)) + 
      geom_line(lwd=1) + theme_nr + ylab("Population") + JAlab
p2 <- ggplot(df, aes(x=Time, y=FracPop, col=Class)) + 
      geom_line(lwd=1) + theme_nr + ylab("Fraction of Population") + JAlab
grid.arrange(p1, p2, nrow=2)

p3 <- ggplot(df, aes(x=Time, y=Infections, col=Class)) +
      geom_line(lwd=1) + theme_nr + ylab("Number of Infections") + JAlab
p4 <- ggplot(df, aes(x=Time, y=InfectPerIndiv, col=Class)) +
      geom_line(lwd=1) + theme_nr + ylab("Infections per Individual") + JAlab
grid.arrange(p3, p4, nrow=2)

ggplot(df, aes(x=Time, y=FracInfected, col=Class)) + 
geom_line(lwd=1) + theme_nr + ylab("Fraction infected") + JAlab

ggplot(subset(df, Class != "B"), aes(x=TotalInfected, y=YearsToMortality, col=Class)) + 
  geom_point(cex=4) + theme_nr + xlab("Number of infected trees") + 
  ylab("Years to death of infected individuals") + JAlab

ggplot(subset(df, Time < 70), aes(x=TotalInfected, y=YearsToInfection, col=Class)) + 
  geom_point(cex=4) + theme_nr + xlab("Number of infected trees") + 
  ylab("Years to Infection of Healthy Individuals") + JAlab

TODO:

Add resprouting.

noamross/age-infects documentation built on May 23, 2019, 9:30 p.m.