On the mathematical modeling of schistosomiasis transmission dynamics with heterogeneous intermediate host

Abstract

Schistosomiasis is a neglected disease affecting almost every region of the world, with its endemicity mainly experience in sub-Saharan Africa. It remains difficult to eradicate due to heterogeneity associated with its transmission mode. A mathematical model of Schistosomiasis integrating heterogeneous host transmission pathways is thus formulated and analyzed to investigate the impact of the disease in the human population. Mathematical analyses are presented, including establishing the existence and uniqueness of solutions, computation of the model equilibria, and the basic reproduction number (R₀). Stability analyses of the model equilibrium states show that disease-free and endemic equilibrium points are locally and globally asymptotically stable whenever R₀ < 1 and R₀>1, respectively. Additionally, bifurcation analysis is carried out to establish the existence of a forward bifurcation around R₀ = 1. Using Latin-hypercube sampling, global sensitivity analysis was performed to examine and investigate the most significant model parameters in R₀ which drives the infection. The sensitivity analysis result indicates that the snail's natural death rate, cercariae, and miracidia decay rates are the most influential parameters. Furthermore, numerical simulations of the model were done to show time series plots, phase portraits, and 3-D representations of the model and also to visualize the impact of the most sensitive parameters on the disease dynamics. Our numerical findings suggest that reducing the snail population will directly reduce Schistosomiasis transmission within the human population and thus lead to its eradication.

1. Introduction

Schistosomiasis which is also known as Bilharzia, Bilharziasis, and snail fever, is one of the neglected tropical diseases (NTD) common in some tropical and subtropical countries, in areas where there is no access to safe drinking water, adequate hygiene, and appropriate sanitation (sewage disposal). It is a water-based disease that is transmitted through contact with water infested with cercariae and it affects about 240 million people worldwide and more than 700 million people live in an area where the disease is endemic [1]. Schistosomiasis infection is also known as a parasite disease that is caused by schistosomiasis organisms (blood fluke). The two main forms of schistosomiasis infection are urogenital and intestinal which are caused by five schistosomiasis species responsible for human or animal infections namely, Schistosoma mansoni, Schistosoma intercalatum, Schistosoma haematobium, Schistosoma mekonai, and Schistosoma japonicum [1]. There are two stages of schistosomiasis infection namely; acute and chronic stages. Acute schistosomiasis infection referred to as Katayama Syndrome is a clinical manifestation of infection caused by trematode worms [2].

When non-immune persons are exposed to schistosomiasis infection, there is an incubation period symptoms of infection include abdominal pain, blood in stool or urine, wasting, anemia, dysuria [2]. Chronic schistosomiasis infection is caused by schistosome, which has severe long-term implications [3]. The symptoms of the chronic stage are growth retardation, renal failure, ascites, grand mal epilepsy, decreased fitness, anemia, portal hypertension, transverse myelitis, bladder carcinoma, and infertility [2]. Normally, snails harbor cercariae and shed the infective stage (cercariae) of the parasite (schistosome) in the freshwater [4]. Hence, there is an increase in the prevalence of Schistosomiasis as a result of freshwater habitats such as dams and irrigation, which enhances snail habitation and thus, increases transmission. Despite the progressive achievements in schistosomiasis control over the years, Schistosomiasis is still one of the most devastating parasitic diseases ravaging the entire world. Guiro et al. [5] proposed a mathematical model for schistosomiasis control, considering two population hosts, humans, and snails with delays. Their results suggest that an effective education campaign and reasonable coverage level of drug treatment will help in the control of Schistosomiasis. A mathematical model for the human-cattle-snail transmission of Schistosomiasis in Hubei Province, China was studied by Chen et al. [6]. Their results show that to control or eliminate Schistosomiasis in Hubei Province, China. Environmental factors need to be considered to lower (narrow) cattle-snail transmission. Since chemotherapy alone cannot end the spread of the causative parasite (schistosome), additional interventions and control strategies must be combined to lower re-infection, reduce the prevalence, and move toward the elimination of Schistosomiasis [7]. Many researchers have considered the use of molluscicides and different kinds of environmental modifications to control snails to reduce or eliminate schistosomiasis prevalence [7–9]. Interested readers can also see the following references [10–20] for other published works that have used mathematical models to gain insight into the transmission dynamics of Schistosomiasis and other infectious diseases in the human population and various suggestion for their control/elimination.

Given the aforementioned above, our study presents a mathematical modeling work that examines qualitative mathematical analysis for the prevention and control of the spread of Schistosomiasis through the control of snails. This is in addition to chemotherapy, molluscicides, and other environmental factors already studied in the previous research works and also the inclusion of the acute and chronic infected populations (compartments). The present study also investigates if mammals such as cattle contribute to the spread of Schistosomiasis infection within the population. The model described in this work considers three population hosts that is: humans, cattle, and snails.

