In cmhoove14/DDNTD: Density Dependence in Neglected Tropical Diseases

knitr::opts_chunk$set(echo = TRUE)

devtools::load_all()

Basic Schistosomiasis Model

Beginning with a basic schistosomiasis model with $S-E-I$ infection dynamics and logistic population growth among the intermediate host snail population and human infection modeled via the state variable $W$ representing the mean worm burden across the human population, assumed to be negative binomially distributed with clumping parameter $\kappa$. We have four ODEs:

r latexImg("\\frac{dS}{dt}=f_N(S+E)-\\mu_N S-\\Lambda S")

r latexImg("\\frac{dE}{dt}=\\Lambda S-(\\mu_N+\\sigma)E")

r latexImg("\\frac{dI}{dt}=\\sigma E - \\mu_I I")

r latexImg("\\frac{dW}{dt}=\\lambda-(\\mu_W+\\mu_H)W")

Where $N=S+E+I$, $\Lambda$ is the man-to-snail force of infection (FOI), and $\lambda$ is the snail-to-man FOI, more details on these below.

Parameter symbology, values and definitions are shown in table 1.

tab1 <- data.frame(val = as.character(c(10,
                                        50,
                                        60,
                                        40,
                                        10,
                                        3.3,
                                        1.5,
                                        60,
                                        4820,
                                        6.25e5,
                                        5e-3,
                                        1.8e-3)), 
                   units = c("$Sd^{-1}$",
                             "$Nm^{-2}$",
                             "$Nd^{-1}$",
                             "$Ed^{-1}$",
                             "$Id^{-1}$",
                             "$Wy^{-1}$",
                             "$Hm^{-2}$",
                             "$Hy^{-1}$",
                             "$WI^{-1}d^{-1}$",
                             "$ES^{-1}d^{-1}$",
                             "unitless",
                             "unitless"),
                   des = c("Snail fecundity rate",
                           "Snail environmental carrying capacity",
                           "Snail mortality rate",
                           "Pre-patent period",
                           "Excess mortality of infected snails",
                           "Adult parasite mortality rate",
                           "Human host population",
                           "Human host mortality rate",
                           "Man-to-snail transmission parameter",
                           "Snail-to-man transmission parameter",
                           "Adult parasite density dependent fecundity parameter",
                           "Density dependent parasite establishment acquired immunity parameter"))

rownames(tab1) <- c("$f_N$",
                    "$K$",
                    "$\\mu_N$",
                    "$\\sigma$",
                    "$\\mu_I$",
                    "$\\mu_W$",
                    "$H$",
                    "$\\mu_H$",
                    "$\\lambda$",
                    "$\\beta$",
                    "$\\alpha$",
                    "$\\xi$")

knitr::kable(tab1, row.names = TRUE, 
             col.names = c("Value", "Units", "Description"), 
             format = "markdown", escape = FALSE, 
             caption = "Parameter values and descriptions used in the model")

$R_{eff}$ Derivation

Dimensionality reduction: fast snail infection dynamics

We begin with the assumption that the rate of change of the intermediate host infection dynamics is fast compared to the adult parasites and therefore reaches an equilibirum, i.e. r latexImg("\\frac{dS}{dt}=\\frac{dE}{dt}=\\frac{dI}{dt}=0"). We then solve for the equilibirum infected snail population, $I^*$:

Solve for $I^*$

r latexImg("\\sigma E^* - \\mu_I I^*=0")

r latexImg("I^*=\\frac{\\sigma E^*}{\\mu_I}")

Solve for $E^*$

r latexImg("\\Lambda S^*-(\\mu_N+\\sigma)E^*=0")

r latexImg("E^*=\\frac{\\Lambda S^*}{\\mu_N+\\sigma}")

Solve for $S^*$

First, we write the equilibrium total snail population, $N^$, in terms of $S^$ by substituting for $E^$ and $I^$ from $N^=S^+E^+I^$: r latexImg("N^*=S^*+\\frac{\\Lambda S^*}{\\mu_N+\\sigma}+\\frac{\\sigma E^*}{\\mu_I}")

