/
parasitoidhost.rmd
237 lines (188 loc) · 8.17 KB
/
parasitoidhost.rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
---
title: "Parasitoid-host models with embedded ShinyApps"
author: Ottar N. Bjørnstad
output: html_document
runtime: shiny
---
Version 0.5-8 August 27, 2022
https://github.com/objornstad/
This Rmarkdown of Nicholson-Bailey (Nicholson and Bailey 1935) and negative-binomial (May 1978) parasitoid-host model was written by Ottar N. Bjørnstad and is released with a CC-BY-NC license for anyone to improve and re-share (acknowledging origin). Please email me a copy if you make an update (onb1 at psu dot edu). The app was originally developed as part of the epimdr-package (https://cran.r-project.org/package=epimdr; Bjørnstad 2019).
The app requires the shiny package to be installed to run. "Run document" in Rstudio will launch the App.
```{r, include=FALSE}
using<-function(...) {
libs<-unlist(list(...))
req<-unlist(lapply(libs,require,character.only=TRUE))
need<-libs[req==FALSE]
if(length(need)>0){
install.packages(need)
lapply(need,require,character.only=TRUE)
}
}
using("shiny")
```
The Nicholson-Bailey (1935) model assumes a population of randomly searching parasitoids, $P$, parasitizing immature stages of their host, $H$. The host growth rate is $R$ and the parasitoid search efficiency is $a$. Given the random search, the fraction of hosts escaping parasitization is $\mbox{exp}(- a P)$ leading to two equations:
$\begin{aligned}
H_{t+1} & = \underbrace{R H_t}_{\mbox{reproduction}} \underbrace{e^{-a P_t}}_{\mbox{not parasitized}}\\
P_{t+1} & = R H_t \underbrace{(1-e^{-a P_t}) }_{\mbox{parasitized}}
\end{aligned}$
Nicholson and Bailey (1935) showed that the equilibrium $\{H^* = \frac{\mbox{log}(R)}{a (R-1)}, P^* = \frac{\mbox{log}(R)}{a}\}$" is an unstable focus leading to increasingly violent consumer-resource cycles.
The Nicholson-Bailey shiny app is:
```{r, echo=FALSE}
# This creates the NB User Interface (UI)
ui = pageWithSidebar(
headerPanel("Nicholson-Bailey Model"),
sidebarPanel(
sliderInput("R", "Growth rate (R):", 1.1,
min = 1, max = 2, step=.01),
sliderInput("a", "Search efficiency (a):", 0.03,
min = 0, max = .5),
numericInput("P0", "Initial parasitoid:", 10,
min = 1, max = 100),
numericInput("H0", "Initial host:", 20,
min = 1, max = 100),
numericInput("Tmax", "Tmax:", 100,
min = 1, max = 500)
),
mainPanel(tabsetPanel(
tabPanel("Simulation", plotOutput("plot1", height = 500)),
tabPanel("Phase plane", plotOutput("plot2", height = 500)),
tabPanel("Details",
withMathJax(
helpText("MODEL:"),
helpText("Host $$H_t = R H_{t-1} (1 - \\mbox{exp}(- a P_{t-1}))$$"),
helpText("Parasitoid $$P_t = R H_{t-1} \\mbox{exp}(- a P_{t-1})$$"),
helpText("Equilibria $$H^* = \\frac{\\mbox{log}(R)}{a (R-1)},
P^* = \\frac{\\mbox{log}(R)}{a}$$")),
helpText("REFERENCE: Nicholson AJ, Bailey VA (1935) The balance of animal populations.
Proceedings of the Zoological Society of London 3: 551-598")
))
)
)
# This creates the 'behind the scenes' code (Server)
server = function(input, output) {
NB = function(R, a, T = 100, H0 = 10, P0 = 1){
#T is length of simulation (number of time-steps)
#H0 and P0 are initial numbers
#we provide default parameters except for R and a
H=rep(NA, T) #host series
P=rep(NA, T) #parasitoid series
H[1] = H0 #Initiating the host series
P[1] = P0 #Initiating the host series
for(t in 2:T){
H[t] = R * H[t-1] * exp(- a * P[t-1])
P[t] = R * H[t-1] * (1-exp(- a * P[t-1]))
if(P[t-1]==0) break
} #end of loop
#the two vectors of results are stored in a "list"
res= data.frame(H = H, P = P)
#the list is passed out of this function
return(res)
} #end of function
output$plot1 <- renderPlot({
sim= NB(R=input$R, a=input$a, H0=input$H0, P0=input$P0, T=input$Tmax)
time = 1:input$Tmax
plot(time, sim$H, type= "b",xlab = "Generations", ylab = "Abundance",
ylim = range(sim, na.rm=TRUE))
points(time, sim$P, type = "b", col = "red")
legend("topleft",
legend=c("H", "P"),
lty=c(1,1),
pch=c(1,1),
col=c("black", "red"))
})
output$plot2 <- renderPlot({
sim= NB(R=input$R, a=input$a, H0=input$H0, P0=input$P0, T=input$Tmax)
time = 1:input$Tmax
Hstar=log(input$R)/(input$a*(input$R-1))
Pstar=log(input$R)/input$a
plot(sim$H, sim$P, type= "b",xlab = "Host", ylab = "Parasitoid")
points(Hstar, Pstar, col=2, pch=19)
})
}
shinyApp(ui, server, options = list(height = 640))
```
May (1978) proposed that heterogeneities in searching can stabilize parasitoid host dynamics. Assuming a negative binomial (rather than Poisson) attack probability the model is:
$\begin{aligned}
H_{t+1} & = R H_t (1+\frac{a P_t}{k})^{-k}\\
P_{t+1} & = R H_t (1-(1+\frac{a P_t}{k})^{-k}),
\end{aligned}$
where $k$ is the aggregation parameter ($k \rightarrow \infty$, Neg Bin $\rightarrow$ Poisson). He coined the CV$^2$-rule which says that if the coefficient-of-variation in attack rate is greater than 1, the parasitoid-host dynamics stabilizes; The CV for the negative binomial being $1/\sqrt{k}$. The equilibrium $P^* = k (R^{1/k}-1) /a, H^*=P^* R / (R-1)$ is either stable or unstable depending on whether k is smaller or greater than 1.
```{r, echo=FALSE}
# This creates the User Interface (UI)
ui = pageWithSidebar(
headerPanel("May's Parasitoid-host Model"),
sidebarPanel(
sliderInput("R", "Growth rate (R):", 1.1,
min = 1, max = 2, step=.01),
sliderInput("a", "Search efficiency (a):", 0.1,
min = 0, max = .5),
sliderInput("k", "aggregation (k):", 1.5,
min = 0.1, max = 3, step=0.1),
numericInput("P0", "Initial parasitoid:", 10,
min = 1, max = 100),
numericInput("H0", "Initial host:", 20,
min = 1, max = 100),
numericInput("Tmax", "Tmax:", 100,
min = 1, max = 500)
),
mainPanel(tabsetPanel(
tabPanel("Simulation", plotOutput("plot1", height = 500)),
tabPanel("Phase plane", plotOutput("plot2", height = 500)),
tabPanel("Details",
withMathJax(
helpText("MODEL:"),
helpText("Host $$H_t = R H_{t-1} (1 + a P_{t-1})^k$$"),
helpText("Parasitoid $$P_t = R H_{t-1} (1-(1 + a P_{t-1})^k)$$"),
helpText("REFERENCE: May RM (1978) Host-parasitoid systems in patchy
environments: a phenomenological model. J Anim Ecol 47: 833-843")
)
)
)
)
)
server = function(input, output) {
NB = function(R, a, k, T = 100, H0 = 10, P0 = 1){
#T is length of simulation (number of time-steps)
#H0 and P0 are initial numbers
#we provide default parameters except for R and a
H=rep(NA, T) #host series
P=rep(NA, T) #parasitoid series
H[1] = H0 #Initiating the host series
P[1] = P0 #Initiating the host series
for(t in 2:T){
H[t] = R * H[t-1] * (1+ a * P[t-1]/k)^(-k)
P[t] = R * H[t-1] * (1-(1+ a * P[t-1]/k)^(-k))
if(P[t-1]==0) break
} #end of loop
#the two vectors of results are stored in a "list"
res= data.frame(H = H, P = P)
#the list is passed out of this function
return(res)
} #end of function
output$plot1 <- renderPlot({
sim= NB(R=input$R, a=input$a, k=input$k, H0=input$H0, P0=input$P0, T=input$Tmax)
time = 1:input$Tmax
plot(time, sim$H, type= "b",xlab = "Generations", ylab = "Abundance",
ylim = range(sim, na.rm=TRUE))
points(time, sim$P, type = "b", col = "red")
legend("topleft",
legend=c("H", "P"),
lty=c(1,1),
pch=c(1,1),
col=c("black", "red"))
})
output$plot2 <- renderPlot({
sim= NB(R=input$R, a=input$a, k=input$k, H0=input$H0, P0=input$P0, T=input$Tmax)
time = 1:input$Tmax
plot(sim$H, sim$P, type= "b",xlab = "Host", ylab = "Parasitoid")
#Pstar=input$k*(input$R^(1/input$k)-1)/input$a
#Hstar=Pstar*input$R/(input$R-1)
#points(Hstar, Pstar, col=2, pch=19)
})
}
shinyApp(ui, server, options = list(height = 660))
```
References:
Bjørnstad, O.N. (2018) Epidemics: models and data using R. Springer.
May, R.M. (1978) Host-parasitoid systems in patchy environments: a phenomenological model. J Anim Ecol 47: 833-843
Nicholson A.J. and Bailey V.A. (1935) The balance of animal populations. Proceedings of the Zoological Society of London 3: 551-598.