Mathematical Global Dynamics and Control Strategies on Echinococcus multilocularis Infection.

Echinococcus multilocularis, a major cause of echinococcosis in human, is a parasitic sylvatic disease between two major hosts in a predator-prey relation. A new model for the transmission dynamics of Echinococcus multilocularis in the population of red foxes and voles with environment as a source of infection is formulated and rigorously analyzed. The model is used to access the impact of treatment on red foxes and environmental disinfection as control strategies on the disease dynamics. The control reproduction number is computed and is used to rigorously prove the local and global dynamics of models' equilibria. Using available data on Echinococcus, elasticity indices and partial rank correlation coefficients of control reproduction number and cumulative new cases in red foxes and voles are computed. Parameters that have high influence locally and globally are identified. Numerical experiments indicate that administering disinfection of environment only induces more positive impact than applying treatment only on red foxes in controlling the infection. Generally, interventions towards treating red foxes and environmental disinfection could be sufficient in tackling transmission of disease in the populations.

1. Introduction

Echinococcus multilocularis (EM) is a parasitic taeniid tapeworm and one of the six species of the Echinococcus genus. The disease is sylvatic and has indirect life cycle between two major hosts in a predator-prey interaction [1-3]. Adult EM inhabits the small intestine of canines (such as red foxes) which are regarded as definitive hosts and produces eggs that are released to the environment, typically through faeces. After oral ingestion of eggs by rodents (such as voles), regarded as intermediate hosts, a larval stage (metacestode) develops in any of the internal organs (liver, kidney, heart, etc). The mature metacestodes are capable of producing numerous protoscoleces, each having the potential to develop into an adult EM when a definitive host preyed on an intermediate host, and the cycle continues [1, 2, 4, 5]. In some parts of the world, other wild canids (such as coyotes, raccoon dogs, and wolves) can also serve as definitive hosts while other animals (like sheep, domestic dogs, cats, and rats) can be regarded as intermediate hosts [1, 4, 6]. The parasite (EM) causes alveolar echinococcosis in humans regarded as accidental hosts and is characterized with a tumour-like, destructive growth with the potential of causing high fatality rate [1, 4, 7]. The disease was initially confined to a certain part of the globe; however, as researches indicated, it spread all over the globe especially in rural nomadic communities that are economically less privileged, geographically and/or behaviourally detached to certain extent from healthcare systems [8-12]. Furthermore, the disease has high prevalence in red foxes population (1%–76.7%) [3, 5, 13–16] and low prevalence in rodents (voles) (0.4%–30%) [3, 5, 13, 15]. Due to the prevalence of the disease in intermediate hosts and the difficulties involved in treating definitive hosts in a given community, control can take longer times and in many cases may last indefinitely.

Modelling approaches, and in particular mathematical modelling, can give an insight into the biology and epidemiology of diseases in terms of revealing facts on data gaps, understanding the interaction between organisms, predicting future, control, and quantifying certain unmeasurable quantities such as force of infection, basic reproduction number, etc [2]. Mathematical models of Echinococcus multilocularis have been developed in literature to attend to must of the above assertions. For instance, in [5], Ishikawa et al. proposed a mathematical model that described the transmission of Echinococcus multilocularis in both the definitive (foxes) and intermediate host (voles) populations in Hokkaido. They quantitatively studied seasonal transition in the prevalence of EM in foxes and the risk of infection with human alveolar echinococcosis. Roberts and Aubert [3] developed a simple mathematical model in an attempt to determine the likely effect of combining treatment for EM infection in red foxes and voles in France with the existing vaccination campaign against rabies. It was shown that if the prevalence of the EM in foxes is low in a particular region, the parasite can be eradicated or controlled. Wang et al. [17] proposed a model for the transmission of echinococcosis in dogs, livestock, and human populations to explore effective control and preventive measures in Xianjing. Basic reproduction number, on which the dynamics of model was completely determined, has been estimated, and sensitivity analysis was carried out based on data relevant to the study area. In [18], a model that took into account, the contribution of domestic and stray dogs on the transmission of the parasite in humans was proposed. The global dynamics revealed that, without disposing the stray dogs, the disease became endemic even if the domestic dogs are controlled.

