Mathematical model of voluntary vaccination against schistosomiasis

Simplification of the ODE

As in Gao et al. (2011), we simplify the ODEs by introducing the following dimensionless variables: S₁ = Λ₁S₁, I₁ = Λ₁I₁, V₁ = Λ₁V₁, S₂ = Λ₂S₂, I₂ = Λ₂I₂, M = M₀M, P = P, $d_1={\mu }_1+{\delta }_1+\eta ,d_2={\mu }_2+\theta ,d_3={\mu }_2+{\delta }_2+\theta ,\gamma =\frac{{\Lambda }_1{\gamma }_1}{M_0},d_4={\mu }_4+\tau ,\delta ={\gamma }_2{\Lambda }_2,{\alpha }_2=\frac{ϵ}{M_0}$. For simplicity, we discard the bold notation and we obtain the following system. (10)${\displaystyle \frac{d{S}_1}{dt}=1-\frac{{\beta }_{1}P{S}_1}{1+{\alpha }_{1}P}-{\mu }_1{S}_1+\eta I_1-v{S}_1+\omega V_1,}$ (11)${\displaystyle \frac{d{I}_1}{dt}=\frac{{\beta }_{1}P{S}_1}{1+{\alpha }_{1}P}-d_1{I}_1,}$ (12)${\displaystyle \frac{d{V}_1}{dt}=v{S}_1-\omega V_1-{\mu }_1{V}_1,}$ (13)${\displaystyle \frac{dM}{dt}=\gamma I_1-{\mu }_{3}M,}$ (14)${\displaystyle \frac{d{S}_2}{dt}=1-\frac{{\beta }_{2}M{S}_2}{1+{\alpha }_2{M}^2}-d_2{S}_2,}$ (15)${\displaystyle \frac{d{I}_2}{dt}=\frac{{\beta }_{2}M{S}_2}{1+{\alpha }_2{M}^2}-d_3{I}_2,}$ (16)${\displaystyle \frac{dP}{dt}=\delta I_2-d_{4}P.}$

Disease-free equilibrium

The equilibria of the dynamics Eqs. (10)–(16) are obtained by setting the time derivatives to 0 and solving the following system of algebraic equations. (17)${\displaystyle 0=1-\frac{{\beta }_{1}P{S}_1}{1+{\alpha }_{1}P}-{\mu }_1{S}_1+\eta I_1-v{S}_1+\omega V_1,}$ (18)${\displaystyle 0=\frac{{\beta }_{1}P{S}_1}{1+{\alpha }_{1}P}-d_1{I}_1,}$ (19)${\displaystyle 0=v{S}_1-\omega V_1-{\mu }_1{V}_1,}$ (20)${\displaystyle 0=\gamma I_1-{\mu }_{3}M,}$ (21)${\displaystyle 0=1-\frac{{\beta }_{2}M{S}_2}{1+{\alpha }_2{M}^2}-d_2{S}_2,}$ (22)${\displaystyle 0=\frac{{\beta }_{2}M{S}_2}{1+{\alpha }_2{M}^2}-d_3{I}_2,}$ (23)${\displaystyle 0=\delta I_2-d_{4}P.}$ There are two equilibria of the dynamics: the disease-free equilibrium and the endemic equilibrium.

Setting I₁ = I₂ = M = P = 0, the system Eqs. (17)–(23) reduces to (24)${\displaystyle 0=1-{\mu }_1{S}_1-v{S}_1+\omega V_1,}$ (25)${\displaystyle 0=1-d_2{S}_2,}$ (26)${\displaystyle 0=v{S}_1-\omega V_1-{\mu }_1{V}_1.}$ By Eq. (25), $S_2^0=\frac{1}{d_2}$. Adding Eq. (24) and Eq. (26) gives 1 = μ₁(S₁ + V₁). Thus, the disease-free equilibrium is given by (27)${\displaystyle S_1^0=\frac{{\mu }_1+\omega }{{\mu }_1\left({\mu }_1+v+\omega \right)},}$ (28)${\displaystyle S_2^0=\frac{1}{d_2},}$ (29)${\displaystyle V_1^0=\frac{v}{{\mu }_1\left({\mu }_1+v+\omega \right)}.}$