The paper is organized as follows: Section 2 gives the Model formulation, Section 3 the Mathematical model analysis is presented, Section 4 presents the Sensitivity analysis of parameters in the model, Section 5 presents the Results of the numerical simulations and finally, Section 6 gives a brief Discussion and Conclusion.

2. Model formulation

The life cycle of the schistosome parasite is formulated mathematically. The populations considered are the human population, other mammal population, and snail population with the Cercariae C(t) and Miracida M(t). Cercariae C(t) is a larva worm that the infected snails shed into the aquatic environment, while the Miracida M(t) are eggs that the infected mammals (human and other mammals) shed in the stream when they come to the river for human activities such as fishing, swimming, fetching water, and/or drinking of water by other mammals. In this study, other mammals include cattle, mice, and dogs that the schistosome parasite can infect. The total human population, N_h(t), at time, t, is divided into the sub-populations; susceptible, S_h(t), acute infected, I_A(t), and chronic infected, I_ch(t) humans. The other mammal population, N_a(t), is subdivided into susceptible mammals, S_a and infected mammals, I_a, while the snail population, N_s(t) is made up of susceptible snails, S_s(t) and infected snails I_s(t). Human recruitment and natural death rates are assumed to be Λ_h and μ_h, respectively. The susceptible humans, S_h(t), become infected through contact with fresh water contaminated with Cercariae from infected snails. The force of infection is given by

where β_h is the human transmission rate, α_h is the adult parasite production rate for humans, C₀, is the half-saturation constant of Cercariae within the environment and ε₁ is the limitation of the growth velocity of the pathogens for Cercariae. Susceptible humans become acutely infected upon infection and enter into class, I_A(t). In the acutely infected sub-population, a proportion of them may recover with partial immunity depending on the mass drug administration program within the community and return to the susceptible class at a rate, ψ. On the other hand, some human may develop a chronic infection at a rate, k and move to class, I_ch. Chronically infected humans may succumb to infection at a rate, δ_h. Shedding of infection within the environment by acute infected and chronic infected humans is assumed to occur at rate N_Eτ₁γ_h and N_Eτ₁γ_hσ, respectively. Here, N_E is the number of eggs shed by mammals (both human and other mammals), τ₁ is the probability of eggs developing into Miracidium, γ_h is the infected human shedding rate and σ is a parameter that influences the shedding rate of the chronically infected humans. For the other mammals, the recruitment rate is assumed to be Λ_a while the natural death rate is μ_a. The susceptible mammals are infected through contact with fresh water contaminated with Cercariae from infected snails and the force of infection is given by

where β_a is the other mammal transmission rate, α_a, is the adult parasite production probability for other mammals. The other mammals shed the eggs that hatch into the water at a rate, N_Eτ₁γ_a where γ_a is the other mammals shedding rate. Death due to infection among other infected mammals is assumed to be at a rate, δ_a. Furthermore, susceptible snails are recruited at a rate, Λ_s and death due to natural causes of snails is at a rate, μ_s. Susceptible snails become infected by contact with miracidia from the viral shedding of infected humans and other mammals and the force of infection is given by

where β_s, is the snail transmission rate, α_s is the Miracidial penetration probability for snails, ε₂, is the limitation of the growth velocity for the Miracidium and M₀, is the half-saturation constant for Miracida within the environment. Death due to infection of snails is assumed to be at a rate, δ_s. Infected snails shed larva worm in the environment at a rate, τ₂γ_s where τ₂ is the density of Cercariae for snails and γ_s is the snail shedding rate. Finally, the decay rates for the Cercariae and Miracida are μ_c, μ_M, respectively. The study assumes that there is no immigration of infectious humans, other mammals, and snails within their populations. Also, susceptible humans, other mammals and snails are recruited at a constant rate. The dynamics of infection are presented in the model flow diagram in Figure 1. This leads to the following system of non-linear differential equations:

with initial conditions

Figure 1

All the parameters are assumed to be non-negative over the modeling time frame.

3. Mathematical methods

3.1. Invariant region

A biologically feasible population model such as model system (4) should always have positive solutions. To show that all solutions of model system (4) are non-negative, we let N_h, N_a and N_s be the total human population, other mammals population and snail populations, respectively, where, N_h(t) = S_h(t)+I_A(t)+I_ch(t), N_a(t) = S_a(t)+I_a(t), and N_s(t) = S_s(t)+I_s(t), respectively with the initial conditions N_h(0), N_a(0), N_s(0), M(0), and C(0).

Adding the human, other mammals and snail populations give the respective total populations as and . In the absence of infection which implies that the disease-related death rates for chronically infected human, δ_h, other infected mammals, δ_a and infected snail, δ_s are negligible, we assume that the total populations, N_h, N_a and N_s are asymptotically stable and yield

Solving (5) using the inequality theorem by Birkhoff and Rota [21] and integrating with the initial conditions N_h(0), N_a(0), N_s(0), C(0), and M(0) gives