It is a known fact that environment aids the transmission dynamics of Echinococcus multilocularis, as discussed in the above reviews. However, this important component is neglected in the modelling process of the disease. In this paper, we establish a noble mathematical model for the transmission dynamics of Echinococcus multilocularis within the population of foxes as definitive hosts and the voles as intermediate hosts with concentration of parasites in the environment as a source of infection for intermediate hosts. Considering the high prevalence of EM in red fox population (1–76.7%) compared to rodent population as reported in previous studies, we incorporate treatment as control strategy on the infected red foxes. As a result of the treatment, we assume that recovered red foxes will acquire long-time immunity and as such will not return to susceptible. Furthermore, disinfection or cleaning of environment to reduce the concentration of the disease is also incorporated as a second control strategy. We carry out rigorous analysis on the computation of the basic control reproduction number, a threshold quantity used, to determine the existence and stability dynamics of equilibria. Furthermore, using data available from literature, we conduct elasticity indices of parameters on the control reproduction number and global sensitivity analysis using partial rank correlation coefficients of control reproduction number and cumulative new infections on the two hosts populations. Numerical simulations are conducted to support analytic results and effects of control strategies on the model. This paper is organized as follows. The model formulation, equations, and flow diagram are presented in Section 2. Basic properties of the model on existence, uniqueness, positivity, and boundedness of solutions are discussed in Section 3. Furthermore, existence and global stabilities with systematic calculation of control reproduction number are presented in Section 4. Numerics which comprise of elasticity index, global sensitivity analysis, numerical simulations, and effects of control strategies are presented in Section 5. Finally, we present concluding remarks in Section 6.

2. Model Formulation

The total population of red foxes, which is assumed constant (birth and death rates, μ, are assumed equal) in the environment at time t, denoted by N(t) is divided into susceptible (S(t)), exposed (E(t)), infected (I(t)), and recovered (R(t)) subpopulations so that

The susceptible population is increased by recruitment of foxes by birth or immigration at rate μN and is decreased when it preyed with searching efficiency s on an infected vole containing protoscoleces in hydatid cysts [5] with probability p of becoming infectious. The exposed fox population is increasing by the number of susceptibles that preyed on infected voles and is decreasing by progression to infected population and natural deaths, at rates α and μ, respectively. The infected fox population is increased by progression of exposed foxes and decreased as a result of treatment and natural deaths at rates ξ and μ, respectively. The population of recovered is increased by the treated infected foxes and decreased by natural death at rates ξ and μ, respectively. We assume here that treated red foxes have either permanent or long-lasting immunity to the parasite and hence will not return to susceptible population.

Similarly, the total population of voles, also assume constant (birth and death rates, μ, are assumed equal) in the environment at time t, denoted by N(t) is subdivided into susceptible (S(t)), exposed (E(t)), and infected (I(t)) subpopulations so that

The population of susceptible voles is increased by birth at rate μN and decreased by infection from the concentration of parasites in the environment at the rate β and natural deaths at rate μ, which is also applicable to all the subpopulations of voles. Furthermore, the concentration of Echinococcus in the environment B(t) is increased by shedding of the parasites by infected foxes at rate η and decreased by disinfection or cleaning of environment at rate μ. Based on the above descriptions, the model can be described completely by the following system of ordinary differential equations, which follow from the schematic diagram shown in Figure 1:with initial conditions

Figure 1

Flow diagram of the transmission dynamics of models (3)–(10).

3. Basic Properties of the Model

The model (3)–(10) monitors the dynamics of red foxes and voles populations, and all its associated parameters are assumed nonnegative. Hence, we now present the basic results for the properties of the model.

Theorem 1 .

The following region is positively invariant for the model (3): Ω=Ω × Ω × Ω ⊂ ℝ+4 × ℝ+3 × ℝ+, where Ω={(S, E, I, R) ⊂ ℝ+4 : S+E+I+R=N}, Ω={(S, E, I) ⊂ ℝ+3 : S+E+I=N} and Ω={B ⊂ ℝ+ : B ≤ (ηN/μ)}.

Proof

The detailed proof of Theorem 1 is presented in Appendix A.

Equations (A.14)–(A.18) establish the boundedness of total populations for red foxes, voles, and concentration of parasites, respectively, and by extension verifies the boundedness of subpopulations. Thus, region Ω is positively invariant. Hence, in this region, the model (3)–(10) is considered to be mathematically and epidemiologically well posed, and therefore, the dynamics of the model can be studied in Ω.

4. Existence and Stability of Equilibria

4.1. Disease-Free Equilibrium (DFE)

The DFE of the model is obtained by equating the right-hand sides of model equations (3)–(10) to zero as follows:

4.2. Calculation of Control Reproduction Number