Effective reproduction number

The effective reproduction number, $\mathcal{R}$, is found by the next generation matrix method (van den Driessche & Watmough, 2002).

There are four compartments carrying infections, I₁, M, I₂, P and we will keep them in this order. The rate of new infections is given by (30)${\displaystyle \mathcal{F}={\left[\begin{array}{c}\hfill \frac{{\beta }_{1}P{S}_1}{1+{\alpha }_{1}P},0,\frac{{\beta }_{2}M{S}_2}{1+{\alpha }_2{M}^2},0\hfill \end{array}\right]}^T.}$ Differentiating $\mathcal{F}$ at the disease-free equilibrium, we obtain (31)${\displaystyle F=\left[\begin{array}{cccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {\beta }_1{S}_1^0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\beta }_2{S}_2^0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right].}$ The other transmissions in the system are given by (32)${\displaystyle \mathcal{V}={\left[\begin{array}{c}\hfill -d_1{I}_1,\gamma I_1-{\mu }_{3}M,-d_3{I}_2,\delta I_2-d_{4}P\hfill \end{array}\right]}^T.}$ Differentiating $\mathcal{V}$ at the disease-free equilibrium, we obtain (33)${\displaystyle V=\left[\begin{array}{cccc}\hfill -d_1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill \gamma \hfill & \hfill -{\mu }_3\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -d_3\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \delta \hfill & \hfill -d_4\hfill \end{array}\right].}$ Thus, (34)${\displaystyle V^{-1}=\left[\begin{array}{cccc}\hfill -\frac{1}{d_1}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill -\frac{\gamma }{d_1{\mu }_3}\hfill & \hfill -\frac{1}{{\mu }_3}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -\frac{1}{d_3}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -\frac{\delta }{d_3{d}_4}\hfill & \hfill -\frac{1}{d_4}\hfill \end{array}\right],}$ and (35)${\displaystyle F{V}^{-1}=\left[\begin{array}{cccc}\hfill 0\hfill & \hfill 0\hfill & \hfill -\frac{S_1^0{\beta }_1\delta }{d_3{d}_4}\hfill & \hfill -\frac{S_1^0{\beta }_1}{d_4}\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill -\frac{S_2^0{\beta }_2\gamma }{d_1{\mu }_3}\hfill & \hfill -\frac{S_2^0{\beta }_2}{{\mu }_3}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right].}$ The largest eigenvalue of FV⁻¹ is (36)${\displaystyle \mathcal{R}=\rho \left(F{V}^{-1}\right)=\sqrt{\frac{S_1^0{S}_2^0{\beta }_1{\beta }_2\delta \gamma }{d_1{d}_3{d}_4{\mu }_3}}}$ (37)${\displaystyle =\sqrt{\frac{\left({\mu }_1+\omega \right){\beta }_1{\beta }_2\delta \gamma }{d_1{d}_2{d}_3{d}_4{\mu }_1{\mu }_3\left({\mu }_1+\nu +\omega \right)}}.}$

The disease-free equilibrium is locally asymptotically stable if $\mathcal{R}<1$ and the endemic equilibrium is stable if $\mathcal{R}>1$ (van den Driessche & Watmough, 2002).

Endemic equilibrium