r latexImg("N^*=S^*+\\frac{\\Lambda S^*}{\\mu_N+\\sigma}+\\frac{\\sigma\\Lambda S^*}{\\mu_I(\\mu_N+\\sigma)}")

r latexImg("N^*=S^*\\Big(1+\\frac{\\Lambda }{\\mu_N+\\sigma}+\\frac{\\sigma\\Lambda }{\\mu_I(\\mu_N+\\sigma)}\\Big)")

Then, with r latexImg("0=f_N\\big(1-\\frac{N^*}{K}\\big)\\big(S^*+E^*\\big)-\\mu_N S^*-\\Lambda S^*")

we substitute for $E^$ and $N^$ and divide by $S^*$ to arrive at:

r latexImg("f_N\\Bigg(1-\\frac{S^*\\big(1+\\frac{\\Lambda }{\\mu_N+\\sigma}+\\frac{\\sigma\\Lambda }{\\mu_I(\\mu_N+\\sigma)}\\big)}{K}\\Bigg)\\Bigg(1+\\frac{\\Lambda }{\\mu_N+\\sigma}\\Bigg)-\\mu_N-\\Lambda =0")

Solving for $S^*$:

r latexImg("1-\\frac{S^*\\big(1+\\frac{\\Lambda }{\\mu_N+\\sigma}+\\frac{\\sigma\\Lambda }{\\mu_I(\\mu_N+\\sigma)}\\big)}{K}=\\frac{\\mu_N+\\Lambda }{f_N\\Big(1+\\frac{\\Lambda }{\\mu_N+\\sigma}\\Big)}")

r latexImg("\\frac{S^*\\Big(1+\\frac{\\Lambda }{\\mu_N+\\sigma}+\\frac{\\sigma\\Lambda }{\\mu_I(\\mu_N+\\sigma)}\\Big)}{K}=1-\\frac{\\mu_N+\\Lambda }{f_N\\Big(1+\\frac{\\Lambda }{\\mu_N+\\sigma}\\Big)}")

r latexImg("S^*=\\Bigg(1-\\frac{\\mu_N+\\Lambda }{f_N\\Big(1+\\frac{\\Lambda }{\\mu_N+\\sigma}\\Big)}\\Bigg)\\Bigg(\\frac{K}{\\Big(1+\\frac{\\Lambda }{\\mu_N+\\sigma}+\\frac{\\sigma\\Lambda }{\\mu_I(\\mu_N+\\sigma)}\\Big)}\\Bigg)")

Now want to try and simplify this. To start, we'll assign r latexImg("C_1=\\frac{\\Lambda }{\\mu_N+\\sigma}") to give:

r latexImg("S^*=\\Bigg(1-\\frac{\\mu_N+\\Lambda }{f_N\\big(1+C_1\\big)}\\Bigg)\\Bigg(\\frac{K}{\\big(1+C_1+\\frac{\\sigma C_1}{\\mu_I}\\big)}\\Bigg)")

Then distribute:

r latexImg("S^*=\\frac{K}{\\big(1+C_1+\\frac{\\sigma C_1}{\\mu_I}\\big)}-\\frac{K\\big(\\mu_N+\\Lambda \\big)}{f_N\\big(1+C_1\\big)\\big(1+C_1+\\frac{\\sigma C_1}{\\mu_I}\\big)}")

Then multiply the LHS by r latexImg("\\frac{f_N\\big(1+C_1\\big)}{f_N\\big(1+C_1\\big)}")

r latexImg("S^*=\\frac{K\\big(f_N\\big(1+C_1\\big)\\big)-K\\big(\\mu_N+\\Lambda \\big)}{f_N\\big(1+C_1\\big)\\big(1+C_1+\\frac{\\sigma C_1}{\\mu_I}\\big)}")

r latexImg("S^*=\\frac{K\\big(f_N+f_N C_1-\\mu_N-\\Lambda \\big)}{f_N\\big(1+C_1\\big)\\big(1+C_1+\\frac{\\sigma C_1}{\\mu_I}\\big)}")