As t → ∞, the total populations for humans, other mammals and snails approach , and , respectively. Since I_A, I_ch ≤ N_h(t), I_a ≤ N_a(t), and I_s ≤ N_s(t) it implies that the Cercariae (C) and Miracidia (M) populations can be written as

Applying the inequality theorem by Birkhoff and Rota [21] with the initial conditions C(0) = C₀₀ and M(0) = M₀₀ yields

and

As t → ∞, the Cercariae and Miracidia populations tends to respectively.

This shows that all the feasible solutions N_h(t), N_a(t), N_s(t), C(t) and M(t) are bounded and the region,

is positively invariant. Thus, we conclude that for all t>0, model system (4) is biologically feasible and well-posed in .

3.2. Positivity of the solution

We show that the solutions of the model system (4) are positive for all time by proving the following theorem:

Theorem 3.1. With positive initial conditions, the solution set of model system (4), [S_h(t),I_A(t),I_ch(t), S_a(t), I_a(t), M(t), S_s(t), I_s(t), C(t)], is positive for all time, t>0.

Proof. From the first equation of the model Equation (4), we have

Integrating with initial condition S_h(0), we obtain

with

In a similar way with and

Therefore, we conclude that the solution set [S_h(t), I_A(t), I_ch(t), S_a(t), I_a(t), M(t), S_s(t), I_s(t), C(t)] of model Equation (4) is positive for all t>0. This shows that the model is well-posed and biologically meaningful since the sub-population cannot be negative.

3.3. Existence of disease-free equilibrium

We examine the stability of the disease-free equilibrium state of the schistosomiasis model system (4) in terms of basic reproduction number. The disease-free equilibrium, E₀, is the equilibrium state in the absence of schistosomiasis disease in the environment. At equilibrium state,

This is solved simultaneously for the schistosomiasis model Equation (4) to yield the disease-free equilibrium,

3.3.1. The basic reproduction number

The basic reproduction number of the schistosomiasis model system (4) denoted by R₀, is a threshold parameter that is interpreted biologically as the average number of infected individuals created by one infected individual introduced in a susceptible population during the course of its infectious period [22, 23]. We compute R₀ using the next-generation matrix approach by Van den Driessche and Watmough [22]. Applying the approach, let and denote the Jacobian matrices for the rates of appearance of new infections and the transfer of infections in and out of the infected compartments, respectively at DFE, E₀, so that

where

The basic reproduction number, R₀, for the model system (4) is defined as the spectral radius of the operator, , where is the inverse matrix of . This gives

for

Here, the reproduction numbers, R_A, R_ch, R_m, and R_s, are the reproduction numbers for the acutely infected humans, chronic infected humans, other infected mammals and infected snails contribute, respectively. This shows that the basic reproduction number, R₀ of the schistosomiasis model system (4) comprises of four parts, that is, R_A, R_ch, R_m, R_s. Importantly, whenever R₀ < 1, it implies that the disease will fizzle out in the population and whenever R₀>1, it means that the disease will persist in the population.

3.4. Local stability of the disease-free equilibrium

We prove the local stability of the disease-free equilibrium, E₀, for the Schistosomiasis model system (4) using the Jacobian method.

Theorem 3.2. The disease-free equilibrium state, E₀, of the schistosomiasis model system (4) is locally asymptotically stable when R₀ < 1 and unstable forR₀>1.

Proof. We compute the Jacobian matrix of the model system (4) at the equilibrium state, E₀, as follows.

where f, g, h, q, A, B, D are defined in (12). The eigenvalues of (15) are λ₁ = −μ_h, λ₂ = −μ_a, λ₃ = −μ_s and the roots of the characteristics equation

where

With the definitions of reproduction numbers in (14), we have by Routh-Hurwitz criteria in Kim et al. [24] and the conditions in Heffernan et al. [25] that the polynomial in (16) has negative real roots provided that P>0, Q>0, R>0, U>0, V>0, W>0, PQ>R, QU>W, UV>RW and this happens when R₀ < 1. With negative real eigenvalues, it implies that the disease-free equilibrium state, E₀, for Schistosomiasis model system (4) is locally asymptotically stable if R₀ < 1 otherwise it is unstable.

3.5. Existence of endemic equilibrium

For the endemic equilibrium, let

in model system (4). We have at equilibrium state that

Substituting , and C^* into and and solve simultaneously, we obtain

with as solution of the polynomial

where

Solving (20) for , it can be seen that one of the solutions is which represents the disease-free equilibrium state, while the other solution of the polynomial (20) is the endemic equilibrium state of the model system (4) that is given by

and it is positive whenever R₀>1. We further use Descarte's rule of signs to discuss the existence of possible positive roots of Equation (18) especially when R₀ < 1. The results are summarized in Table 1. Hence, the following theorem.

Theorem 3.3. System (4) has the following endemic equilibrium:

• If R₀>1, then the system has one positive unique endemic equilibrium state.

• If R₀ < 1, then the system has zero or two endemic equilibria.