Here we find solutions of the system Eqs. (17)–(23) for the endemic equilibrium with the pathogen still present in the environment. By Eq. (19), (40)${\displaystyle V_1=\frac{v{S}_1}{\omega +{\mu }_1}.}$ Adding Eqs. (17)–(19) yields (41)${\displaystyle 1={\mu }_1\left(S_1+V_1\right)-\left(d_1-\eta \right)I_1.}$ Thus, (42)${\displaystyle S_1=\frac{1-\left(d_1-\eta \right)I_1}{{\mu }_1\left(1+\frac{v}{\omega +{\mu }_1}\right)}=\left(1-\left(d_1-\eta \right)I_1\right)S_1^0.}$ By Eq. (20), (43)${\displaystyle M=\frac{\gamma I_1}{{\mu }_3}.}$ By Eq. (21), (44)${\displaystyle 1+{\alpha }_2{M}^2={\beta }_{2}M{S}_2+d_2\left(1+{\alpha }_2{M}^2\right)S_2.}$ Thus, (45)${\displaystyle \frac{S_2}{1+{\alpha }_2{M}^2}=\frac{1}{{\beta }_{2}M+d_2\left(1+{\alpha }_2{M}^2\right)}=\frac{{\mu }_3^2}{{\beta }_2\gamma {\mu }_3{I}_1+d_2\left({\mu }_3^2+{\alpha }_2{\gamma }^2{I}_1^2\right)}.}$ By Eq. (22), (46)${\displaystyle I_2=\frac{{\beta }_{2}M}{d_3}\frac{S_2}{1+{\alpha }_2{M}^2}=\frac{{\beta }_2\gamma {\mu }_3{I}_1}{d_3\left[{\beta }_2\gamma {\mu }_3{I}_1+d_2\left({\mu }_3^2+{\alpha }_2{\gamma }^2{I}_1^2\right)\right]}.}$ By Eq. (23), (47)${\displaystyle P=\frac{\delta I_2}{d_4}=\frac{\delta {\beta }_2\gamma {\mu }_3{I}_1}{d_3{d}_4\left[{\beta }_2\gamma {\mu }_3{I}_1+d_2\left({\mu }_3^2+{\alpha }_2{\gamma }^2{I}_1^2\right)\right]}.}$ Plugging Eqs. (42) and (47) into (18), or, equivalently, into d₁I₁ + α₁Pd₁I₁ = β₁PS₁, and then simplifying, yields the following cubic equation for I₁ (48)${\displaystyle I_1^{\ast }\left(a_1{I}_1^{\ast 2}+a_2{I}_1^{\ast }+a_3\right)=0,}$ where (49)${\displaystyle a_1=\frac{d_1{d}_2{d}_3{d}_4{\alpha }_2{\gamma }^2}{S_1^0},}$ (50)${\displaystyle a_2={\beta }_2\gamma {\mu }_3\left(\delta {\beta }_1\left(d_1-\eta \right)+\frac{d_1{d}_3{d}_4}{S_1^0}+\frac{d_1{\alpha }_1\delta }{S_1^0}\right),}$ (51)${\displaystyle a_3=\frac{d_1{d}_2{d}_3{d}_4{\mu }_3^2}{S_1^0}-\delta {\beta }_1{\beta }_2\gamma {\mu }_3.}$

Note that (52)${\displaystyle a_3={\beta }_1{\beta }_2\delta \gamma {\mu }_3\left(\frac{1}{{\mathcal{R}}^2}-1\right).}$

Thus, a₃ < 0 if and only if $\mathcal{R}>1$. Consequently, Eq. (48) has a unique positive root (53)${\displaystyle I_1^{\ast }=\frac{-a_2+\sqrt{a_2^2-4a_1{a}_3}}{2a_1}}$ if and only if $\mathcal{R}>1$. Once $I_1^{\ast }$ is given by Eq. (53) as a solution of Eq. (48), the other compartments are given by (54)${\displaystyle S_1^{\ast }=\left(1-\left(d_1-\eta \right)I_1^{\ast }\right)S_1^0,}$ (55)${\displaystyle V_1^{\ast }=\frac{v{S}_1^{\ast }}{\omega +{\mu }_1},}$ (56)${\displaystyle M^{\ast }=\frac{\gamma I_1^{\ast }}{{\mu }_3},}$ (57)${\displaystyle S_2^{\ast }=\frac{1+{\alpha }_2{M}^{\ast 2}}{{\beta }_2{M}^{\ast }+d_2\left(1+{\alpha }_2{M}^{\ast 2}\right)},}$ (58)${\displaystyle I_2^{\ast }=\frac{{\beta }_2{M}^{\ast }}{d_3}\frac{S_2^{\ast }}{1+{\alpha }_2{M}^{\ast 2}},}$ (59)${\displaystyle P^{\ast }=\frac{\delta I_2^{\ast }}{d_4}.}$