The basic reproduction number in epidemic models is an important threshold value that quantifies the infection risk in order to effectively control the disease. Furthermore, it plays a vital role in stability analysis of equilibria of the models. It can be derived using the next-generation matrix approach [19]. However, when there is intervention, it is referred as the control reproduction number. For detailed computation of the control reproduction number, refer Appendix B. Therefore, the basic control reproduction number, denoted by ℛ, is given by the following equation:

It is worth stating that ℛ is the basic control reproduction number, which represents the number of secondary infection cases generated by introducing at least one infective agent into the population that is assumed wholly susceptible. This number is obtained from the contribution of average number of secondary infections through fox-to-environment-to-vole transmission , voles-to-fox transmission , and environment-to-voles transmission (ℛ=η/Kμ) as a result of one infectious subject during its infectious period.

4.2.1. Stability of DFE

Theorem 2 .

The disease-free equilibrium is locally asymptotically stable when ℛ < 1 and unstable if ℛ > 1.

The local asymptotic stability of DFE is established using Theorem 2 in [19].

The epidemiological implication of the result in Theorem 2 is that the Echinococcus m. can be eliminated from the populations when ℛ < 1 if the initial sizes of subpopulations are within the basin of attraction of the DFE. In order to guarantee the total elimination of the disease irrespective of the initial population started with in Ω, it is necessary to prove the global asymptotic stability (in Ω) of the DFE, which is shown below. Here, we use the matrix-theoretic method as described in [20].

Theorem 3 .

The DFE of the model (3)–(10) given by (12) is globally asymptotically stable (GAS) in Ω when ℛ < 1. If ℛ > 1, the DFE is unstable, the system is uniformly persistent, and there exists at least one endemic equilibrium in the interior of Ω.

Let F and V matrices be given as above. The matrix V−1F is computed to be

Observe that the matrix V−1F is reducible (second, fifth, and sixth columns are the nonzero columns); hence, Theorem 2.2 of [20] is not applicable, and instead we use the conditions of Theorem 2.1 in [20] to construct the Lyapunov function. Define X=(E, I, R, E, I, B) and Y=(S, S) so that the dynamics of infected compartments can be expressed as follows:wheresince S0 ≥ S and S0 ≥ S inside Ω.

We define the Lyapunov function as L=uV−1X, where u is the left eigenvector of the matrix V−1F with respect to the eigenvalue ℛ. Thus,

Multiplying and equating (17), we have

Clearly, it can be seen that u1=0,  u3=0, and u4=0, and we can deduce that

If we take u6=1, one solution of u is as follows:

The Lyapunov function L is now given as

Differentiating L along the solutions of infected compartments in (3)–(10) gives

Using (13) and (21), simplifying and rearranging, we have

From (23), (K/(K+B)) < 1 and (S/N) ≤ 1, and then (K/(K+B))(S/N) < 1 and (S/N) ≤ 1. Therefore, if ℛ < 1, it implies that (dL/dt) ≤ 0. Furthermore, (dL/dt)=0 implies either E=0, I=0, R=0, E=0, I=0, or B=0. Thus, the largest invariant set where (dL/dt)=0 is the singleton (ℰ0). Therefore, by LaSalle's invariance principle [21], ℰ0 is GAS in Ω if ℛ < 1. Furthermore, if ℛ > 1 in (23), the first term is positive, while the second and third terms will be zero in Ω when I=B=0; therefore, (dL/dt) > 0, and hence, ℰ0 is unstable. Using the argument in Theorem 2.2 of Shuai and van den Driessche [20], it can be shown that the instability of ℰ0 when ℛ > 1 implies that the system is uniformly persistent in Ω, thus implying the existence of at least one positive endemic equilibrium.

4.3. Existence and Global Stability of Endemic Equilibrium

The existence of endemic equilibrium follows from the argument in Theorem 3. In the presence of disease in the community, the endemic equilibrium ℰ1 is obtained by setting the right-hand sides of equations (3)–(10) to zero, and thus,where

Theorem 4 .

If ℛ > 1, then the unique endemic equilibrium ℰ1 of model (3)–(10) is globally asymptotically stable (GAS) in Ω.

To prove the uniqueness and global stability of ℰ1, we apply the method of graph-theoretic as described in Section 3 of [20]. Detailed proof of the theorem is also given in Appendix C.

5. Numerics: Elasticity Indices, Numerical Simulations, and Control Strategies