In case (4), there is possibility of backward bifurcation of the model system (4) as it may contain endemic equilibrium with DFE when R₀ < 1. This is further determine by carrying out a bifurcation analysis in the following subsection.

3.6. Bifurcation analysis

The stability of endemic equilibrium state is established using Center Manifold Theory by Chavez et al. [26]. The theory shows the direction of bifurcation, whether it is a forward or backward bifurcation. A forward bifurcation indicates that the disease-free and endemic equilibrium states are locally asymptotically stable if R₀ < 1 and R₀>1, respectively, which implies the existence of global stability of the equilibrium states. A backward bifurcation means there is a co-existence of disease-free and endemic equilibrium when R₀ < 1. Therefore, we apply the theory in Chavez et al. [26] by re-writing the state variables as follows

such that the model system (4) can be rewritten as

where is the right hand side of model system (4).

Let the bifurcation parameter be at R₀ = 1 so that the Jacobian matrix, J(E₀), of (15) will have a simple zero eigenvalue and negative real eigenvalues. The left and right eigenvectors associated with a simple zero eigenvalue of the Jacobian matrix (15) are w = (w₁, w₂, w₃, w₄, w₅, w₆, w₇, w₈, w₉) and v = (v₁, v₂, v₃, v₄, v₅, v₆, v₇, v₈, v₉), respectively, where

We state and prove the following theorem

Theorem 3.4. Model system (4) undergoes a forward bifurcation atR₀ = 1.

Proof. To determine the direction of the bifurcation, Theorem 4.1 in Castillo-Chavez and Song [27] is used to compute the bifurcation coefficients, m and n, of model system (4). The non-zero associated partial derivatives of F at the disease-free equilibrium, E₀, are

Now applying the concept of Castillo-Chavez and Song [27], the value of m gives

Upon substituting the eigenvectors and partial derivatives, we obtain

where

Also for n, we have

A forward bifurcation exists since m < 0 in (24) and n>0 in (26). Therefore, a forward bifurcation exists for the model (4) at R₀ = 1. This completes the proof.

With the existence of forward bifurcation, a unique endemic equilibrium, E₁, exists and is locally asymptotically stable whenever R₀>1 and unstable when R₀ < 1. Hence, there is possibility of global stability of the equilibrium states. A forward bifurcation plot is shown in Figure 2.

Figure 2

3.7. Global stability of the model

In this subsection, the global behavior of the disease-free and endemic equilibria of the system (4) is established by constructing a suitable Lyapunov functions and applying LaSalle's invariance principle.

3.7.1. Global stability of disease-free equilibrium, E₀

Theorem 3.5. The disease-free equilibrium, E₀, of the schistosomiasis model system (4) is globally asymptotically stable whenR₀ < 1 otherwise unstable.

Proof. Constructing a Lyapunov function using a matrix-theoretic method by [28] gives

where the letters, A, B, and D are defined in (12).

Taking the derivative of Equation (27) along the trajectory yields

Expanding Equation (28) and simplifying in terms of reproduction numbers (R_A, R_ch, R_m, R_s) gives

With , , , we have

Since R_A+R_ch+R_m ≤ 1, R_s ≤ 1 from the definition of the reproduction number in (13) and (14), this yields L′ ≤ 0 if R₀ < 1. Hence by LaSalle's invariance principle, the disease-free equilibrium is globally asymptotically stable if R₀ ≤ 1. This completes the proof.

3.7.2. Global stability of endemic equilibrium, E₁

To prove the stability of the endemic equilibrium (E₁) we state and prove the following theorem

Theorem 3.6. The endemic equilibrium, E₁, of the schistosomiasis model system (4) is globally asymptotically stable wheneverR₀>1 and ψ = 0.

Proof. The global stability of the endemic equilibrium is proved by constructing the Volterra-Lyapunov function

where , , A, B, and D are defined in (12) and (18).

Taking the time derivative of along the solutions of the schistosomiasis model Equation (4) yields

At endemic equilibrium state, we have

Substituting Equation (32) into (31) and expanding and simplifying gives

Using the assumption that ψ = 0 and , , we have by Arithmetic-Geometric Theorem that

This implies that L′ ≤ 0 in the region, Z. Applying LaSalle's invariance principle, the endemic equilibrium state of the schistosomiasis model system (4) is globally asymptotically stable in Z if R₀>1 and ψ = 0. This means that with different initial conditions at ψ = 0, the global stability of the schistosomiasis model system (4) will always converge to the disease-free and endemic equilibrium states whenever R₀ < 1 and R₀>1, respectively. The numerical simulations are carried out in the next section.

4. Sensitivity analysis

4.1. Parameter values

To perform numerical simulations, parameters for the model are mainly sourced from literature, as shown in Table 2. Parameters that are represented by probabilities are assumed to be within the range of [0 − 1]. Any parameters that are outside the prescribed range are properly sampled from a uniform distribution within ranges defined in Table 2.

4.2. Sensitivity analysis results

