Middle of the epidemiological model

tpt_model2.qmd

@@ -32,7 +32,8 @@ You can ask report errors or for additional analysis [here](https://github.com/T
 
 ## Packages
 
-Required: `deSolve`, `tibble`, `purrr`, `dplyr`, `magrittr`, `bbmle`, `RColorBrewer`
+Required: `deSolve`, `tibble`, `purrr`, `dplyr`, `magrittr`, `bbmle`, `RColorBrewer`,
+`parallel`
 
 ```{r}
 library(tibble)
@@ -95,151 +96,18 @@ polygon2 <- function(x, y, col = 4, alpha = .2, ...) {
 }
 ```
 
-
-## Epidemiological models
-
-
-### A minimalist SID
-
-![Flow diagram of a minimalist SID model.](SIDbis.png){width=460 #fig-tinySID}
-
-#### Dynamics
-
-$$
-\begin{align} 
-  \frac{dS}{dt} &= \sigma I + \alpha D - (1 - \pi) \beta DS \\
-  \frac{dI}{dt} &= (1 - \pi) \beta DS - (\delta + \sigma) I   \\
-  \frac{dD}{dt} &= \delta I - \alpha D
-\end{align}
-$$
-
-R code:
-
 ```{r}
-sid_dyn <- function(N, I0 = 0, D0, beta, sigma, pi, delta, alpha, times) {
-  ode2(c(S = N - I0 - D0, I = I0, D = D0),
-       times,
-       function(time, state, pars) {
-         with(as.list(c(state, pars)), {
-           infection <- (1 - pi) * beta * D * S
-           dS  <- sigma * I + alpha * D - infection
-           dI  <- infection - (delta + sigma) * I
-           dD  <- delta * I - alpha * D
-           list(c(dS, dI, dD))
-         })
-       },
-       c(beta = beta, sigma = sigma, pi = pi, delta = delta, alpha = alpha))
+polygon3 <- function(x, y1, y2, col = 4, ...) {
+  polygon(c(x, rev(x)), c(y1, rev(y2)), col = col, border = NA, ...)
 }
 ```
 
-#### Equilibrium
-
-$$
-\begin{align} 
-  S^* &= N - I^* - D^*             \\
-  I^* &= \frac{\alpha}{\delta} D^* \\
-  D^* &= \frac{\delta}{\alpha + \delta} N -
-         \frac{\alpha}{\alpha + \delta} \frac{\delta + \sigma}{(1 - \pi)\beta}
-\end{align}
-$$
-
-R code:
-
 ```{r}
-sid_equ <- function(N, beta, sigma, pi, delta, alpha) {
-  Dstar <- delta * N / (alpha + delta) - alpha * (delta + sigma) /
-           ((alpha + delta) * (1 - pi) * beta)
-  Istar <- alpha * Dstar / delta
-  c(S = N - Istar - Dstar, I = Istar, D = Dstar)
+mclapply2 <- function(...) {
+  parallel::mclapply(..., mc.cores = parallel::detectCores() - 1)
 }
 ```
 
-A function that plots the dynamics and equilibrium values of the infection prevalence
-and the disease incidence and prevalence:
-
-```{r}
-plot_dyn_equ <- function(
-    N     = 1e5,
-    D0    = 350,
-    beta  = .000001,
-    sigma = 0,
-    pi    = 0,
-    delta = .0001,
-    alpha = 1 / 70,
-    times = seq2(0, 25000)) {
-  sims <- sid_dyn(N = N, I0 = 0, D0 = D0, beta = beta, sigma = sigma, pi = pi,
-                  delta = delta, alpha = alpha, times = times)
-
-  equs <- sid_equ(N = N, beta = beta, sigma = sigma, pi = pi,
-                  delta = delta, alpha = alpha)
-  
-  plotl2 <- function(...) plotl(..., xlab = "time (year)")
-  
-  opar <- par(mfrow = c(1, 3), cex = 1, plt = c(.25, .95, .25, .9))
-  with(sims, {
-    plotl2(time, I, col = 2, ylab = "infection prevalence")
-    abline(h = equs["I"], col = 2)
-    plotl2(time, D, col = 3, ylab = "disease prevalence")
-    abline(h = equs["D"], col = 3)
-    plotl2(time, delta * I, col = 4, ylab = "yearly disease incidence")
-    abline(h = delta * equs["I"], col = 4)
-  })
-  par(opar)
-}
-```
-
-Let's do some verification:
-
-```{r fig.height = 2.7}
-plot_dyn_equ()
-```
-
-
-#### Calibration
-
-In absence of TPT (i.e. $\pi = 0$), let's further assume that $\sigma = 0$. Let's call
-$d$ the yearly disease incidence. Then,
-
-$$
-\begin{align} 
-  \delta &= \frac{d^*}{I^*}             \\
-  \alpha &= \frac{D^*}{d^*}   \\
-  \beta  &= \frac{N}{\alpha} - \frac{\alpha + \delta}{\alpha \delta} D^*
-\end{align}
-$$
-
-Corresponding R code:
-
-```{r}
-param_val <- function(inf_prev, dis_prev, dis_incd, N = 1e5) {
-  delta <- dis_incd / inf_prev
-  alpha <- dis_incd / dis_prev
-  c(alpha = alpha,
-    beta  = alpha * delta / (delta * N - (alpha + delta) * dis_prev),
-    delta = delta)
-}
-```
-
-In Vietnam:
-
-* prevalence of infection: somewhere between 10 and 30% of the population?
-* prevalence of disease: 176 / 100,000
-* yearly incidence of disease: 322 / 100,000
-
-```{r}
-(p_val <- param_val(2e4, 176, 322))
-```
-
-Let's see:
-
-```{r fig.height = 2.7}
-plot_dyn_equ(D0    = 176,
-             beta  = p_val[["beta"]],
-             delta = p_val[["delta"]],
-             alpha = p_val[["alpha"]],
-             times = seq2(0, 3000))
-```
-
 
 ## Prophylaxis model
 
@@ -531,7 +399,7 @@ ts <- seq2(0, 30)
 plot4(ts, hill(ts, 10, 1, 1), ylim = 0:1, ylab = "treatment efficacy")
 ```
 
-## Full model
+### Full prophylaxis
 
 The $\pi$ proportion of the epidemiological model can now be expressed as:
 
@@ -561,53 +429,190 @@ pi_value <- function(
 }
 ```
 
+Let's try it:
+
+```{r include=FALSE, eval=FALSE}
+pi_value(tau = .9, pd = 10, n = 1e4, le = 1e6) |> 
+  density() |> 
+  with({
+    plot2(x, y, xlab = "proportion of new infections sterilized",
+          ylab = "density of probability", xlim = 0:1)
+    polygon2(x, y)}
+  )
+```
+
 ```{r}
-pi_values <- pi_value()
+tau_vals <- c(.2, .5, .75, .9, 1)
+
+out1 <- tau_vals |>
+  map(~ pi_value(tau = .x, pd = 10, n = 1e4, le = 1e6)) |> 
+  map(density)
+
+plot2(NA, xlim = 0:1, ylim = c(0, max(unlist(map(out1, extract2, "y")))),
+      xlab = "proportion of new infections sterilized", ylab = "density of probability")
+
+walk2(out1, seq_along(tau_vals), ~ with(.x, {
+  polygon2(x, y, col = .y)
+  lines2(x, y, col = .y)
+}))
+
+add_legend_sterilization <- function(...) {
+  legend2("topright", legend = paste0(rev(100 * tau_vals), "%"),
+          col = rev(seq_along(tau_vals)),
+          title = "percentage of recent infections identified")
+}
+
+add_legend_sterilization()
 ```
 
-Once we have the values of $\pi$
+A function that plot a quantile distribution along the $y$-axis:
 
 ```{r}
-N <- 1e5
-I0 <- 2e4
-D0 <- 176
-d0 <- 322
-p_val <- param_val(I0, D0, d0, N)
-dynamics_dfs <- map(pi_values, sid_dyn,
-                    N = N, I0 = I0, D0 = D0,
-                    beta = p_val[["beta"]],
-                    sigma = 0,
-                    delta = p_val[["delta"]],
-                    alpha = p_val[["alpha"]], times = seq2(0, 50))
+quantplot <- function(x, y, probs = c(.025, seq(.05, .95, .05), .975),
+                      col = adjustcolor(4, .07)) {
+  y <- map_dfr(y, quantile, probs)
+  nbcol <- ncol(y)
+  nb <- (nbcol - (nbcol %% 2)) / 2 - 1
+  walk(0:nb, ~ polygon3(x, y[[1 + .x]], y[[nbcol - .x]], col = col))
+}
 ```
 
+Let's generate the distribution of the proportion of new infections sterilized for
+various proposed durations of treatment (about 1'40"):
+
 ```{r}
-plot(NA, xlim = c(0, 50), ylim = c(0, 350))
-walk(dynamics_dfs, ~ with(.x, lines(time, I * p_val[["delta"]], col = adjustcolor(4, .05))))
+pd_vals <- seq(1, 40, le = 512)
+out <- map(tau_vals,
+           ~ mclapply2(pd_vals,
+                       function(x) pi_value(pd = x, tau = .x, n = 1e4, le = 1e6)))
 ```
 
+This gives:
+
 ```{r}
-simulations <- function(
-    tau = .95,
-    pd = 20,
-    times = seq2(0, 50),
-    N = 1e5, I0 = 2e4, D0 = 176, d0 = 322,
-    Xv  = 15,             hv = 4,
-    Xf  = 200,  Yf = .15, hf = 7,
-    Xmp = .875, Ym = 1,   hm = 50,
-    Xe  = 10,   Ye = 1,   he = 1,
-    by = pd / (le - 1), le = 512,
-    n = 1000) {
-  p_val <- param_val(I0, D0, d0, N)
-  pi_value(tau, pd, Xv, hv, Xf, Yf, hf, Xmp, Ym, hm, Xe, Ye, he, by, le, n) |> 
-    map(sid_dyn, N = N, I0 = I0, D0 = D0, beta = p_val[["beta"]], sigma = 0,
-        delta = p_val[["delta"]], alpha = p_val[["alpha"]], times = times)
+plot(NA, xlim = c(0, 45), ylim = 0:1,
+     xlab = "proposed duration of treatment (days)",
+     ylab = "proportion of new infections sterilized")
+
+walk2(out,
+      seq_along(out),
+      ~ quantplot(pd_vals, .x,  probs = seq(.025, .975, le = 75),
+                  col = adjustcolor(.y, .03)))
+
+add_legend_sterilization()
+```
+
+
+## Epidemiological consequences
+
+### A minimalist SID
+
+![Flow diagram of a minimalist SID model.](SIsIfD.png){width=320 #fig-tinySID}
+
+#### Dynamics
+
+$$
+\begin{align} 
+  \frac{dS}{dt} &=  \gamma_s I_s + \gamma_f I_f + \alpha D - (1 - p\pi)\beta D S\\
+  \frac{dI_s}{dt} &=  (1 - p) \beta D S - (\gamma_s + \delta_s) I_s \\
+  \frac{dI_f}{dt} &=  (1-\pi) p \beta D S - (\gamma_f + \delta_f) I_f  \\
+  \frac{dD}{dt} &= \delta_s I_s + \delta_f I_f - \alpha D
+\end{align}
+$$
+
+with $N = S + I_s + I_f + D$. R code:
+
+```{r}
+sid_dyn <- function(N = 1e5, Is0 = 0, If0 = 0, D0,
+                    pi, p, beta, gammas, gammaf, deltas, deltaf, alpha, times) {
+  ode2(c(S = N - Is0 - If0 - D0, Is = Is0, If = If0, D = D0),
+       times,
+       function(time, state, pars) {
+         with(as.list(c(state, pars)), {
+           infection <- beta * D * S
+           dS  <- gammas * Is + gammaf * If + alpha * D - (1 - p * pi) * infection
+           dIs <- (1 - p) * infection - (gammas + deltas) * Is
+           dIf <- (1 - pi) * p * infection - (gammaf + deltaf) * If
+           dD  <- deltas * Is + deltaf * If - alpha * D
+           list(c(dS, dIs, dIf, dD))
+         })
+       },
+       c(pi = pi, p = p, beta = beta, gammas = gammas, gammaf = gammaf,
+         deltas = deltas, deltaf = deltaf, alpha = alpha))
 }
 ```
 
+#### Equilibrium
+
+$$
+\begin{align} 
+  S^*   &= \frac{(\gamma_f + \delta_f)(\gamma_s + \delta_s)}
+                {(1 - p) \gamma_f \delta_s +
+                 (1 - \pi) p \gamma_s \delta_f +
+                 (1 - p \pi) \delta_f \delta_s} \frac{\alpha}{\beta}                  \\
+  I_s^* &=  \frac{(1 - p) \beta}{\gamma_s + \delta_s} D^* S^*                         \\
+  I_f^* &=  \frac{(1 - \pi) p \beta}{\gamma_f + \delta_f} D^* S^*                     \\
+  D^*   &= \frac{N - S^*}
+                {1 + (1 - p \beta S^*) / (\gamma_s + \delta_s) +
+                     ([1 - \pi] p \beta S^*) / (\gamma_f + \delta_f)}
+\end{align}
+$$
+
+R code:
+
 ```{r}
-system.time(
-out <- map(c(3, 5, 10, 15, 20, 25, 30, 35, 40), ~ simulations(pd = .x))
-)
+sid_equ <- function(N = 1e5, pi, p, beta, gammas, gammaf, deltas, deltaf, alpha) {
+  S <- (gammaf + deltaf) * (gammas + deltas) * alpha / 
+    (((1 - p) * gammaf * deltas +
+        (1 - pi) * p * gammas * deltaf +
+        (1 - p * pi) * deltaf * deltas) * beta)
+  D <- (N - S) /
+    (1 + (1 - p * beta * S) / (gammas + deltas) +
+       ((1 - pi) * p * beta * S) / (gammaf + deltaf))
+  betaDS <- beta * D * S
+  c(S  = S,
+    Is = (1 - p) * betaDS / (gammas + detlas),
+    If = (1 - pi) * p * betaDS / (gammaf + deltaf),
+    D  = D)
+}
 ```
 
+A function that plots the dynamics and equilibrium values of the infection prevalence
+and the disease incidence and prevalence:
+
+```{r}
+plot_dyn_equ <- function(
+    N = 1e5,
+    Is0 = 0,
+    If0 = 0,
+    D0 = 350,
+    pi = 0,
+    p = .9,
+    beta = 1e-6,
+    gammas = 0,
+    gammaf = 0,
+    deltas = 0,
+    deltaf = 1e-4,
+    alpha = 1 / 70,
+    times = seq2(0, 25000)) {
+  sims <- sid_dyn(N, Is0, If0, D0,
+                  pi, p, beta, gammas, gammaf, deltas, deltaf, alpha, times)
+
+  equs <- sid_equ(N, pi, p, beta, gammas, gammaf, deltas, deltaf, alpha)
+  
+  plotl2 <- function(...) plotl(..., xlab = "time (year)")
+  
+  opar <- par(mfrow = c(1, 4)) #, cex = 1, plt = c(.25, .95, .25, .9))
+  with(sims, {
+    plotl2(time, Is, col = 2, ylab = "infection prevalence")
+    abline(h = equs["Is"], col = 2)
+    plotl2(time, If, col = 3, ylab = "infection prevalence")
+    abline(h = equs["If"], col = 3)
+    plotl2(time, D, col = 3, ylab = "disease prevalence")
+    abline(h = equs["D"], col = 3)
+    plotl2(time, deltas * Is + deltaf * If, col = 4, ylab = "yearly disease incidence")
+    abline(h = deltas * equs["Is"] + deltaf * equs["If"], col = 4)
+  })
+  par(opar)
+}
+```