Mathematical model of voluntary vaccination against schistosomiasis
 Published
 Accepted
 Received
 Academic Editor
 Dominic Thorrington
 Subject Areas
 Parasitology, Epidemiology, Infectious Diseases, Statistics
 Keywords
 Schistosoma mansoni, Vaccination game, Nash equilibrium, Neglected Tropical Diseases
 Copyright
 © 2024 Lopez et al.
 Licence
 This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, reproduction and adaptation in any medium and for any purpose provided that it is properly attributed. For attribution, the original author(s), title, publication source (PeerJ) and either DOI or URL of the article must be cited.
 Cite this article
 2024. Mathematical model of voluntary vaccination against schistosomiasis. PeerJ 12:e16869 https://doi.org/10.7717/peerj.16869
Abstract
Human schistosomiasis is a chronic and debilitating neglected tropical disease caused by parasitic worms of the genus Schistosoma. It is endemic in many countries in subSaharan Africa. Although there is currently no vaccine available, vaccines are in development. In this paper, we extend a simple compartmental model of schistosomiasis transmission by incorporating the vaccination option. Unlike previous models of schistosomiasis transmission that focus on control and treatment at the population level, our model focuses on incorporating human behavior and voluntary individual vaccination. We identify vaccination rates needed to achieve herd immunity as well as optimal voluntary vaccination rates. We demonstrate that the prevalence remains too high (higher than 1%) unless the vaccination costs are sufficiently low. Thus, we can conclude that voluntary vaccination (with or without mass drug administration) may not be sufficient to eliminate schistosomiasis as a public health concern. The cost of the vaccine (relative to the cost of schistosomiasis infection) is the most important factor determining whether voluntary vaccination can yield elimination of schistosomiasis. When the cost is low, the optimal voluntary vaccination rate is high enough that the prevalence of schistosomiasis declines under 1%. Once the vaccine becomes available for public use, it will be crucial to ensure that the individuals have as cheap an access to the vaccine as possible.
Introduction
Human schistosomiasis is a chronic and debilitating neglected tropical disease caused by parasitic flatworms of the genus Schistosoma (Ross et al., 2002). It is endemic in many countries in Africa, South America, and Asia (Madsen & Stauffe, 2022). Worldwide there are an estimated 800 million people at risk of infection (Steinmann et al., 2006); over 230 million people are infected with about 201.5 million living in Africa (Verjee, 2019).
Schistosoma genus consists of 23 species (Littlewood & Webster, 2017); we will focus on S. mansoni which is endemic throughout subSaharan Africa. The life cycle of Schistosoma mansoni is described, for example in McManus et al. (2018). The cycle involves an intermediate freshwater snail host of Biomphalaria species (Habib et al., 2021) and the definitive human host. Eggs are excreted in the human faeces and they hatch upon contact with water. After hatching, the eggs release freeswimming ciliated larvae, miracidia which seek and penetrate snail hosts. Within the snails, the parasites develop into sporocysts which reproduce asexually to produce numerous larvae, called cercariae. The larvae emerge from snails in response to sunlight, and swim seeking human hosts. Once cercariae penetrate the skin of a human host their tails drop off and the larvae transform into schistosomula. They enter blood vessels and migrate to the liver, where they mature into adults. From the liver, the male and female worms migrate in pairs to the bowel. Females produce eggs which are excreted in faeces and the cycle continues.
Schistosomiasis control efforts include the following strategies:

disease treatment largescale mass drug administration (MDA) of praziquantel (PZQ) (Doenhoff et al., 2009),

health education,

snail intermediate host control, and