Sensitivity analysis examines the relative change of the mathematical model output (variable) when the model parameter(s) are altered or changed. It helps in determining and identifying the model parameters that require special attention. To carry out the sensitivity analysis of the basic reproduction number R₀, we adopt the Latin Hypercube Sampling (LHS) scheme and Partial Rank Correlation Coefficients (PRCCs) technique used by Blower and Dowlatabadi [39] to compute and identify the biological implication of each model parameter on R₀. See also the following published articles [40–45] with similar analyses for different epidemic models where the LHS techniques have been applied. We performed 1,000 simulations per run and examined and evaluated the PRCC of the model parameters concerning R₀. The PRCC values indicate the degree of monotonicity between the parameters of the model and R₀. Thus, juxtaposing the PRCC values gives us a clear picture on the contribution of each parameter of the model on R₀. Table 3 and the Tornado plot of Figure 3 give the model parameters with their PRCCs. The parameters with the potential of increasing(decreasing) the value of R₀, are those with positive(negative) PRCCs. Thereby, increasing parameters, β_s, β_h, β_a, γ_h, γ_a, γ_s, τ₁, τ₂, increases the value of R₀ which in turn increase the spread of schistosomiasis disease within the population (see Table 2 and Figure 3). On the other hand, parameters μ_c, μ_m, μ_s with negative PRCCs influence the reduction of the burden of schistosomiasis disease within the population. When these parameters are increased, we observe a reduction in the value of R₀, hence the endemicity of the Schistosomiasis within the population is also reduced. The natural death rate for the snail population is the most sensitive parameter with respect to R₀ with a PRCC value of 0.776046728. The implication of this is that increasing snail death rate reduces the spread of Schistosomiasis within the population. To examine the significance of the model parameters, we compute the p-value of PRCC of each model parameter using the Fisher Transformation method as shown in Table 3. Table 3 reveals that the parameters, β_s, β_h, γ_h, γ_s, μ_c, μ_m, μ_s, τ₁, τ₂ have significant p-values while the parameters, β_a, γ_a, have insignificant p-values. This implies that the shedding and transmission rates of the other mammals, γ_a and β_a, are non-monotonically related to R₀, though they can still produce changes in the transmission dynamics of schistosomiasis infection.

Figure 3

The pairwise comparison of the significant parameters of Table 3, whose p < 0.05 is given in Tables 5, 6. The pairwise comparison is carried out to observe the model parameters with notable influence on increasing the burden of the schistosomiasis infection in the population, as well as to establish whether there exists any difference between the processes describing the compared parameters. We presented the outcome of the pairwise PRCC comparison between the unadjusted p-values (Table 4) and the false discovery rate (FDR) adjusted p-values in Tables 5, 6. According to Table 6, whenever the p-values of the compared pair of significant parameters are < 0.05, it is interpreted to be significant (TRUE) otherwise it is insignificant (FALSE). It is also important to note from Table 6 that the pairs with significantly different (TRUE) are more sensitive parameters.

The effect of the most sensitive parameters on R₀ is presented in Figures 3, 4 as scatter plots to support results in Table 2. These plots show how parameters, β_s, γ_s, γ_h, τ₁, μ_c, μ_m, μ_s, τ₂ monotonically increase or decrease R₀. The most noticeable result is the scatter plot for snail death μ_s which is the most sensitive parameter. The sensitivity result reveals that the spread of schistosomiasis infection will be controlled and eradicated in the population if more snails die. If the Miracidia and Cercariae (μ_m, μ_c) are cleared within the population, we will observe a reduction in the contact rates of human to Cercariae (β_h), snail to Miracidia (β_s), probability of egg developing into a Miracidium (τ₁) and the density of Cercariae in the environment (τ₂).

Figure 4

5. Numerical simulation results

5.1. Time series trajectories

To support the analytical results presented in this paper, we present numerical simulations in the form of time series plots and phase plots on the dynamics of Schistosomiasis. The simulations presented are produced using MATLAB software. Figures 5A,B, present the time series plots for the acute infected class I_A(t), chronic infected class, I_ch(t), infected animals I_a(t), and infected snails, I_a(t), while varying the snail death rate, μ_s. The snail death rate, μ_s is considered for simulation as it is the most sensitive parameter. It is observed that when the death rate of the snail population is increasing, the peak of I_A(t), I_ch(t), and I_a(t), shifts to the right in a reduced manner. On the other hand, Figures 6A,B present the dynamics of Miracidia, M(t) and Cercariae, C(t). We observe that an increase in the snail death rates leads to a decrease in the population of the Miracidia. Also, the populations of I_s(t) and C(t) reduce drastically as the snail population's death rate increases. This indicates that reducing the snail population reduces the time it will take to infect human and mammal populations, including Cercariae that will be cleared in the population. This eventually reduces the number of Schistosomiasis cases within the population.

Figure 5

Figure 6

5.2. Phase portraits

