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.


Introduction
Echinococcus multilocularis (EM) is a parasitic taeniid tapeworm and one of the six species of the Echinococcus genus. e disease is sylvatic and has indirect life cycle between two major hosts in a predator-prey interaction [1][2][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). e 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]. e 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]. e 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][9][10][11][12]. Furthermore, the disease has high prevalence in red foxes population (1%-76.7%) [3,5,[13][14][15][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. ey 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. e 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. is paper is organized as follows. e 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.

Model Formulation
e total population of red foxes, which is assumed constant (birth and death rates, μ f , are assumed equal) in the environment at time t, denoted by N * f (t) is divided into susceptible (S f (t)), exposed (E f (t)), infected (I f (t)), and recovered (R f (t)) subpopulations so that e susceptible population is increased by recruitment of foxes by birth or immigration at rate μ f N * f 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. e 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 α f and μ f , respectively.
e infected fox population is increased by progression of exposed foxes and decreased as a result of treatment and natural deaths at rates ξ f and μ f , respectively. e population of recovered is increased by the treated infected foxes and decreased by natural death at rates ξ f and μ f , 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, μ v , are assumed equal) in the environment at time t, denoted by N * v (t) is subdivided into susceptible (S v (t)), exposed (E v (t)), and infected (I v (t)) subpopulations so that e population of susceptible voles is increased by birth at rate μ v N * v and decreased by infection from the concentration of parasites in the environment at the rate β v and natural deaths at rate μ v , 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 η f and decreased by disinfection or cleaning of environment at rate μ b . 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 3. Basic Properties of the Model e 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.
e following region is positively invariant for the model (3): Proof. e detailed proof of eorem 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. us, 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 Ω.

Disease-Free Equilibrium (DFE).
e DFE of the model is obtained by equating the right-hand sides of model equations (3)-(10) to zero as follows:

Calculation of Control Reproduction Number.
e 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. erefore, the basic control reproduction number, denoted by R c , is given by the following equation: It is worth stating that R c is the basic control reproduction number, which represents the number of secondary infection cases generated by introducing at least one Computational and Mathematical Methods in Medicine infective agent into the population that is assumed wholly susceptible. is number is obtained from the contribution of average number of secondary infections through fox-to- as a result of one infectious subject during its infectious period.

Stability of DFE
e disease-free equilibrium is locally asymptotically stable when R c < 1 and unstable if R c > 1.

Proof.
e local asymptotic stability of DFE is established using eorem 2 in [19].
□ e epidemiological implication of the result in eorem 2 is that the Echinococcus m. can be eliminated from the populations when R c < 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 matrixtheoretic method as described in [20]. (12) is globally asymptotically stable (GAS) in Ω when R c < 1. If R c > 1, the DFE is unstable, the system is uniformly persistent, and there exists at least one endemic equilibrium in the interior of Ω.

Theorem 3. e DFE of the model (3)-(10) given by
Proof. Let F and V matrices be given as above. e matrix V −1 F is computed to be Observe that the matrix V −1 F is reducible (second, fifth, and sixth columns are the nonzero columns); hence, eorem 2.2 of [20] is not applicable, and instead we use the conditions of eorem 2.1 in [20] to construct the so that the dynamics of infected compartments can be expressed as follows: where Multiplying and equating (17), we have Clearly, it can be seen that u 1 � 0, u 3 � 0, and u 4 � 0, and we can deduce that If we take u 6 � 1, one solution of u is as follows: e 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 v /N * v ) ≤ 1, and then erefore, by LaSalle's invariance principle [21], E 0 is GAS in Ω if R c < 1. Furthermore, if R c > 1 in (23), the first term is positive, while the second and third terms will be zero in Ω when I v � B � 0; therefore, (dL/dt) > 0, and hence, E 0 is unstable. Using the argument in eorem 2.2 of Shuai and van den Driessche [20], it can be shown that the instability of E 0 when R c > 1 implies that the system is uniformly persistent in Ω, thus implying the existence of at least one positive endemic equilibrium.