Finding optimal individual vaccination strategy

To find a Nash equilibrium, we have to solve (60)${\displaystyle {\pi }_{NV}\left(\nu \right)={\pi }_V\left(\nu \right)+c}$ where π_NV and π_V are given by Eqs. (8) and (9), respectively. Rearranging Eq. (60) yields (61)${\displaystyle \frac{\frac{{\beta }_1{P}^{\ast }}{1+{\alpha }_1{P}^{\ast }}}{\frac{{\beta }_1{P}^{\ast }}{1+{\alpha }_1{P}^{\ast }}+{\mu }_1}=\frac{c}{1-\frac{\omega }{\omega +{\mu }_1}}.}$ We solve it for P^∗ to get (62)${\displaystyle P^{\ast }=\frac{c{\mu }_1}{{\beta }_1-\frac{{\beta }_1\omega }{\omega +{\mu }_1}-c{\beta }_1-c{\mu }_1{\alpha }_1}.}$ Since P^∗ is given by Eq. (59), we get (63)${\displaystyle I_2^{\ast }=\frac{d_{4}c{\mu }_1}{\delta \left({\beta }_1-\frac{{\beta }_1\omega }{\omega +{\mu }_1}-c{\beta }_1-c{\mu }_1{\alpha }_1\right)}.}$

From now on, we will use previous calculations to express ν in terms of $I_2^{\ast }$ given by Eq. (63). By Eqs. (56), (57), and (58), $I_1^{\ast }$ can be expressed in terms of I₂ as (64)${\displaystyle I_1^{\ast }=\frac{\left(1-d_3{I}_2^{\ast }\right){\beta }_2\pm \sqrt{\left(\right.{\beta }_2\left(1-d_3{I}_2^{\ast }\right){\left)\right.}^2-4\left(\frac{I_2^{\ast }{d}_2{d}_3{\alpha }_2\gamma }{{\mu }_3}\right)\left(\frac{I_2^{\ast }{d}_2{d}_3{\mu }_3}{\gamma }\right)}}{2\frac{I_2^{\ast }{d}_2{d}_3{\alpha }_2\gamma }{{\mu }_3}}.}$ By Eq. (48), we can express $S_1^0$ in terms of $I_1^{\ast }$ as (65)${\displaystyle S_1^0=\frac{d_1{d}_2{d}_3{d}_4{\alpha }_2{\gamma }^2{I}_1^{\ast 3}+d_1{\beta }_2\gamma {\mu }_3\left(d_3{d}_4+{\alpha }_1\delta \right)I_1^{\ast 2}+d_1{d}_2{d}_3{d}_4{\mu }_3^2{I}_1^{\ast }}{{\beta }_1{\beta }_2\delta \gamma {\mu }_3\left(I_1^{\ast }-\left(d_1-\eta \right)I_1^{\ast 2}\right)},}$ and, by Eq. (27), we can express ν in terms of $S_1^0$ as (66)${\displaystyle v=\frac{{\mu }_1+\omega -S_1^0{\mu }_1^2-S_1^0{\mu }_1\omega }{S_1^0{\mu }_1}.}$ Hence, the Nash equilibrium is given by Eq. (66), where $S_1^0$ is given in Eq. (65), $I_1^{\ast }$ is given by Eq. (64), and $I_2^{\ast }$ is given by Eq. (58).