In this section, we use the parameter values in Table 1 with the aim of illustrating the theoretical results and quantifying the control measures for Echinococcus multilocularis.

Table 1

Description of model parameters, indicating baselines, ranges, and references.

Parameter	Description	Baseline value	Value range with time unit 1 day
μ v	Vole per-capita birth rate = death rate	1/(30 × 6)	5 − 7 months [3]
ξ f	Rate of treatment for red fox	0.01	Assumed
β v	Transmission rate from environment to voles	0.004	Assumed
s	Searching efficiency of red fox	0.00002	Assumed
μ f	Red fox per-capita birth rate = death rate	1/(365 × 3.5)	2 − 5 years [3]
p	Probability that an infected vole preyed on infects a red fox	0.07	[5]
α v	Incubation rate in voles	1/28	25 − 31 days [15]
α f	Incubation rate in red fox	1/75	60 − 90 days [15]
μ b	Rate of cleaning/disinfection of environment	1/31	Assumed
η f	Rate of red fox contribution of E. multilocularis to environment	0.02	Assumed
K	Concentration of E. multilocularis in the feaces that yields 50% chance of catching Echinococcus	52	Assumed

5.1. Elasticity Indices

As evident from the expression of basic control reproduction number, ℛ1 in (13), it is interesting to know qualitatively and estimate quantitatively how perturbations of associated parameters have influence on ℛ. In order to achieve this, we determine the normalized forward sensitivity index as introduced in Chitnis et al. [22], otherwise called elasticity indices [23] of parameters on ℛ. This quantity ϒ for ℛ with respect to a parameter p is defined as follows:

In Table 2, we compute the elasticity indices of ℛ with respect to the parameters at static baseline values as indicated in Table 1 and arranged in the descending order of magnitudes. The computations indicate equal influence of six parameters associated with incidences for transmission of the parasites (β, s, p, η), disinfection rate (μ), and concentration of Echinococcus in the environment (K). Most importantly, treatment of red foxes (ξ) has the second largest value followed by incubation rates of voles and foxes in that order.

Table 2

Elasticity indices of ℛ relative to model parameters and ranges.

Parameter	Formula ((∂ℛc/∂p)(p/ℛc))	Baseline value	Range	Elasticity index
β v	1/3	0.004	0.00040 − 0.008	0.33333
s	−(1/3)	0.00002	(2 × 10−6) − (4 × 10−5)	−0.33333
μ b	−(1/3)	1/31	(1/62) − (1/3.1)	−0.33333
η f	1/3	0.02	0.002 − 0.04	0.33333
p	−(1/3)	0.07	0.007, 0.14	−0.33333
K	−(1/3)	52	26 − 78	−0.33333
ξ f	−(1/3)(ξf/(μf+ξf))	0.01	0.001 − 0.02	−0.30913
α v	(1/3)(μv/(μv+αv))	1/28	(1/56) − (1/2.8)	0.044872
α f	(1/3)(μf/(μf+αf))	1/75	(1/75 × 2) − (1/7.5)	0.018484

5.2. Global Sensitivity Analysis

From our result in Section 5.1, it is obvious that the local sensitivity analysis on ℛ could not explicitly differentiate the most influential parameters and thus the need for global sensitivity analysis. We adapt the approach in [24] to analyze the global sensitivity of the parameters on both ℛ and cumulative new cases in rodents and red foxes, respectively. Using the method of partial rank correlation coefficients (PRCC), as described and implemented in [25], we carry out the global sensitivity analysis of 9 parameters on the control reproduction number ℛ and cumulative new infections in the populations of both red foxes and rodents. The main objective is to determine the most influential parameters for the purpose of control and extent of infectivity in the two populations. To compute the PRCC values, we used the MatLab R2017b with ranges of parameters in Table 2 divided into 1000 sample sizes, and the results are displayed in Figure 2. The parameter with the PRCC value far away from zero indicates the more influential parameter is on both ℛ and cumulative new cases.

Figure 2

Global sensitivity analysis displaying the partial rank correlation coefficients (PRCC) of (a) control reproduction number ℛ, (b) cumulative new cases in red foxes, and (c) cumulative new cases in voles populations.