Existence and Global Stability of Endemic Equilibrium.
e existence of endemic equilibrium follows from the argument in eorem 3. In the presence of disease in the community, the endemic equilibrium E 1 is obtained by setting the right-hand sides of equations (3)-(10) to zero, and thus, where Theorem 4. If R c > 1, then the unique endemic equilibrium E 1 of model (3)-(10) is globally asymptotically stable (GAS) in Ω.
Proof. To prove the uniqueness and global stability of E 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.

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.

Elasticity Indices.
As evident from the expression of basic control reproduction number, R c 1 in (13), it is interesting to know qualitatively and estimate quantitatively how perturbations of associated parameters have influence on R c . 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 R c . is quantity Υ for R c with respect to a parameter p is defined as follows:

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

Global Sensitivity
Analysis. From our result in Section 5.1, it is obvious that the local sensitivity analysis on R c 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 R c 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 R c and cumulative new infections in the populations of both red foxes and rodents. e 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. e parameter with the PRCC value far away from zero indicates the more influential parameter is on both R c and cumulative new cases.
In Figure 2(a), the global sensitivity of parameters on R c is depicted. It can be seen that the rates of cleaning/disinfecting the environment (μ b ) and rate of treating red foxes (ξ f ) have the most global influence on R c , followed by rate of red foxes contribution of E. multilocularis to the environment (η f ). e 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 (α f ) 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 (β v ) is the most global influential parameter, followed by the incubation rate in rodents (α v ). 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. Figure 3 depicts the global stability of disease-free equilibrium as proved in eorem 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. e parameter values in Table 1 are used so that the control reproduction number R c � 0.97 < 1. It can be seen that all disease compartments (E * f , I * f , R * f , E * v , I * v , and B * ) converge asymptotically to zero while the noninfected compartments (S * f and S * v ) converge to their respective total populations.

Numerical Simulations.
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 η f � 0.2 depicting GAS of endemic equilibrium as proved in eorem 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 R c � 1.94 > 1. Hence, the disease will persist in the community.

