-
Notifications
You must be signed in to change notification settings - Fork 16
/
temporal-models.qmd
246 lines (177 loc) · 10.2 KB
/
temporal-models.qmd
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
238
239
240
241
242
243
244
245
246
# Population models
Mathematical models can be fitted to the DPC data to express epidemic progress in terms of rates and absolute/relative quantities. The latter can be accomplished using population dynamics (or growth-curve) models for which the estimated parameters are usually meaningful biologically and appropriately describe epidemics that do not decrease in disease intensity. By fitting an appropriate model to the progress curve data, another set of parameters is available to the researcher when attempting to represent, understand or compare epidemics.
The family of models that describe the growth of epidemics, hence population dynamics model, are known as deterministic models of continuous time [@chapter2007b]. These models are usually fitted to DPC data to obtain two or more biologically meaningful parameters. Here, these models and their formulations are shown using R scripts to simulate the theoretical curves for each model.
## Non-flexible models
These population dynamics models require at least two parameters, hence they are known as non-flexible, as opposed to the flexible ones for which there are at least one additional (third) parameter.
Following the convention proposed by [@chapter2007b] in their book "The study of plant disease epidemics":
- time is represented by $t$
- disease intensity by $y$
- the rate of change in $y$ between two time units is represented by $\frac{dy}{dt}$
Now we can proceed and learn which non-flexible models exist and for which situation they are more appropriate.
### Exponential
The differential equation for the exponential model is given by
$\frac{dy}{dt} = r_E.y$,
where $r_E$ is the apparent infection rate (subscript E for this model) (sensu Vanderplank) and $y$ is the disease intensity. Biologically, this formulation suggests that diseased plants, or $y$, and $r_E$ at each time contribute to disease increase. The value of $\frac{dy}{dt}$ is minimal when $y = 0$ and increases exponentially with the increase in $y$.
The integral for the exponential model is given by
$y = y_0 e^{r_Et}$,
where $y0$ is and $r$ are obtained via estimation. Let's simulate two curves by varying $r$ while fixing $y0$ and varying the latter while fixing $r_E$. We produce the two plots in *ggplot* and add the predicted curve using the \`stat_function\`. But first, we need to define values for the two model parameters. Further modifications to these values will be handled directly in the simulation (e.g. doubling infection rate, reducing initial inoculum by half, etc.).
```{r}
#| warning: false
#| message: false
library(tidyverse) # essential packages
library(cowplot)
library(r4pde)
theme_set(theme_r4pde()) # set global theme
```
```{r}
y0 <- 0.001
r <- 0.06
tmax <- 60 # maximum duration t of the epidemics
dat <- data.frame(t = seq(1:tmax), y = seq(0:1)) # define the axes
```
In the plot below, note that the infection rate in one curve was doubled ($r$ = 0.12)
```{r}
#| warning: false
#| label: fig-exp1
#| fig-cap: "Exponential curves with two rates of infection (0.06 and 0.12) and the same initial inoculum (0.001) "
dat |>
ggplot(aes(t, y)) +
stat_function(fun = function(t) y0 * exp(r * t), linetype = 1) +
stat_function(fun = function(t) y0 * exp(r * 2 * t), linetype = 2) +
ylim(0, 1) +
theme_r4pde()+
labs(x = "Time")
```
Now the inoculum was increased five times while using the same doubled rate.
```{r}
#| warning: false
#| label: fig-exp2
#| fig-cap: "Exponential curves with the same rate of infection (0.12) and single and five times the initial inoculum (0.001)"
dat |>
ggplot(aes(t, y)) +
stat_function(fun = function(t) y0 * exp(r * 2 * t), linetype = 1) +
stat_function(fun = function(t) y0 * 5 * exp(r * 2 * t), linetype = 2) +
ylim(0, 1) +
theme_r4pde()+
labs(x = "Time")
```
### Monomolecular
The differential of the monomolecular model is given by
$\frac{dy}{dt} = r_M (1-y)$
where now the $r_M$ is the rate parameter of the monomolecular model and $(1-y)$ is the proportion of non-infected (healthy) individuals or host tissue. Note that $\frac{dy}{dt}$ is maximum when $y = 0$ and decreases when $y$ approaches 1. Its decline is due to decrease in the proportion of individuals or healthy sites with the increase in $y$. Any inoculum capable of infecting the host will more likely land on infected individuals or sites.
The integral of the monomolecular model is given by
$\frac{dy}{dt} = 1 - (1-y)e^{-r_Mt}$
This model commonly describes the temporal patterns of the monocyclic epidemics. In those, the inoculum produced during the course of the epidemics do not contribute new infections. Therefore, different from the exponential model, disease intensity $y$ does not affect the epidemics and so the absolute rate is proportional to $(1-y)$.
Let's simulate two monomolecular curve with different rate parameters where one is one third of the other.
```{r}
#| warning: false
#| label: fig-mono1
#| fig-cap: "Monomolecular curves with two rates of infection (0.06 and 0.02) and the same initial inoculum (0.001) "
dat |>
ggplot(aes(t, y)) +
stat_function(fun = function(t) 1 - ((1 - y0) * exp(-r * t))) +
stat_function(fun = function(t) 1 - ((1 - y0) * exp(-(r / 3) * t))) +
labs(x = "Time") +
theme_r4pde()+
annotate(geom = "text", x = 35, y = 0.77, label = "r = 0.06") +
annotate(geom = "text", x = 50, y = 0.55, label = "r = 0.02")
```
Now inoculum was increased 100 times with the reduced rate.
```{r}
#| warning: false
#| label: fig-mono2
#| fig-cap: "Monomolecular curves with one rate (0.06) and the initial inoculum increased 100 times"
dat |>
ggplot(aes(t, y)) +
stat_function(fun = function(t) 1 - ((1 - y0) * exp(-r / 2 * t))) +
stat_function(fun = function(t) 1 - ((1 - (y0 * 100)) * exp(-r / 2 * t))) +
theme_r4pde()+
labs(x = "Time") +
annotate(geom = "text", x = 35, y = 0.77, label = "y0 = 0.1") +
annotate(geom = "text", x = 45, y = 0.65, label = "y0 = 0.001")
```
### Logistic
The logistic model is a more elaborated version of the two previous models as it incorporates the features of them both. Its differential is given by
$\frac{dy}{dt} = r_L. y . (1 - y)$,
where $r_L$ is the infection rate of the logistic model, $y$ is the proportion of diseased individuals or host tissue and $(1-y)$ is the proportion of non-affected individuals or host area.
Biologically, $y$ in its differential equation implies that $\frac{dy}{dt}$ increases with the increase in $y$ (as in the exponential) because more disease means more inoculum. However, $(1-y)$ leads to a decrease in $\frac{dy}{dt}$ when $y$ approaches the maximum $y=1$, because the proportion of healthy individuals or host area decreases (as in the monomolecular). Therefore, $\frac{dy}{dt}$ is minimal at the onset of the epidemics, reaches a maximum when $y/2$ and declines until $y=1$.
The integral is given by
$y = \frac{1}{1 + (1-y_0).e^{-r.t}}$,
where $r_L$ is the apparent infection rate of the logistic model and $y0$ is the disease intensity at $t=0$. This model provides a good fit to polycyclic epidemics.
Let's check two curves where in one the infection rate is double while keeping the same initial inoculum.
```{r}
#| warning: false
#| label: fig-log1
#| fig-cap: "Logistic curves with two rates of infection (0.18 and 0.024) and the same initial inoculum (0.001) "
dat |>
ggplot(aes(t, y)) +
stat_function(
linetype = 2,
fun = function(t) 1 / (1 + ((1 - y0) / y0) * exp(-r * 2 * t))
) +
stat_function(fun = function(t) 1 / (1 + ((1 - y0) / y0) * exp(-r * 4 * t))) +
labs(x = "Time") +
theme_r4pde()+
annotate(geom = "text", x = 41, y = 0.77, label = "r = 0.18") +
annotate(geom = "text", x = 50, y = 0.10, label = "r = 0.024")
```
Now the inoculum is reduced 10 times for a same infection rate.
```{r}
#| warning: false
#| label: fig-log2
#| fig-cap: "Logistic curves with a single rate of infection (0.24) and two initial inoculum (0.001 and 0.0001) "
dat |>
ggplot(aes(t, y)) +
stat_function(
linetype = 2,
fun = function(t) 1 / (1 + ((1 - (y0 / 10)) / (y0 / 10)) * exp(-r * 3 * t))
) +
stat_function(fun = function(t) 1 / (1 + ((1 - y0) / y0) * exp(-r * 3 * t))) +
labs(x = "Time") +
theme_r4pde()+
annotate(geom = "text", x = 35, y = 0.77, label = "y0 = 0.001") +
annotate(geom = "text", x = 50, y = 0.10, label = "y0 = 0.0001")
```
### Gompertz
The Gompertz model is similar to the logistic and also provides a very good fit to several polycyclic diseases. The differential equation is given by
$\frac{dy}{dt} = r_G.[ln(1) - ln(y)]$
Differently from the logistic, the variable representing the non-infected individuals or host area is $-ln(y)$. The integral equation is given by
$y = e^{(ln(y0)).{e^{-r_G.t)}}}$,
where $r_G$ is the apparent infection rate for the Gompertz models and $y_0$ is the disease intensity at $t = 0$.
Let's check curves for two rates.
```{r}
#| warning: false
#| label: fig-gomp1
#| fig-cap: "Gompertz curves with two rates of infection (0.12 and 0.03) and the same initial inoculum (0.001)"
dat |>
ggplot(aes(t, y)) +
stat_function(
linetype = 2,
fun = function(t) exp(log(y0) * exp(-r/2 * t))
) +
stat_function(fun = function(t) exp(log(y0) * exp(-r*2 * t))) +
labs(x = "Time") +
theme_r4pde()+
annotate(geom = "text", x = 41, y = 0.77, label = "r = 0.12") +
annotate(geom = "text", x = 50, y = 0.10, label = "r = 0.03")
```
And those when inoculum was reduced one thousand times.
```{r}
#| warning: false
#| label: fig-gomp2
#| fig-cap: "Gompertz curves with a single rate of infection (0.12) and two levels of initial inoculum (0.001 and 0.00001)"
dat |>
ggplot(aes(t, y)) +
stat_function(
linetype = 2,
fun = function(t) exp(log(y0) * exp(-r*2 * t))
) +
stat_function(fun = function(t) exp(log(y0/1000) * exp(-r*2 * t))) +
labs(x = "Time") +
theme_r4pde()+
annotate(geom = "text", x = 15, y = 0.77, label = "y0 = 0.001") +
annotate(geom = "text", x = 25, y = 0.10, label = "y0 = 0.00001")
```
## Interactive application
A shiny app was developed to demonstrate these four models interactively. Click on the image below to get access to the app.
[{#fig-models fig-align="center" width="533"}](https://edelponte.shinyapps.io/epidemics/)