Figures 7, 8 present the phase portraits of the infection sub-populations when R₀>1 and R₀ < 1, respectively. When R₀ < 1, it can be seen that solutions tend to disease-free equilibrium over time. This means that Schistosomiasis infection can be contained within the population over time provided that R₀ < 1. On the other hand when R₀>1, the solutions lead to an endemic equilibrium as seen by a decrease in the susceptible population and an increase in the infected population. It is also observed that as Miracidia initially increase, Cercariae increases as well and reaches a peak. Thereafter, Miracidia decreases gradually as Cercariae continues to increase. This means that Schistosomiasis infection will persist within the population whenever R₀>1.

Figure 7

Figure 8

5.3. Effects of parameters on R₀

Using contour plots, Figures 9A,B depict that R₀ attains a lower value when the probability of eggs developing into Miracidium is minimal and the snail death rate is high as shown by the contour plot and the 3-D representation. Figure 9C shows R₀ vs. varying the natural decay rates of Miracidium, (μ_m) and Cercariae (μ_c). It can be seen that if the decay rates are effectively increased, the reproduction number, R₀ can be reduced to below one, implying that the disease-free equilibrium can be achieved. The implication of Figures 9A–C is that there is a possibility of approaching a disease-free equilibrium if these parameters are controlled accordingly. However, this may be difficult to attain due to the lack of proper control measures that can simultaneously increase the decay rates. Generally, it can be seen that the control of snails is very important in the reduction of the spread of Schistosomiasis within the population. Therefore, measures need to be put in place so that all breeding areas for snails are identified and dealt with appropriately without harming other aquatic organisms.

Figure 9

6. Conclusion

In this paper, a deterministic mathematical model of transmission dynamics of Schistosomiasis disease that incorporates acute and chronic infected humans is formulated. The main objective of the model is to study the impact of snail control on the dynamics of Schistosomiasis. The model analysis comprised of establishing the invariant region, the positivity of solutions, the existence of disease-free and endemic equilibrium, bifurcation analysis and computation of the basic reproduction number, R₀ using the next-generation approach. We proved that whenever R₀ < 1, the disease-free equilibrium is globally asymptotically stable, which means that the disease will fizzle out within the population over time. The endemic equilibrium was found to be globally asymptotically stable implying that if not controlled, Schistosomiasis will persist within the population. We also carried out a sensitivity analysis of the parameters governing the basic reproduction number, R₀, using the Latin Hypercube Sampling (LHS) scheme-Partial Rank Correlation Coefficient (PRCC) technique. The results of the sensitivity analysis revealed that the spread of schistosomiasis infection can be controlled and eradicated in the population by increasing the natural death rate of the snail population. It also established that there is need to increase the decay rates of Miracidia and Cercariae within the environment by applying a control measure. Therefore, the results of our study effectively conclude that the most effective method to control the transmission dynamics and spread of schistosomiasis disease is by rapidly killing of infected snails and applying control measures that will clear the parasite, miracidia and cercariae, in the environment.

One limitation of our model is that it was not fitted to real-life data. Instead, our model was parameterized using accurate parameter estimates from existing literature on Schistosomiasis dynamics with inclusive heterogeneity in the disease transmission. Secondly, the assumption of mammals could be made specific by considering only cattle in the modeling framework. To bypass this complexity, we carried out a global sensitivity analysis to vary the parameter values and establish the uncertainty in the model parameters. Notwithstanding this limitation, we maintain that our findings remain valid and applicable if followed by policymakers to help eradicate the spread of Schistosomiasis.

Future work may look at an optimal control analysis of the model presented in the presence of different control strategies. The model may also be extended to include infection dynamics in the presence of snail age stratification.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

• 1.

  WHO. WHO Guideline on Control and Elimination of Human Schistosomiasis. WHO (2022).

  • Google Scholar

• 2.

  RossAGVickersDOldsGRShahSMMcManusDP. Katayama syndrome. Lancet Infect Dis. (2007) 7:218–24. 10.1016/S1473-3099(07)70053-1

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 3.

  ShiffC. The importance of definitive diagnosis in chronic schistosomiasis, with reference to Schistosoma haematobium. J Parasitol Res. (2012) 2012:761269. 10.1155/2012/761269

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 4.

  Costa de LimaeMPochaRFilhocouraPKatzN. A 13-years follow-up of treatment and snail control in the endemic area for Schistosoma mansoni in Brazil, incidence of infection and re-infection. Bull World Health Organ. (1993) 71:197–205.

  • Pubmed Abstract
  • Google Scholar

• 5.

  GuiroAOuaroSTraoreA. Stability analysis of a schistosomiasis model with delays. Adv Differ Equat. (2013) 2013:1–15. 10.1186/1687-1847-2013-303

  • CrossRef
  • Google Scholar