In Figure 2(a), the global sensitivity of parameters on ℛ is depicted. It can be seen that the rates of cleaning/disinfecting the environment (μ) and rate of treating red foxes (ξ) have the most global influence on ℛ, followed by rate of red foxes contribution of E. multilocularis to the environment (η). The global sensitivity of parameters on the cumulative number of new cases for red foxes is also displayed in Figure 2(b) which indicates that the incubation rate in red foxes (α) has the highest global influence, followed by the rate of searching efficiency of red foxes (s) and probability that an infected vole preyed on infects a red fox (p) in that order. Lastly, in Figure 2(c), the global sensitivity of parameters on cumulative new infection cases in rodents indicates that the transmission rate from environment to rodents (β) is the most global influential parameter, followed by the incubation rate in rodents (α). From the global sensitivity analysis, for control purposes, it can be suggested that more emphasis should be given to cleaning/disinfecting the environment, for example, by removing carcass and administering praziquantel to red foxes.

5.3. Numerical Simulations

Figure 3 depicts the global stability of disease-free equilibrium as proved in Theorem 3 with different initial conditions, where the numbers of foxes, voles, and concentration of E. multilocularis converges asymptotically to the equilibrium point using different initial conditions. The parameter values in Table 1 are used so that the control reproduction number ℛ=0.97 < 1. It can be seen that all disease compartments (E, I, R, E, I,  and B) converge asymptotically to zero while the noninfected compartments (S and S) converge to their respective total populations.

Figure 3

Time evolution of (a) red foxes, (b) voles, and (c) concentration of E. multilocularis using different initial conditions when there is no infection with parameter values in Table 1, which gives the control reproduction number ℛ=0.97.

In Figure 4, the time evolution for number of red foxes, voles, and the concentration of Echinococcus multilocularis for model (3)–(10) is illustrated using the parameter values in Table 1, except for η=0.2 depicting GAS of endemic equilibrium as proved in Theorem 4. It can be observed that populations of foxes, voles, and concentration of parasite converge asymptotically to their respective endemic equilibrium points irrespective of the initial population started with. Here, the control reproduction number ℛ=1.94 > 1. Hence, the disease will persist in the community.

Figure 4

Time evolution of (a) red foxes, (b) voles, and (c) concentration of E. multilocularis with different initial conditions when there is infection using parameter values in Table 1, except for ξ=0.2 with control reproduction number ℛ=1.94.

5.4. Effects of Control Strategies on ℛ

In this section, the effects of treating red foxes (ξ) and the cleaning or disinfecting the environment (μ) are going to be explored. So here it will be interesting to see how values of infectious red foxes and/or the ℛ will change as these parameters are varied when other parameters are fixed at the baseline values.

The effect of treatment-only on red foxes is illustrated in Figure 5(a) using the baseline parameter values in Table 1 except for μ=0. The value of ξ is varied from ξ=0.01/102 to ξ=0.01; as a result, the number of infectious red foxes is reduced from 126 to 19, respectively, as displayed by the graphs. Similarly, the effect of controlling the red foxes by cleaning and/or disinfecting environment-only is shown in Figure 5(b) with fixed parameter values except ξ=0, and the value of μ varied from μ=(1/31)/102 to μ=(1/31). It can be seen that the cumulative number of infectious red foxes decrease from 134 with ℛ=10.01 to 8, ℛ=2.16, respectively. It is evident that the implementation of either of the two control strategies may not be adequate in eradicating the parasite completely from the community. Therefore, when the control strategies are administered simultaneously, as depicted in Figure 5(c), the cumulative number of infectious red foxes decreases from 125 with ℛ=9.60 to 0, with ℛ=0.80, respectively. Hence, combining the two control strategies is more effectively followed by environmental cleaning/disinfection and treatment of red foxes. The later results agree with our elasticity indices in Section 5.1 for the two parameters (ξ and μ).

Figure 5

Numerical simulations displaying effects of control strategies on the cumulative number of infectious red foxes, using parameter values in Table 1 with varying values of μ and ξ, except in (a) effect of treating red foxes only (μ=0), (b) effect of environmental cleaning only (ξ=0), and (c) effect of combining controls.