This is basically where we left it in the Elimination Feasibility paper, but I think it might be possible to simplify a bit more by first factoring out $f_N$ from the numerator and then canceling:

r latexImg("S^*=\\frac{K\\cancel{f_N}\\big(1+C_1-\\frac{\\mu_N}{f_N}-\\frac{\\Lambda }{f_N}\\big)}{\\cancel{f_N}\\big(1+C_1\\big)\\big(1+C_1+\\frac{\\sigma C_1}{\\mu_I}\\big)}")

r latexImg("S^*=\\frac{K\\big(1+C_1-\\frac{\\mu_N}{f_N}-\\frac{\\Lambda }{f_N}\\big)}{\\big(1+C_1\\big)\\big(1+C_1+\\frac{\\sigma C_1}{\\mu_I}\\big)}")

Then distribute in the denominator:

r latexImg("S^*=\\frac{K\\big(1+C_1-\\frac{\\mu_N}{f_N}-\\frac{\\Lambda }{f_N}\\big)}{\\frac{2\\sigma C_1^2}{\\mu_I}+\\frac{3\\sigma C_1}{\\mu_I}+1}")

And then factor out r latexImg("\\frac{\\sigma C_1}{\\mu_I}") from the denominator:

r latexImg("S^*=\\frac{K\\big(1+C_1-\\frac{\\mu_N}{f_N}-\\frac{\\Lambda }{f_N}\\big)}{\\frac{\\sigma C_1}{\\mu_I}\\big(2C_1+3+\\frac{\\mu_I}{\\sigma C_1}\\big)}")

Now with the rate of change of the mean worm burden:

r latexImg("\\frac{dW}{dt}=\\lambda I^*-(\\mu_W+\\mu_H)W")

And: r latexImg("I^*=\\Big(\\frac{\\sigma}{\\mu_I}\\Big)\\Big(\\frac{\\Lambda }{\\mu_N+\\sigma}\\Big)S^*") and r latexImg("C_1=\\frac{\\Lambda }{\\mu_N+\\sigma}")

We get:

r latexImg("\\frac{dW}{dt}=\\lambda\\Big(\\frac{\\sigma C_1}{\\mu_I}\\Big)S^*-(\\mu_W+\\mu_H)W")

r latexImg("\\frac{dW}{dt}=\\lambda\\Big(\\cancel{\\frac{\\sigma C_1}{\\mu_I}}\\Big)\\Big(\\frac{K\\big(1+C_1-\\frac{\\mu_N}{f_N}-\\frac{\\Lambda }{f_N}\\big)}{\\cancel{\\frac{\\sigma C_1}{\\mu_I}}\\big(2C_1+3+\\frac{\\mu_I}{\\sigma C_1}\\big)}\\Big)-(\\mu_W+\\mu_H)W")

r latexImg("\\frac{dW}{dt}=\\Big(\\frac{K\\lambda\\big(1+C_1-\\frac{\\mu_N}{f_N}-\\frac{\\Lambda }{f_N}\\big)}{\\big(2C_1+3+\\frac{\\mu_I}{\\sigma C_1}\\big)}\\Big)-(\\mu_W+\\mu_H)W")

Now factor out $(\mu_W+\mu_H)W$ to get:

r latexImg("\\frac{dW}{dt}=\\Big((\\mu_W+\\mu_H)W\\Big)\\Big(\\frac{K\\lambda\\big(1+C_1-\\frac{\\mu_N}{f_N}-\\frac{\\Lambda }{f_N}\\big)}{\\big(2C_1+3+\\frac{\\mu_I}{\\sigma C_1}\\big)\\big((\\mu_W+\\mu_H)W\\big)}-1\\Big)")

Given the definition of $R_{eff}$ as the number of adult worms produced by a single adult worm over its lifespan, we interpret $(\mu_W+\mu_H)$ as the mean lifespan and therefore have:

r latexImg("\\frac{dW}{dt}=(\\mu_W+\\mu_H)W(R_{eff}-1)")

and therefore:

r latexImg("R_{eff}=\\frac{K\\lambda\\big(1+C_1-\\frac{\\mu_N}{f_N}-\\frac{\\Lambda }{f_N}\\big)}{\\big(2C_1+3+\\frac{\\mu_I}{\\sigma C_1}\\big)\\big((\\mu_W+\\mu_H)W\\big)}")

This gives us our first key analytic finding, the identification of the breakpoint. If $R_{eff}<1$, r latexImg("\\frac{dW}{dt}<0") which implies that an $R_{eff}<1$ causes a worm population on the way to extinction (e.g. below the breakpoint)

Now plugging back in for $C_1$:

r latexImg("R_{eff}=\\frac{K\\lambda\\big(1+\\frac{\\Lambda }{\\mu_N+\\sigma}-\\frac{\\mu_N}{f_N}-\\frac{\\Lambda }{f_N}\\big)}{\\big(\\frac{2\\Lambda }{\\mu_N+\\sigma}+3+\\frac{\\mu_I(\\mu_N+\\sigma)}{\\sigma\\Lambda }\\big)\\big((\\mu_W+\\mu_H)W\\big)}")

Then factor out r latexImg("\\frac{1}{f_N(\\mu_N+\\sigma)}") from the numerator, factor out r latexImg("\\frac{1}{\\mu_N+\\sigma}") from the denominator, and separate the pesky $(\mu_W+\mu_H)W$

r latexImg("R_{eff}=\\frac{K\\lambda\\big(\\frac{1}{f_N(\\mu_N+\\sigma)}\\big)\\big(f_N(\\mu_N+\\sigma)+f_N\\Lambda -\\mu_N(\\mu_N+\\sigma)-\\Lambda (\\mu_N+\\sigma)\\big)}{\\frac{1}{\\mu_N+\\sigma}\\big(2\\Lambda+3(\\mu_N+\\sigma)+\\frac{\\mu_I(\\mu_N+\\sigma)^2}{\\sigma\\Lambda}\\big)}\\times\\frac{1}{\\Big((\\mu_W+\\mu_H)W\\Big)}")

The r latexImg("\\frac{1}{\\mu_N+\\sigma}") in both the numerator and denominator cancel, the remaing r latexImg("\\frac{1}{f_N}") can be moved to the denominator leaving:

r latexImg("R_{eff}=\\frac{K\\lambda\\big(f_N(\\mu_N+\\sigma)+f_N\\Lambda -\\mu_N(\\mu_N+\\sigma)-\\Lambda (\\mu_N+\\sigma)\\big)}{f_N\\big(2\\Lambda+3(\\mu_N+\\sigma)+\\frac{\\mu_I(\\mu_N+\\sigma)^2}{\\sigma\\Lambda}\\big)}\\times\\frac{1}{\\Big((\\mu_W+\\mu_H)W\\Big)}")

Some more fenagling of the numerator (distribute and then factor out $(\mu_N+\sigma+\Lambda)$) gives:

r latexImg("R_{eff}=\\frac{K\\lambda\\Big(\\mu_N+\\sigma+\\Lambda\\Big)\\Big(f_N-\\mu_N-\\frac{\\Lambda\\sigma}{(\\mu_N+\\sigma+\\Lambda)}\\Big)}{f_N\\big(2\\Lambda+3(\\mu_N+\\sigma)+\\frac{\\mu_I(\\mu_N+\\sigma)^2}{\\sigma\\Lambda}\\big)}\\times\\frac{1}{\\Big((\\mu_W+\\mu_H)W\\Big)}")

So this tells us for $R_{eff}>0$, r latexImg("f_N>\\mu_N+\\frac{\\Lambda\\sigma}{\\mu_N+\\sigma+\\Lambda}"). I think this makes sense and is interpretable as the intermediate host birth rate, $f_N$, has to be greater than the intermediate host baseline mortality rate, $\mu_N$, and excess mortality from infection, r latexImg("\\frac{\\Lambda\\sigma}{\\mu_N+\\sigma+\\Lambda}"), else the intermediate host population would go extinct. Let's keep this in mind but continue on.