water, sanitation and hygiene (WASH) programs (Tchuenté et al., 2017).
Successes in Japan, China, Egypt and in some subSaharan African countries such as Cameroon, Angola, Burkina Faso, Central African Republic, Chad, Congo, Mali, Senegal and Uganda demonstrate that control with progression towards elimination is possible (Rollinson et al., 2013). MDA by PZQ is a costeffective ‘preventive chemotherapy’ and it is currently the strategy of choice and endorsed by WHO (Tchuenté et al., 2017; WHO, 2021). However, this strategy is unsustainable in the long term and interruptions in these MDA programs can lead to rebounds of egg count (Ross et al., 2017). Vaccines are being developed, but none are available yet (Molehin, McManus & You, 2022; Molehin, 2020; Molehin et al., 2016). The vaccine development faces many challenges, including the complexity of the schistosome life cycle, the parasite’s ability to evade the immune system and the lack of adequate animal models for test trials (Molehin, McManus & You, 2022). Furthermore, there is a limited economic incentive to advance novel vaccine platforms as the disease affects the poorest regions of the world (Molehin, McManus & You, 2022).
Mathematical modeling plays a crucial and integral part of disease control and elimination (Anderson & May, 1992; Behrend et al., 2020). Many models exist for schistosomiasis transmission and control, including Woolhouse, Hasibeder & Chandiwana (1996), Spear et al. (2002), Chiyaka & Garira (2009), Zhou et al. (2013), Mbah et al. (2014), Stylianou et al. (2017), Lo et al. (2018), Gurarie et al. (2018), Kadaleka, Abelman & Tchuenche (2021), Kadaleka et al. (2021), Kadaleka, Abelman & Tchuenche (2022), Madubueze et al. (2022) and Ronoh et al. (2021). In Collyer et al. (2019) and Kura et al. (2020), the authors modeled the impact of schistosomiasis vaccine. They found that in high transmission settings, MDA alone is unable to achieve the WHO goals of morbidity control and elimination as a public health problem. However, vaccination is able to achieve both goals in combination with MDA. Other models focus on snail intermediate hosts (Woolhouse, 1991; Woolhouse & Chandiwana, 1990; Feng, Li & Milner, 2002; Allen & Victory Jr, 2003; Zhao & Milner, 2008; Mangal, Paterson & Fenton, 2008; Anderson, Loker & Wearing, 2021). In French et al. (2010), the authors fitted a model to data from a largescale administration of PZQ in Uganda.
The aim of this paper is to focus on incorporating human behavior and voluntary individual vaccination against schistosomiasis. We want to determine whether voluntary vaccination alone could eliminate schistosomiasis as a public health concern, i.e., decrease the prevalence of high intensity infections under 1% (WHO, 2022). We extend a compartmental model presented in Gao et al. (2011) which investigated the effect of MDA on schistosomiasis transmission. Inspired by Stylianou et al. (2017), Kura et al. (2019), we assume the vaccination is already available. We focus on what happens when MDA is no longer in place; similarly to modeling the postMDA development in other NTDs such as trachoma (Barazanji et al., 2023), lymphatic filariasis (Rychtář & Taylor, 2022), or yaws (Kimball et al., 2022). Even if the vaccination is incorporated into existing pediatric vaccine programs and made mandatory by the government, it does not automatically mean that the population would adhere to the mandates. Vaccine hesitancy and avoidance is a real concern in the US (Tolsma, 2015), Europe (Reczulska, Tomaszewska & Raciborski, 2022) as well as Africa (Anjorin et al., 2021). For example, Central Africa has a significantly lower COVID19 vaccine acceptance rate (less than 35%) than Southern Africa (about 75%) (Anjorin et al., 2021). There is a conflict between individual freedom and interests and the public health benefits (Paplicki et al., 2018). The vaccination, if adopted by enough people in the population, produces herd immunity and decreases the disease prevalence. This benefit can be enjoyed even by those not vaccinated (Serpell & Green, 2006). Thus, vaccination programs are prone to freeriding (Ibuka et al., 2014) because individuals maximize their selfinterests (such as avoiding the costs associated with vaccination), rather than the interests of the entire group (Maskin, 1999).
We apply the game theory framework popularized in Bauch & Earn (2004). The framework has been applied to many diseases; see Wang et al. (2016), Verelst, Willem & Beutels (2016) and Chang et al. (2020) for recent reviews. As argued in Wang et al. (2016), epidemics models incorporating human behavior provide more insight and better predictions. Thus, the game theory models have been applied to study the prevention and elimination of many NTDs, mpox (formerly monkeypox, Bankuru et al. (2020), Augsburger et al. (2022), Augsburger et al. (2023), chikungunya (Klein et al., 2020), typhoid fever (AcostaAlonzo et al., 2020), Chagas disease (Han et al., 2020), visceral leishmaniasis (Fortunato et al., 2021), lymphatic filariasis (Rychtář & Taylor, 2022), rabies (Campo et al., 2022), yellow fever (Caasi et al., 2022), or zika (Angina et al., 2022).
In the ideal case, the interests of the individual—to minimize one’s costs, or to maximize one’s benefits—align with the interest of the entire population—to reduce the prevalence of the disease below a certain threshold such as 1% for children age 5–14 (WHO, 2022). If this is the case, by behaving optimally (in their own sense), the individuals will behave optimally from the public health perspective. Thus, the individuals will more likely adhere to the mandatory vaccination policy and contribute to disease elimination as the public health concern. However, because interests can differ, a behavior that is optimal from the perspective of an individual may not be optimal from the perspective of the group and vice versa. To avoid confusion, in the rest of the paper, when we say “optimal”, we will mean optimal from the individual perspective, unless specified otherwise.
Material and methods
We introduce a mathematical model for voluntary vaccination against schistosomiasis. First, we incorporate a possible vaccination into a compartmental model of schistosomiasis transmission developed by Gao et al. (2011). Then, we add the game theory component that will allow us to investigate individuals’ optimal vaccination decisions.
Compartmental model
The human population is divided into susceptible (S_{1}), infectious (I_{1}) and vaccinated (V_{1}). The snail population is divided into susceptible (S_{2}) and infected (I_{2}). The schistosomiasis pathogen is divided into the snailpenetrating stage miracidia (M), and the humanpenetrating stage, cercariae (P).
Human individuals are born susceptible to schistosomiasis at a rate Λ_{1}. Susceptible individuals become infected through contact with freeliving cercariae present in contaminated fresh water. Because of saturating and crowding effect, we use a Holling type II incidence rate $\frac{{\beta}_{1}P}{1+{\alpha}_{1}P}$ (Holling, 1959; Real, 1977), where β_{1} is the rate of transmission in small concentrations of P and α_{1} is a scaling factor.
As in Gao et al. (2011), we assume that the infected humans are treated at rate η, returning to the susceptible population; without treatment the individuals stay infected.
Susceptible individuals are vaccinated at a rate ν. Vaccinated individuals are assumed immune against the disease. They lose their vaccineinduced immunity at a rate ω and become susceptible again. Infected humans may get vaccinated as well. From a practical standpoint, individuals with low intensity of infection will likely consider themselves susceptible and would vaccinate. Nevertheless, we assume that the vaccine does not work in these instances and the individuals stay infected. Infected humans release parasite eggs giving rise to the population of miraccidia M at rate γ_{1}; we ignore the egg hatching period.
Susceptible snails are born at rate Λ_{2}. They become infected at a rate $\frac{{\beta}_{2}M{S}_{2}}{{M}_{0}+\u03f5{M}^{2}}$ which is a Holling Type III incidence rate (Holling, 1959; Real, 1977), where β_{2} is the rate transmission in small concentrations of M and M_{0} and ϵ are scaling factors. Infected snails give rise to the population of cercariae P at a rate, γ_{2}.
For simplicity, we assume that the risk of contracting schistosomiasis after the age ${\mu}_{1}^{1}$ is negligible. Thus, all humans are removed from the population at risk at a rate μ_{1} as they age. The infected cases also suffer from the diseaserelated death rate δ_{1}; so they are removed from the population at a total rate μ_{1} + δ_{1}. The susceptible snails die at a rate μ_{2} + θ, where μ_{2} is the natural death rate and θ is the elimination rate of snails. Infected snails die at a rate μ_{3} + δ_{2} + θ, where δ_{2} is the diseaserelated death rate of snails. miracidia (M) die at a rate μ_{3}. The death rate of cercariae population P is μ_{4} + τ where μ_{4} is the natural death rate and τ is the elimination rate. We ignore the negligible removal rates of miracidia and cercariae due to human and snail infections.
The transmission dynamic is illustrated in Fig. 1. The notation is summarized in Table 1.
Symbol  Meaning  Value  Range  Source 

Λ_{1}  Birth rate (humans)  0.031  [0.02, 0.04]  World Bank (2022) 
${\mu}_{1}^{1}$  Max age of people at risk  20  [15, 25]  Jordan (1972) 
μ_{2}  Natural death rate (snails)  1.85  [1.5, 2.4]  Appleton (1977) 
μ_{3}  Natural death rate (miracidia)  1460  [1100, 1750]  Maldonado, AcostaMatienzo et al. (1948) 
μ_{4}  Natural death rate (cercariae)  830  [500, 1100]  Whitfield et al. (2003) 
γ_{1}  Miracidia production rate  1.1 × 10^{5}  [10^{5}, 2 × 10^{5}]  Alwan & LoVerde (2021) 
γ_{2}  Cercariae production rate  1.55 × 10^{5}  [0.9 × 10^{5}, 2.2 × 10^{5}]  Gabrielli & Garba Djirmay (2022) 
δ_{1}  Disease related mortality rate (humans)  10^{−4}  [0, 10^{−2}]  WHO (2021) 
η  MDA treatment rate of humans  0  –  Assumed 
τ  Elimination rate of cercariae  0  –  Assumed 
θ  Elimination rate of snails  0  –  Assumed 
ν  Vaccination rate  variable  [0, 0.1]  Assumed 
ω  Vaccine waning rate  1/6.5  [1/8, 1/5]  Zhang et al. (2014) 
δ_{2}  Disease related mortality rate (snails)  0.25  [0, 0.5]  Fitted 
β_{1}  Human infection rate by cercariae  0.0013  [0.001, 0.0015]  Fitted 
α_{1}  Scaling factor for human infection rate  0.0315  [0.01, 0.05]  Fitted 
β_{2}  Snails infection rate by miracidia  12.71  [10, 15]  Fitted 
M_{0}  Scaling factor for snail infection rate  3500  [3000, 5000]  Fitted 
ɛ  Scaling factor for snail infection rate  1.689  [1, 2]  Fitted 
Λ_{2}  Birth rate (snails)  10  [5, 15]  Fitted 
c  Cost of vaccine relative to cost of schistosomiasis  0.02  [0, 0.1]  Assumed 
d_{1}  Rate out of I_{1}  μ_{1} + δ_{1} + η  
d_{2}  Rate out of S_{2}  μ_{2} + θ  
d_{3}  Rate out of I_{2}  μ_{2} + δ_{2} + θ  
d_{4}  Rate out of P  μ_{4} + θ  
γ  Auxiliary variable  $\frac{{\Lambda}_{1}{\gamma}_{1}}{{M}_{0}}$  
δ  Auxiliary variable  γ_{2}Λ_{2}  
α_{2}  Auxiliary variable  $\frac{\varepsilon}{{M}_{0}}$ 
The model yields the following differential equations. (1)$\frac{d{S}_{1}}{dt}={\Lambda}_{1}\frac{{\beta}_{1}P{S}_{1}}{1+{\alpha}_{1}P}{\mu}_{1}{S}_{1}+\eta {I}_{1}v{S}_{1}+\omega {V}_{1},$ (2)$\frac{d{I}_{1}}{dt}=\frac{{\beta}_{1}P{S}_{1}}{1+{\alpha}_{1}P}\left({\mu}_{1}+{\delta}_{1}+\eta \right){I}_{1},$ (3)$\frac{dM}{dt}={\gamma}_{1}{I}_{1}{\mu}_{3}M,$ (4)$\frac{d{S}_{2}}{dt}={\Lambda}_{2}\left({\mu}_{2}+\theta \right){S}_{2},$ (5)$\frac{d{I}_{2}}{dt}=\frac{{\beta}_{2}M{S}_{2}}{{M}_{0}+\u03f5{M}^{2}}\left({\mu}_{2}+{\delta}_{2}+\theta \right){I}_{2},$ (6)$\frac{dP}{dt}={\gamma}_{2}{I}_{2}\left({\mu}_{4}+\tau \right)P,$ (7)$\frac{d{V}_{1}}{dt}=v{S}_{1}\omega {V}_{1}{\mu}_{1}{V}_{1}.\phantom{\rule{0.25em}{0ex}}$
Game theory component
We add a game theory component to study individual vaccination based on the framework introduced in Bauch & Earn (2004) for childhood diseases and used in many other settings, including for the recent COVID19 (Choi & Shim, 2021) or mpox outbreaks (Augsburger et al., 2022).
The game is played by susceptible individuals. As in Bauch & Earn (2004), we assume the players are rational, act in their own best selfinterests, and have complete information about the schistosomiasis epidemic. The individuals decide whether to vaccinate or stay unvaccinated. The payoff is a function of the individual’s vaccination status and the vaccination status of the rest of the population. The payoff incorporates the cost of the vaccination (relative to the cost of the infection which can be assumed as 1 unit), c, the risk of getting infected, π_{NV} if not vaccinated and π_{V} if vaccinated. To evaluate the probability of getting infected, we assume that the epidemics reached a steady state with P^{∗} cercariae given later by Eq. (59); P^{∗} depends on ν, the vaccination rate in the population, but not on the decision of the focal individual. The probability that an unvaccinated individual becomes infected (8)$\pi NV\left(\nu \right)=\frac{\frac{{\beta}_{1}P\ast}{1+{\alpha}_{1}{P}^{\ast}}}{\frac{{\beta}_{1}P\ast}{1+{\alpha}_{1}{P}^{\ast}}+{\mu}_{1}},$ where $\frac{{\beta}_{1}P\ast}{1+{\alpha}_{1}{P}^{\ast}}+{\mu}_{1}$ is a rate at which susceptible individuals with no intention to vaccinate leave the Susceptible compartment and $\frac{{\beta}_{1}P\ast}{1+{\alpha}_{1}{P}^{\ast}}$ is the rate at which they acquire the infection. Similarly, the probability of infection of a vaccinated individual is given by (9)$\pi V\left(\nu \right)=\frac{\omega}{\omega +{\mu}_{1}}{\pi}_{NV}\left(\nu \right),$ where the first term represents the probability of the vaccine waning during the individual’s lifetime.
The solution of the game, called the Nash equilibrium, is the populationlevel vaccination rate—denoted ν_{NE}—at which no individual can increase their own benefits by deviating from the population strategy. It follows that either (1) ν_{NE} = 0 when π_{NV}(0) ≤ π_{V}(0) + c, i.e., when the expected cost of not vaccinating is smaller than the expected cost of vaccinating in a population where nobody else vaccinates, or (2) ν_{NE} solves π_{NV}(ν) = π_{V}(ν) + c, i.e., when the expected payoffs of notvaccinating or vaccinating are equal. Here, π_{NV}(0) is evaluated from Eq. (8) by substituting ν = 0 for the vaccination rate and solving for the equilibria of the system Eqs. (10)–(16) as done in the later sections; the probability π_{V} is evaluated analogously by Eq. (9); c is the cost of vaccine relative to cost of schistosomiasis infection, i.e., C_{vaccine}/C_{Schistosomiasis}. Thus, while not explicitly written, P^{∗} in Eq. (8) is a function of ν and thus π_{NV}(ν) = π_{V}(ν) + c is an equation involving ν even after we substitute for π_{NV} and π_{V} from Eqs. (8) and (9).
Analysis
Simplification of the ODE
As in Gao et al. (2011), we simplify the ODEs by introducing the following dimensionless variables: S_{1} = Λ_{1}S_{1}, I_{1} = Λ_{1}I_{1}, V_{1} = Λ_{1}V_{1}, S_{2} = Λ_{2}S_{2}, I_{2} = Λ_{2}I_{2}, M = M_{0}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{\u03f5}{{M}_{0}}$. For simplicity, we discard the bold notation and we obtain the following system. (10)$\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)$\frac{d{I}_{1}}{dt}=\frac{{\beta}_{1}P{S}_{1}}{1+{\alpha}_{1}P}{d}_{1}{I}_{1},$ (12)$\frac{d{V}_{1}}{dt}=v{S}_{1}\omega {V}_{1}{\mu}_{1}{V}_{1},$ (13)$\frac{dM}{dt}=\gamma {I}_{1}{\mu}_{3}M,$ (14)$\frac{d{S}_{2}}{dt}=1\frac{{\beta}_{2}M{S}_{2}}{1+{\alpha}_{2}{M}^{2}}{d}_{2}{S}_{2},$ (15)$\frac{d{I}_{2}}{dt}=\frac{{\beta}_{2}M{S}_{2}}{1+{\alpha}_{2}{M}^{2}}{d}_{3}{I}_{2},$ (16)$\frac{dP}{dt}=\delta {I}_{2}{d}_{4}P.$
Diseasefree 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)$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)$0=\frac{{\beta}_{1}P{S}_{1}}{1+{\alpha}_{1}P}{d}_{1}{I}_{1},$ (19)$0=v{S}_{1}\omega {V}_{1}{\mu}_{1}{V}_{1},$ (20)$0=\gamma {I}_{1}{\mu}_{3}M,$ (21)$0=1\frac{{\beta}_{2}M{S}_{2}}{1+{\alpha}_{2}{M}^{2}}{d}_{2}{S}_{2},$ (22)$0=\frac{{\beta}_{2}M{S}_{2}}{1+{\alpha}_{2}{M}^{2}}{d}_{3}{I}_{2},$ (23)$0=\delta {I}_{2}{d}_{4}P.$ There are two equilibria of the dynamics: the diseasefree equilibrium and the endemic equilibrium.
Setting I_{1} = I_{2} = M = P = 0, the system Eqs. (17)–(23) reduces to (24)$0=1{\mu}_{1}{S}_{1}v{S}_{1}+\omega {V}_{1},$ (25)$0=1{d}_{2}{S}_{2},$ (26)$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 = μ_{1}(S_{1} + V_{1}). Thus, the diseasefree equilibrium is given by (27)${S}_{1}^{0}=\frac{{\mu}_{1}+\omega}{{\mu}_{1}\left({\mu}_{1}+v+\omega \right)},$ (28)${S}_{2}^{0}=\frac{1}{{d}_{2}},$ (29)${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_{1}, M, I_{2}, P and we will keep them in this order. The rate of new infections is given by (30)$\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 diseasefree equilibrium, we obtain (31)$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)$\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 diseasefree equilibrium, we obtain (33)$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)${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)$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^{−1} is (36)$\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)$=\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 diseasefree equilibrium is locally asymptotically stable if $\mathcal{R}<1$ and the endemic equilibrium is stable if $\mathcal{R}>1$ (van den Driessche & Watmough, 2002).
Critical vaccination rates
The value of v needed to eliminate schistosomiasis can be found by solving (38)$\mathcal{R}=\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)}}<1$ for ν. It follows that whenever ν > ν_{HI}, where (39)${\nu}_{HI}=max\left\{0,\left({\mu}_{1}+\omega \right)\left(\frac{{\beta}_{1}{\beta}_{2}\delta \gamma}{{d}_{1}{d}_{2}{d}_{3}{d}_{4}{\mu}_{1}{\mu}_{3}}1\right)\right\},$ then $\mathcal{R}<1$ and the disease can be eliminated.
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)${V}_{1}=\frac{v{S}_{1}}{\omega +{\mu}_{1}}.$ Adding Eqs. (17)–(19) yields (41)$1={\mu}_{1}\left({S}_{1}+{V}_{1}\right)\left({d}_{1}\eta \right){I}_{1}.$ Thus, (42)${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)$M=\frac{\gamma {I}_{1}}{{\mu}_{3}}.$ By Eq. (21), (44)$1+{\alpha}_{2}{M}^{2}={\beta}_{2}M{S}_{2}+{d}_{2}\left(1+{\alpha}_{2}{M}^{2}\right){S}_{2}.$ Thus, (45)$\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)${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)$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_{1}I_{1} + α_{1}Pd_{1}I_{1} = β_{1}PS_{1}, and then simplifying, yields the following cubic equation for I_{1} (48)${I}_{1}^{\ast}\left({a}_{1}{I}_{1}^{\ast 2}+{a}_{2}{I}_{1}^{\ast}+{a}_{3}\right)=0,$ where (49)${a}_{1}=\frac{{d}_{1}{d}_{2}{d}_{3}{d}_{4}{\alpha}_{2}{\gamma}^{2}}{{S}_{1}^{0}},$ (50)${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)${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)${a}_{3}={\beta}_{1}{\beta}_{2}\delta \gamma {\mu}_{3}\left(\frac{1}{{\mathcal{R}}^{2}}1\right).$
Thus, a_{3} < 0 if and only if $\mathcal{R}>1$. Consequently, Eq. (48) has a unique positive root (53)$I}_{1}^{\ast}=\frac{{a}_{2}+\sqrt{{a}_{2}^{2}4{a}_{1}{a}_{3}}}{2{a}_{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)${S}_{1}^{\ast}=\left(1\left({d}_{1}\eta \right){I}_{1}^{\ast}\right){S}_{1}^{0},$ (55)${V}_{1}^{\ast}=\frac{v{S}_{1}^{\ast}}{\omega +{\mu}_{1}},$ (56)${M}^{\ast}=\frac{\gamma {I}_{1}^{\ast}}{{\mu}_{3}},$ (57)${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)${I}_{2}^{\ast}=\frac{{\beta}_{2}{M}^{\ast}}{{d}_{3}}\frac{{S}_{2}^{\ast}}{1+{\alpha}_{2}{M}^{\ast 2}},$ (59)${P}^{\ast}=\frac{\delta {I}_{2}^{\ast}}{{d}_{4}}.$
Finding optimal individual vaccination strategy
To find a Nash equilibrium, we have to solve (60)${\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)$\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)${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)${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_{2} as (64)${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)${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)$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).
Model calibration
We focus on transmission of S. mansoni and we locate as many parameters specific to this species as possible. However, since S. haematobium is also endemic in subSaharan Africa, some parameter estimates are based on that species or simply schistosoma species in general; we specifically say so if it is the case. We perform sensitivity and uncertainty analysis to account for possible discrepancies in parameter values.
For birth rate, we will use a country in subSaharan Africa, like Zimbabwe where schistosomiasis in general is endemic (Midzi et al., 2014). In Zimbabwe, the birth rate is 31 births per 1,000 people per year (World Bank, 2022), i.e., Λ_{1} = 0.031.
The egg output of cases infected by S. haematobium (Bradley & McCullough, 1973) as well as the length of water contact (Jordan, 1972) varies by age and there is a sharp drop off after the age 20 for both measures (Kura et al., 2021). We will thus assume the same is true for S. mansoni and consider the aging rate μ_{1} = 1/20.
We will consider snails of the Planorbidae family, especially Biomphalaria species, as they are one are a common intermediate host of schistosomiasis (Gabrielli & Garba Djirmay, 2022). Their life span ranges between 5 to 8 months (Appleton, 1977) and we use the average death rate μ_{2} = 12/6.5 ≈ 1.85 per year. The longevity of S. mansoni miracidia is relatively small, about 56 h and no more than 9 h (Maldonado, AcostaMatienzo et al., 1948). We will thus use μ_{3} = 365/(6/24) = 1, 460. Similarly, S. mansoni cercariae live on average about 10.5 h with a range from 817 h (Whitfield et al., 2003) and so we set (μ_{4} = 365 × 24/10.5 ≈ 830). We note that the cercariae may survive up to 72 h (Nelwan, 2019).
S. mansoni females release about 300 eggs per day (Alwan & LoVerde, 2021; Mooee, Sandgeound et al., 1956); we will thus use γ_{1} = 300∗365 ≈ 1.1 × 10^{5}.
The number of S. mansoni cercariae produced daily is 250–600 (Gabrielli & Garba Djirmay, 2022). We will thus use γ_{2} = 425∗365 ≈ 1.55 × 10^{5}.
We estimate the disease related mortality as δ_{1} = 1/10^{4} based on 2016 global schistosomiasis data of 24,000 deaths and 240 million infections (Gabrielli & Garba Djirmay, 2022; WHO, 2021). This is in general agreement with Kheir et al. (1999) who estimated the annual mortality between 50/10^{5} and 1/1000 (or higher for specific kinds of infections.
There is currently no vaccine (Molehin, McManus & You, 2022) for humans. Nevertheless, based on phase 1 clinical trials in baboons, the longevity of one of the tested vaccines is 5–8 years (Zhang et al., 2014). We thus set the vaccine waning rate to be ω = 1/6.5. The vaccine reduces the parasitic female load by about 90%, but for simplicity we will assume a complete protection.
For the purpose of the model, we will assume η = 0 because PZQ helps to control morbidity by killing adult schistosomes but it is ineffective against juvenile worms (McManus et al., 2018; Hagan et al., 2004). We also assume c = 0.02 with the range [0, 0.05], the cost of the vaccine is about 1/50 of the cost of contracting schistosomiasis (and somewhere between 0 and 1/20 of the cost of the disease).
To find the values of other parameters, we set the controls to 0, i.e., set ν = 0, θ = 0, τ = 0, η = 0, and fitted the model predictions to observed data of (a) the reproduction number, ${\mathcal{R}}_{0}\approx \sqrt[4]{4.31}$ based on Woolhouse, Hasibeder & Chandiwana (1996), (b) the proportion of infected individuals, I_{1}/(I_{1} + S_{1}) ≈ 0.227 based on Midzi et al. (2014), and (c) the proportion of infected snails I_{2}/(I_{2} + S_{2}) ≈ 0.018 based on OdongoAginya et al. (2008). We used MATLAB’s builtin optimization procedure fmincon which is a nonlinear programming solver that returns a minimizer of a given function subject to various constraints. We note that Woolhouse, Hasibeder & Chandiwana (1996) estimated the average number of female worms that reach reproductive age produced by a typical female worm over the course of its life by 4.31. Our model has four stages of a parasite (miracidia, parasites in snails, cercaria, and parasites in humans). Thus, the ${\mathcal{R}}_{0}$ derived by the nextgeneration matrix method should satisfy ${\mathcal{R}}_{0}^{4}\approx 4.31$ so that a typical miracidia produces on average 4.31 miracidia during the full cycle.
Results
For the parameter values specified in Table 1, the vaccination rate leading to elimination of schistosomiasis is given by ν_{HI} ≈ 0.23. This is illustrated in Fig. 2. It means that the entire population needs to be vaccinated in about 4.5 years, i.e., slightly more frequently than the minimal vaccination waning rate.
The optimal voluntary vaccination rate is ν_{NE} ≈ 0.16. At this rate, the entire population would be vaccinated in about 6.25 years, just under the assumed waning rate.
The prevalence when individuals use the optimal voluntary vaccination is about 4.7%. We can thus see that after the termination of MDA and other control measures, schistosomiasis would not be eliminated as a public health concern (currently defined as <1% proportion of heavy intensity schistosomiasis infections in school age children WHO (2021)) by optimal voluntary vaccination alone.
Figure 3A shows the dependence of the optimal voluntary vaccination rate ν_{NE} on the relative cost of vaccination, c. Once the cost of vaccination grows above about 0.053, ν_{NE} = 0. It means that if the cost of vaccination is higher than about 1/20 of the cost of schistosomiasis infection, it is not beneficial for the individuals to vaccinate. On the other hand, when the cost of vaccination is very low, then ν_{NE} ≈ ν_{HI}, meaning that schistosomiasis would be very close to elimination.
Similarly, Fig. 3B shows the dependence of the effective reproduction number on c. In agreement with Fig. 3A, when c ≈ 0, $\mathcal{R}\approx 1$ and when c > 0.053, $\mathcal{R}\approx 1.45$. Note that as long as c > 0, $\mathcal{R}>1$, i.e., the optimal voluntary vaccination will never completely eliminate schistosomiasis on its own. Finally, Fig. 3C shows the dependence of the prevalence on c. It follows that as long as c < 0.005, the prevalence is less than 1%, i.e., schistosomiasis would be considered eliminated as a public health concern.
Figure 4 shows how the outcomes depend on the MDA rate, η. The optimal voluntary vaccination rate, ν_{NE} is positive for η < 0.04 while the vaccination rate needed for her immunity, ν_{HI} is positive for η < 0.057. Moreover, the prevalence of schistosomiasis infections when everybody adopts the optimal voluntary vaccination rate is constant (and higher than 4%) for η < 0.04. It follows that when using a combination of MDA and vaccination, the schistosomiasis can be eliminated, but in most cases it would be eliminated by MDA alone.
Uncertainty and sensitivity analysis
We used the Latin hypercube sampling with partial rank correlation coefficient(LHSPRCC) scheme (Blower & Dowlatabadi, 1994; Saltelli et al., 2004) to complete the uncertainty and sensitivity analysis. A full description of this method can be found in Marino et al. (2008).
Figure 5 shows the results of the analysis for $\mathcal{R}$, ν_{HI} and ν_{NE}. The uncertainty shows the distribution of model prediction among all the sampled parameter values. The most frequent values of $\mathcal{R}$ are between 0.6 and 2 with a peak around 1.2; we note that these are for vaccination rates between 0 and 0.1. The most frequent values of ν_{HI} are between 0 and 0.5 with most vaccination rates being below 0.25. Taken together, we can see that schistosomiasis would most likely be eliminated as long as the vaccination rates are 0.25 per year or higher, i.e., one would need to vaccinate the entire population at risk within 4 years. On the other hand, the optimal voluntary vaccination rate peaks between 0 and 0.03 and most values are less than 0.1. There is thus a big difference between the vaccination rate needed to eliminate schistosomiasis and the rate that is optimal for the individual.
The reproduction number is most sensitive to the aging rate μ_{1}, the death rates of snails, μ_{2}, cercariae μ_{4} and miracidia μ_{3}, vaccination rate ν, and the scaling factor for snail infection rate M_{0}; an increase of any of these parameters would cause decrease of $\mathcal{R}$. Similarly, increase of snail birth rate Λ_{2}, cercariae production γ_{2}, miracidia production γ_{1} or human birth rate Λ_{1} would cause $\mathcal{R}$ to decrease. Finally, an increase of β_{1} or β_{2}, i.e., increase of the rates parasites attack humans or snails, causes $\mathcal{R}$ to increase. The sensitivity analysis of ν_{HI} and ν_{NE} follows a very similar pattern. The exception is that the voluntary vaccination rate is most sensitive to the cost c, the higher the cost, the lower the vaccination rate.
The value of ν_{NE} is not as important as the actual prevalence of schistosomiasis when everybody adopts the optimal voluntary strategy. As seen from Fig. 6, in about 35% of the cases, the prevalence is below 1%; however, in about 65% of the cases, the prevalence is higher than 1%, meaning that schistosomiasis would not be eliminated as a public health concern. The prevalence is most sensitive to the cost of the vaccination (relative to the cost of schistosomiasis).
The situation improves significantly when the vaccination is accompanied by other control measures as seen in Fig. 7. In about 75% of the cases, the prevalence is below 1%. The prevalence is most sensitive to the rate of MDA treatment. The dependence on the other controls (elimination of snails or cercariae) is negligible.
Discussion
The big caveat of our quantitative results, though, is that, for simplicity, our model did not incorporate several important feature of schistosomiasis. First, the age is an important factor influencing the water contact and infection rates (Kura et al., 2021), but we considered it only marginally. To incorporate the agedependent water contact properly, we would have to stratify the human population by age groups. This stratification would also allow better tracking of the prevalence of the infections amongst school age children, which is crucial for the WHO’s elimination goal. The age groups play an important role from the logistical standpoint as well. Like MDA which is administered mostly to school age children (King et al., 2011), the vaccines would have to be administered before age 5 by incorporating into existing pediatric vaccine programs. Due to waning protection, the vaccination would have to be administered every 5 or so years. However, these aspects were not addressed by our model and thus more modeling effort need to be done to properly understand the effects of the vaccine.
Second, we assumed the vaccine offers 100% protection while the real efficacy will be likely around 90% (Zhang et al., 2018). Nevertheless, based on modeling of imperfect vaccine done for example in Reluga & Galvani (2011), Augsburger et al. (2023), Augsburger et al. (2022), as long as the vaccine is 85% or more effective, there are no big differences in model outcomes between perfect and imperfect vaccines. Furthermore, usage of S. mansonionly vaccine would likely not be acceptable in subSaharan Africa as there are regions where both S. mansoni and S. haematobium are endemic. A model that accounts for both species at the same time would be needed to better understand what to do in those regions.
Third, individuals eventually reach immunity (Kura et al., 2021; Wilkins et al., 1984) and this was omitted in our model that concentrated on the young population only. While the recovered compartment should be added to the later iterations of the model, we believe its addition would not significantly alter the results.
Our model can be further improved in several other ways. The underlying compartmental model can be made more realistic by (a) adding “exposed” compartments to human and intermediate hosts (such as in Anderson, Loker & Wearing, 2021), (b) considering the fact that infected humans release eggs rather than miracidia, and most importantly(c) specifically model the parasite load (such as in Woolhouse, Hasibeder & Chandiwana (1996)). Also, schistosomiasis endemicity exhibits a great variation when even neighboring villages show vastly different levels of parasite loads (Carabin et al., 2005). The distribution of schistosoma infections are highly over dispersed among hosts, even within age groups (Bundy, 1988); this can have implications on how effective the vaccination program is in reality. Incorporating some sort of structural modeling network to epidemics, for example as done in Hadjichrysanthou & Sharkey (2015) would be helpful. The game theory part of the model can be extended as follows. We assumed that every individual has the same risk of infection. However, the risk varies by age and by their behavioral pattern (MBra et al., 2018). Individuals thus have different risk perceptions (Poletti, Ajelli & Merler, 2011) and also base their decision on different social aspects (Xia & Liu, 2013). Therefore, it is often beneficial to use multiagentsimulation (MAS) methodology (Iwamura & Tanimoto, 2018; Kabir & Tanimoto, 2019; Kuga, Tanimoto & Jusup, 2019; Kabir & Tanimoto, 2020) which allows more flexibility and realism. Furthermore, our model assumed the risk of contracting the disease to be the only cost associated with notvaccination. If the vaccine is made mandatory, there can also be penalties for vaccine avoidance, possibly shrinking the gap between optimal voluntary vaccination level and the level required to achieve elimination. Finally, we assumed all individuals have perfect and full information. This is unlikely to happen in reality. However, the people will look up to their local leadership for advice and support. It is thus critical for the success of the vaccination campaign that the local leaders receive proper information about the disease and the available prevention methods.
Conclusions
We extended the compartmental model of schistosomiasis transmission (Gao et al., 2011) by adding the possibility of vaccination (Molehin, McManus & You, 2022; Stylianou et al., 2017) and applied the gametheoretic framework (Bauch & Earn, 2004). Unlike previous models of schistosomiasis transmission that focused on control and treatment at the population level, our model focuses on incorporating human behavior and voluntary individual vaccination.
We identified vaccination rates needed to achieve the herd immunity as well as optimal (from the individuals’ perspective) voluntary vaccination rates. We evaluated the prevalence of schistosomiasis for the scenario when everyone uses the optimal vaccination rates. We demonstrated that the prevalence remains too high (higher than 1%) unless the vaccination costs are sufficiently low. Thus, we can conclude that the voluntary vaccination alone may not be sufficient to eliminate schistosomiasis as a public health concern. When combining vaccination with MDA, the elimination is feasible; however, in such scenarios, the elimination would be possible by MDA alone.
We calibrated our model based on the data from literature. However, especially data related to transmission rates were lacking and we thus had to fit our model numerically to empirical data. We argue that there is an ongoing need to strengthen data collection and evaluation for decisionmaking (Toor et al., 2021). We also performed uncertainty and sensitivity analysis and showed that the results are relatively robust; the optimal voluntary vaccination (without MDA) will not eliminate schistosomiasis in at least 65% of the scenarios. With MDA, the situation is somewhat better, the elimination would occur in all but 25% of the scenarios.
The cost of the vaccine for the individual was an important factor determining whether or not voluntary vaccination can yield the elimination of schistosomiasis. When the cost is low (e.g., subsidized by governments or international help), the optimal voluntary vaccination rate is high enough that the prevalence of schistosomiasis declines under 1% and the disease is thus eliminated as a public health concern. Once the vaccine becomes available for public use, it will therefore be crucial to ensure that the individuals have cheap access to the vaccine.
Our main finding that voluntary vaccination alone may not be enough to eliminate schistosomiasis is not surprising. These conclusions had been already reached in a general scenario (Geoffard & Philipson, 1997) as well as demonstrated for specific diseases with a high cost of vaccination relative to the cost of the disease such as cholera (Kobe et al., 2018), Hepatitis B (Chouhan et al., 2020; Scheckelhoff et al., 2021), lymphatic filariasis (Rychtář & Taylor, 2022), polio (Cheng et al., 2020), or typhoid fever (AcostaAlonzo et al., 2020).