• 6.

  ChenZZouLShenDZhangWRuanS. Mathematical modelling and control of Schistosomiasis in Hubei Province, China. Acta Trop. (2010) 115:119–25. 10.1016/j.actatropica.2010.02.012

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 7.

  SokolowSHWoodCLJonesIJSwartzSJLopezMHsiehMHet al. Global assessment of schistosomiasis control over the past century shows targeting the snail intermediate host works best. PLoS Neglect Trop Dis. (2016) 10:e0004794. 10.1371/journal.pntd.0004794

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 8.

  McCulloughFGayralPDuncanJChristieJ. Molluscicides in schistosomiasis control. Bull World Health Organ. (1980) 58:681.

  • Google Scholar

• 9.

  LardansVDissousC. Snail control strategies for reduction of schistosomiasis transmission. Parasitol Today. (1998) 14:413–7. 10.1016/S0169-4758(98)01320-9

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 10.

  ZhangHHarvimPGeorgescuP. Preventing the spread of schistosomiasis in Ghana: possible outcomes of integrated optimal control strategies. J Biol Syst. (2017) 25:625–55. 10.1142/S0218339017400058

  • CrossRef
  • Google Scholar

• 11.

  ZhangTZhaoXQ. Mathematical modeling for schistosomiasis with seasonal influence: a case study in Hubei, China. SIAM J Appl Dyn Syst. (2020) 19:1438–71. 10.1137/19M1280259

  • CrossRef
  • Google Scholar

• 12.

  BakareENwozoC. Mathematical analysis of malaria-schistosomiasis coinfection model. Epidemiol Res Int. (2016) 2016:3854902. 10.1155/2016/3854902

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 13.

  AkanniJOOlaniyiSAkinpeluFO. Global asymptotic dynamics of a nonlinear illicit drug use system. J Appl Math Comput. (2021) 66:39–60. 10.1007/s12190-020-01423-7

  • CrossRef
  • Google Scholar

• 14.

  MustaphaUTMusaSSLawanMAAbbaAHincalEMohammedMDet al. Mathematical modeling and analysis of schistosomiasis transmission dynamics. Int J Model Simul Sci Comput. (2021) 12:2150021. 10.1142/S1793962321500215

  • CrossRef
  • Google Scholar

• 15.

  WuYFLiMTSunGQ. Asymptotic analysis of schistosomiasis persistence in models with general functions. J Frankl Instit. (2016) 353:4772–84. 10.1016/j.jfranklin.2016.09.012

  • CrossRef
  • Google Scholar

• 16.

  RonohMChiroveFPedroSATchamgaMSSMadubuezeCEMadubuezeSCet al. Modelling the spread of schistosomiasis in humans with environmental transmission. Appl Math Modell. (2021) 95:159–75. 10.1016/j.apm.2021.01.046

  • CrossRef
  • Google Scholar

• 17.

  QiLCuiJAGaoYZhuH. Modeling the schistosomiasis on the islets in Nanjing. Int J Biomath. (2012) 5:1250037. 10.1142/S1793524511001799

  • CrossRef
  • Google Scholar

• 18.

  AbimbadeSFOlaniyiSAjalaOA. Recurrent malaria dynamics: insight from mathematical modelling. Eur Phys J Plus. (2022) 137:1–16. 10.1140/epjp/s13360-022-02510-3

  • CrossRef
  • Google Scholar

• 19.

  HethcoteHW. The mathematics of infectious diseases. SIAM Rev. (2000) 42:599–653. 10.1137/S0036144500371907

  • CrossRef
  • Google Scholar

• 20.

  AcheampongEOkyereEIddiSBonneyJHAsamoahJKKWattisJAet al. Mathematical modelling of earlier stages of COVID-19 transmission dynamics in Ghana. Results Phys. (2022) 34:105193. 10.1016/j.rinp.2022.105193

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 21.

  BirkhoffGRotaGC. Ordinary Differential Equations. New York, NY: John Wiley & Sons (1978).

  • Google Scholar

• 22.

  Van den DriesschePWatmoughJ. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math Biosci. (2002) 180:29–48. 10.1016/S0025-5564(02)00108-6

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 23.

  DiekmannOHeesterbeekJAPMetzJA. On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations. J Math Biol. (1990) 28:365–82. 10.1007/BF00178324

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 24.

  KimJHSuWSongYJ. On stability of a polynomial. J Appl Math Inform. (2018) 36:231–6. 10.14317/jami.2018.231

  • CrossRef
  • Google Scholar

• 25.

  HeffernanJMSmithRJWahlLM. Perspectives on the basic reproductive ratio. J R Soc Interface. (2005) 2:281–93. 10.1098/rsif.2005.0042

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 26.

  ChavezCCFengZHuangW. On the computation of R0 and its role on global stability. In: Castillo-Chavex C Blower S Driessche P Kirschner D Yakubu A-A, editors. Mathematical Approaches for Emerging and Re-emerging Infection Diseases: An Introduction, Vol. 125. New York, NY: Springer (2002) p. 31–65. 10.1007/978-1-4757-3667-0_13

  • CrossRef
  • Google Scholar