Factoring out $f_N$ from the RHS of the numerator we have:

r latexImg("R_{eff}=\\frac{K\\lambda\\cancel{f_N}\\Big(\\mu_N+\\sigma+\\Lambda\\Big)\\Big(1-\\frac{\\mu_N}{f_N}-\\frac{\\Lambda\\sigma}{f_N(\\mu_N+\\sigma+\\Lambda)}\\Big)}{\\cancel{f_N}\\big(2\\Lambda+3(\\mu_N+\\sigma)+\\frac{\\mu_I(\\mu_N+\\sigma)^2}{\\sigma\\Lambda}\\big)}\\times\\frac{1}{\\Big((\\mu_W+\\mu_H)W\\Big)}")

Now if we distribute and factor out r latexImg("\\frac{1}{\\sigma\\Lambda}") in the denominator we get:

r latexImg("R_{eff}=\\frac{K\\lambda\\Big(\\mu_N+\\sigma+\\Lambda\\Big)\\Big(1-\\frac{\\mu_N}{f_N}-\\frac{\\Lambda\\sigma}{f_N(\\mu_N+\\sigma+\\Lambda)}\\Big)}{\\Big(\\frac{1}{\\sigma\\Lambda}\\Big)\\Big(2\\sigma\\Lambda^2+3\\sigma\\Lambda\\mu_N+3\\sigma^2\\Lambda+\\mu_I\\mu_N^2+2\\mu_I\\mu_N\\sigma+\\mu_I\\sigma^2\\Big)}\\times\\frac{1}{\\Big((\\mu_W+\\mu_H)W\\Big)}")

Now move the $\Lambda$ to the numerator and cancel out a $\sigma$ in the denominator:

r latexImg("R_{eff}=\\frac{K\\lambda\\Lambda\\Big(\\mu_N+\\sigma+\\Lambda\\Big)\\Big(1-\\frac{\\mu_N}{f_N}-\\frac{\\Lambda\\sigma}{f_N(\\mu_N+\\sigma+\\Lambda)}\\Big)}{2\\Lambda^2+3\\Lambda\\mu_N+3\\sigma\\Lambda+\\frac{\\mu_I\\mu_N^2}{\\sigma}+2\\mu_I\\mu_N+\\mu_I\\sigma}\\times\\frac{1}{\\Big((\\mu_W+\\mu_H)W\\Big)}")

Then with some final finagling of the denominator we get:

r latexImg("R_{eff}=\\frac{K\\lambda\\Lambda\\Big(\\mu_N+\\sigma+\\Lambda\\Big)\\Big(1-\\frac{\\mu_N}{f_N}-\\frac{\\Lambda\\sigma}{f_N(\\mu_N+\\sigma+\\Lambda)}\\Big)}{\\Big(\\Lambda(2\\Lambda+3\\mu_N+3\\sigma)+\\mu_I(\\frac{\\mu_N^2}{\\sigma}+2\\mu_N+\\sigma)\\Big)\\Big((\\mu_W+\\mu_H)W\\Big)}")

Which we can also write as:

r latexImg("R_{eff}=\\frac{K\\lambda\\Lambda\\Big(\\mu_N+\\sigma+\\Lambda\\Big)}{\\Big(\\Lambda(2\\Lambda+3\\mu_N+3\\sigma)+\\mu_I(\\frac{\\mu_N^2}{\\sigma}+2\\mu_N+\\sigma)\\Big)\\Big((\\mu_W+\\mu_H)W\\Big)}\\times\\Big(1-\\frac{\\mu_N}{f_N}-\\frac{\\Lambda\\sigma}{f_N(\\mu_N+\\sigma+\\Lambda)}\\Big)")

cmhoove14/DDNTD documentation built on Nov. 23, 2019, 7:04 p.m.