Effects of Control
Strategies on R c . In this section, the effects of treating red foxes (ξ f ) and the cleaning or disinfecting the environment (μ b ) are going to be explored. So here it will be interesting to see how values of infectious red foxes and/or the R c will change as these parameters are varied when other parameters are fixed at the baseline values. e effect of treatment-only on red foxes is illustrated in Figure 5(a) using the baseline parameter values in Table 1 except for μ b � 0.
e value of ξ f is varied from ξ f � 0.01/10 2 to ξ f � 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 environmentonly is shown in Figure 5(b) with fixed parameter values except ξ f � 0, and the value of μ b varied from μ b � (1/31)/10 2 to μ b � (1/31). It can be seen that the cumulative number of infectious red foxes decrease from 134 with R c � 10.01 to 8, R c � 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. erefore, when the control strategies are administered simultaneously, as depicted in Figure 5(c), the cumulative number of infectious red foxes decreases from 125 with R c � 9.60 to 0, with R c � 0.80, respectively. Hence, combining the two control strategies is more effectively followed by environmental cleaning/disinfection and treatment of red foxes. e later results agree with our elasticity indices in Section 5.1 for the two parameters (ξ f and μ b ).
Given that ξ f and μ b are the control parameters in the model, it is important to see how R c varies as the two parameters are varied with others fixed using contour plots.  e objective is to estimate values of ξ f and μ b that will ensure disease eradication (making R c < 1) as stated in eorem 2.
e results are displayed as contour curves of R c as a function of treatment on red foxes (ξ f ) and cleaning/disinfection of environment (μ b ) at fixed baseline values in Figure 6(a). e least values of ξ f and μ b that will ensure parasites eradication are estimated to be 0.025 and 0.01 so that R c � 0.99 or 0.005 and 0.05 with R c � 0.95, respectively. Furthermore, to access the impact of combined treatment and reduced contribution of parasite to environment by red foxes, contour plots of R c as function of the control strategies with varying rate of contribution by red foxes to the environment (η f ) are displayed in Figures 6(b)-6(d).
e 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 (η f � 0.001) is very small to ensure almost total eradication of the parasites, with range of R c ∈ [0.17, 1.17] and mean � 0.67. In Figure 6(c), with high contribution rate (η f � 0.2), the control strategies must also be high to lower the value of R c ∈ [0.99, 6.84] with mean � 3.92. However, when the rate of contribution (η f � 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 R c ∈ [0.79, 5.43] with mean � 3.11.

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. e major findings of the study are outlined as follows: (1) e disease-free equilibrium of the model is obtained and used to systematically determine the basic control reproduction number (R c ). Furthermore, using a matrix-theoretic method, the DFE is globally asymptotically stable whenever R c is less than unity. e implication of this result is that infection of the parasite can be control in the community if R c can be reduced and maintain below unity.
(2) 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 R c is greater than

Concentration of E. multilocularis
Concentration of E. multilocularis

Concentration of E. multilocularis
Concentration of E. multilocularis  Table 1  transmission of the disease (β v , s, p, and η f ), cleaning/disinfection (μ b ), and concentration of parasite in the environment (K), followed by rate of treatment on foxes (ξ f ) and incubation rates of voles (α v ) and foxes (α f ) in that order.
(4) Having noticed that the local sensitivity analysis on R c 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 (μ b ) and rate of treating red foxes have the most global influence on R c , followed by rate of red foxes contribution of E. multilocularis to the environment (η f ). e global sensitivity of parameters on the cumulative number of new cases for red foxes indicates that the incubation rate in red foxes (α f ) 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 (β v ) is the most global influential parameter, followed by incubation rate in rodents (α v ). It is worth remarking that amongst the seven most influential parameters globally (β v , s, μ b , η f , ξ f , p, and α v ), five are associated directly with the red foxes, and this justifies our choice of treatment on the population. is also explains the high prevalence of the disease on red foxes as reported in literature.

A. Proof of Theorem 1
A.1. Existence and Uniqueness of Solutions. Let X � (S f , E f , I f , R f , S v , E v , I v , B) T be a column vector in R 8 + so that Define f(X) � (f 1 (X), f 2 (X), . . . , f 8 (X)) T to be the vector valued function in R 8 + , where f 1 (X), f 2 (X), . . . , f 8 (X) are the righthand sides of model (3)-(10). e model system with initial conditions can therefore be expressed as follows: (A.1) erefore, using a standard theorem of the dynamical system [26], f is the Lipschitz continuous in X. Hence, there exist a unique solution of (3)-(10) for all times t > 0.

A.2. Positivity of Solutions.
Here, we prove that any solution with the nonnegative initial condition that start in R 8 + will remain there (i.e., it is positive at all times t > 0).
Consider equation (1): and using the variable separable method, we have From (7), so that using the same method as above, we get Hence, For the positivity of other variables, we proceed by the method of contradiction as follows: Suppose by contradiction that the conclusion is not true, then there exists a time t such that so that Since d − (5) � 1, c 5 a 58 � c 6 a 65 ⟹ c 6 � (a 65 /a 58 )c 5 . Given d + (6) � 1, c 7 a 76 � c 6 a 65 + c 6 a 68 ⟹ c 7 � ((a 65 (a 65 + a 68 ))/a 58 a 76 ).
From d + (1) � 1, c 7 a 17 � c 2 a 27 ⟹ c 2 � (a 17 /a 21 )c 7 � ((a 65 a 71 (a 65 + a 68 ))/a 21 a 58 a 76 ). (C.14) It can easily be verified that {E * } is the only invariant set in the interior of Ω that can satisfy L ′ � 0 and thus the uniqueness of E * . Using this Lyapunov function in combination with LaSalle's invariance principle [21] completes the proof that E * is unique and globally asymptotically stable in Ω if R c > 1.

Data Availability
e type of data used in our research are numerical and all the data used are obtained from published literature which are adequately cited therein. Where data are not available (not cited), we use reasonable estimate of the value concern.
Conflicts of Interest e authors declare that they have no conflicts of interest.