Given that ξ and μ are the control parameters in the model, it is important to see how ℛ varies as the two parameters are varied with others fixed using contour plots. The objective is to estimate values of ξ and μ that will ensure disease eradication (making ℛ < 1) as stated in Theorem 2. The results are displayed as contour curves of ℛ as a function of treatment on red foxes (ξ) and cleaning/disinfection of environment (μ) at fixed baseline values in Figure 6(a). The least values of ξ and μ that will ensure parasites eradication are estimated to be 0.025 and 0.01 so that ℛ=0.99 or 0.005 and 0.05 with ℛ=0.95, respectively. Furthermore, to access the impact of combined treatment and reduced contribution of parasite to environment by red foxes, contour plots of ℛ as function of the control strategies with varying rate of contribution by red foxes to the environment (η) are displayed in Figures 6(b)–6(d). The figures show remarkable increase in the associated control reproduction number with increase in rate of contribution of parasites to environment by red foxes. In Figure 6(b), low control strategies are needed if the rate of contribution (η=0.001) is very small to ensure almost total eradication of the parasites, with range of ℛ ∈ [0.17, 1.17] and mean = 0.67. In Figure 6(c), with high contribution rate (η=0.2), the control strategies must also be high to lower the value of ℛ ∈ [0.99, 6.84] with mean = 3.92. However, when the rate of contribution (η=0.1) is moderate, in Figure 6(b), the control strategies must be in reciprocal combinations (low treatment rate versus high disinfection rate and vice versa) putting the range of ℛ ∈ [0.79, 5.43] with mean = 3.11.

Figure 6

Contour curves of ℛ as a function of ξ and μ using three different rates of red foxes contribution to environment (η), with (a) parameter values in Table 1, in (b) η=0.001, (c) η=0.10, and (d) η=0.20.

6. Concluding Remarks

A new global deterministic model for the transmission of Echinococcus multilocularis in the population of red foxes and voles with environment as a source of infection is formulated and used to access the impact of control strategies on the disease dynamics. Moreover, sensitivity analysis is carried out to determine the parameters that have influence on the control reproduction number and cumulative new infectious cases of red foxes and rodents. We start by investigating the basic properties of the model to ascertain its worthiness mathematically and epidemiologically. The major findings of the study are outlined as follows:

The disease-free equilibrium of the model is obtained and used to systematically determine the basic control reproduction number (ℛ). Furthermore, using a matrix-theoretic method, the DFE is globally asymptotically stable whenever ℛ is less than unity. The implication of this result is that infection of the parasite can be control in the community if ℛ can be reduced and maintain below unity.

When the control reproduction number exceeds unity, a unique endemic equilibrium exists, and using a graph-theoretic method, it is shown to be globally asymptotically stable if ℛ is greater than unity. Hence, the infection will persist in the community when the control strategies fail to reduce ℛ below unity.

Elasticity indices of ℛ on key parameters are computed, and the results indicate equal influence of six parameters associated with incidence for transmission of the disease (β, s, p,  and η), cleaning/disinfection (μ), and concentration of parasite in the environment (K), followed by rate of treatment on foxes (ξ) and incubation rates of voles (α) and foxes (α) in that order.

Having noticed that the local sensitivity analysis on ℛ could not differentiate explicitly the most influential parameter(s) of the model, a global sensitivity using PRCC is conducted. From the simulations, the two control parameters: rate of cleaning/disinfecting the environment (μ) and rate of treating red foxes have the most global influence on ℛ, followed by rate of red foxes contribution of E. multilocularis to the environment (η). The global sensitivity of parameters on the cumulative number of new cases for red foxes indicates that the incubation rate in red foxes (α) has the highest global influence, followed by rate of searching efficiency of red foxes (s) and probability that an infected vole preyed on infects a red fox (p) in that order. On the contrary, the global sensitivity of parameters on cumulative new infection cases in rodents shows that the transmission rate from environment to rodents (β) is the most global influential parameter, followed by incubation rate in rodents (α).

Numerical simulations with baseline values and varying ξ and μ, have indicated that administering disinfection of environment only induce more positive control impact (leaving only 8 infected red foxes) compared to applying treatment only on red foxes (leaving about 19 infected red foxes). However, administering the two control strategies induce the most positive impact by treating all the infected red foxes.

Contour curves are used to estimate the least values of two control parameters that will ensure disease eradication, i.e., the value of ℛ to be reduced and kept below unity with baseline values of other parameters fixed. These values of ξ and μ are estimated at 0.025 and 0.01 so that ℛ=0.99 or 0.005 and 0.05 with ℛ=0.95, respectively. Moreover, contour plots of ℛ as function of the control strategies with the varying rate of contribution by red foxes to the environment (η) are computed. The results show remarkable increase in the associated control reproduction number with increase in rate of contribution of parasites to environment by red foxes.

It is worth remarking that amongst the seven most influential parameters globally (β, s, μ, η, ξ, p,  and α), five are associated directly with the red foxes, and this justifies our choice of treatment on the population. This also explains the high prevalence of the disease on red foxes as reported in literature.