• 27.

  Castillo-ChavezCSongB. Dynamical models of tuberculosis and their applications. Math Biosci Eng. (2004) 1:361. 10.3934/mbe.2004.1.361

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 28.

  ShuaiZvan den DriesscheP. Global stability of infectious disease models using Lyapunov functions. SIAM J Appl Math. (2013) 73:1513–32. 10.1137/120876642

  • CrossRef
  • Google Scholar

• 29.

  KalindaCMushayabasaSChimbariMJMukaratirwaS. Optimal control applied to a temperature dependent schistosomiasis model. Biosystems. (2019) 175:47–56. 10.1016/j.biosystems.2018.11.008

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 30.

  FengZEppertAMilnerFAMinchellaD. Estimation of parameters governing the transmission dynamics of schistosomes. Appl Math Lett. (2004) 17:1105–12. 10.1016/j.aml.2004.02.002

  • CrossRef
  • Google Scholar

• 31.

  ChiyakaETGARIRAW. Mathematical analysis of the transmission dynamics of schistosomiasis in the human-snail hosts. J Biol Syst. (2009) 17:397–423. 10.1142/S0218339009002910

  • CrossRef
  • Google Scholar

• 32.

  AbokwaraAMadubuezeCE. The role of non-pharmacological interventions on the dynamics of schistosomiasis. J Math Fund Sci. (2021) 53:243–260. 10.5614/j.math.fund.sci.2021.53.2.6

  • CrossRef
  • Google Scholar

• 33.

  KanyiEAfolabiASOnyangoNO. Mathematical modeling and analysis of transmission dynamics and control of schistosomiasis. J Appl Math. (2021) 2021:6653796. 10.1155/2021/6653796

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 34.

  MangalTDPatersonSFentonA. Predicting the impact of long-term temperature changes on the epidemiology and control of schistosomiasis: a mechanistic model. PLoS ONE. (2008) 3:e1438. 10.1371/journal.pone.0001438

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 35.

  LiYTengZRuanSLiMFengX. A mathematical model for the seasonal transmission of schistosomiasis in the lake and marshland regions of China. Math Biosci Eng. (2017) 14:1279. 10.3934/mbe.2017066

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 36.

  XiangJChenHIshikawaH. A mathematical model for the transmission of Schistosoma japonicum in consideration of seasonal water level fluctuations of Poyang Lake in Jiangxi, China. Parasitol Int. (2013) 62:118–26. 10.1016/j.parint.2012.10.004

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 37.

  GaoSJCaoHHHeYYLiuYJZhangXYYangGJet al. The basic reproductive ratio of Barbour's two-host schistosomiasis model with seasonal fluctuations. Parasites Vectors. (2017) 10:1–9. 10.1186/s13071-017-1983-1

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 38.

  DingCQiuZZhuH. Multi-host transmission dynamics of schistosomiasis and its optimal control. Math Biosci Eng. (2015) 12:983. 10.3934/mbe.2015.12.983

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 39.

  BlowerSMDowlatabadiH. Sensitivity and uncertainty analysis of complex models of disease transmission: an HIV model, as an example. Int Stat Rev. (1994) 62:229–43. 10.2307/1403510

  • CrossRef
  • Google Scholar

• 40.

  HerdichoFFChukwuWTasmanH. An optimal control of malaria transmission model with mosquito seasonal factor. Results Phys. (2021) 25:104238. 10.1016/j.rinp.2021.104238

  • CrossRef
  • Google Scholar

• 41.

  ChukwuCNyabadzaF. A theoretical model of listeriosis driven by cross contamination of ready-to-eat food products. Int J Math Math Sci. (2020) 2020:9207403. 10.1155/2020/9207403

  • CrossRef
  • Google Scholar

• 42.

  ChazukaZMoremediGM. In-host dynamics of the human papillomavirus (HPV) in the presence of immune response. In:MondainiRP, editor. Trends in Biomathematics: Chaos and Control in Epidemics, Ecosystems, and Cells: Selected Works from the 20th BIOMAT Consortium Lectures, Rio de Janeiro, Brazil, 2020. Cham: Springer (2021). p. 79–97. 10.1007/978-3-030-73241-7_6

  • CrossRef
  • Google Scholar

• 43.

  ChukwuCWMushanyuJJugaMLFatmawathi. A mathematical model for co-dynamics of listeriosis and bacterial meningitis diseases. Commun Math Biol Neurosci. (2020) 2020:ID83.

  • Google Scholar

• 44.

  MadubuezeCEChazukaZ. An optimal control model for the transmission dynamics of Lassa fever. Res Squ. [Preprint]. 10.21203/rs.3.rs-1513399/v1

  • CrossRef
  • Google Scholar

• 45.

  ChukwuCJugaMChazukaZMushanyuJ. Mathematical analysis and sensitivity assessment of HIV/AIDS-listeriosis co-infection dynamics. Int J Appl Comput Math. (2022) 8:1–21. 10.1007/s40819-022-01458-3

  • Pubmed Abstract
  • CrossRef
  • Google Scholar