A Schistosomiasis Model with Praziquantel Resistance

Abstract

A compartmental model is established for schistosomiasis with praziquantel resistance. The model considers the impact of genetic resistance and drug treatment on the transmission of schistosomiasis. We calculate the basic reproductive number and discuss the existence and stability of disease-free equilibrium, boundary equilibrium, and coexistence equilibrium. Our analysis shows that regardless of whether drug treatment leads to the emergence of resistance, once the impact of genetic resistance is larger, the resistant strain will be dominant, which is detrimental to the control of schistosomiasis. In addition, once the proportion of human with drug-resistant strain produced by drug treatment is larger, the number of human and snails with resistant strain is larger. This is not a good result for drug treatment with praziquantel.

1. Introduction

Currently, treatment of human beings infected by schistosomiasis primarily focuses on chemotherapy with praziquantel (PZQ). PZQ appeared as a new schistosomicidal compound during the 1970s [1]. In recent years, PZQ has become the drug of choice in most endemic areas because of its efficacy, its ease of administration, its tolerable side-effects, and its cost [1]. Although the effectiveness of PZQ against schistosomiasis is well documented, the precise mode of action of the drug has not clearly defined [2]. It is reported that the chemotherapy of many helminth infections is complicated by the occurrence of drug resistance and drug tolerance (a natural resistance) to certain anthelmintics [2]. Not surprisingly, recent epidemiological evidence suggests the emergence of PZQ-resistant⧵tolerant schistosomes [1, 3, 4]. Resistance is defined as a genetically transmitted loss of sensitivity in a parasite population that was previously sensitive to a given drug [2]. Tolerance is an innate insusceptibility of a parasite to a drug, with the caveat that the parasite must not have been previously exposed to the drug [2]. The first report of possible PZQ resistance came from an intensive focus in northern Senegal, where the drug had produced very low cure rates (18–39%) [5, 6]. And snails collected in the area carried schistosoma strains. When tested in the laboratory, those snails had a decreased susceptibility to PZQ [7, 8]. Additional PZQ-resistant evidence was collected in Egypt [4, 9]. Preliminary studies have begun on these isolates to identify genetic, physiological, and morphological characteristics associated with PZQ resistance, and some of these may find use as markers for monitoring whether or not resistance is developing in endemic areas, where the drug is used [10].

Many papers have reported that drug treatment results in the emergence of schistosome resistance to PZQ [11–13]. Drug treatment can remove drug-susceptible parasites in infected human beings, while resistant parasites survive. However, many investigations find that traits of PZQ resistance of Schistosoma mansoni are dominant inheritance [14–16]. The resistant worms can reproduce and pass the resistant genes to the next generation. Furthermore, the resistance of Schistosoma mansoni to PZQ can be expressed in eggs, miracidia, cercariae, adults, and all stages of development [14–16]. In other words, some definitive hosts carrying resistant schistosomes can infect snails and make those snails carry resistant schistosomes. On the contrary, some snails carrying resistant schistosomes can infect definitive hosts and also make those definitive hosts carry resistant schistosomes. It goes on, the control of schistosomiasis will face enormous difficulties. Therefore, it is necessary to study the impact of this schistosome genetic resistance on the transmission of schistosomiasis.

In previous schistosomiasis models, resistant problems have been studied by considering that the resistance of schistosomiasis is due to drug treatment [11–13]. In [11, 13], the authors proposed a multistrain schistosome model including sensitive and resistant parasite strains. Their goal was to infer the impact of drug treatment on the maintenance of schistosome genetic diversity. In their assumptions, the drug-sensitive parasite strain had an additional per capita death rate, 𝜎, due to treatment. For a parasite strain that had developed drug resistance with a resistance level 𝜃 (𝜃>1), this treatment-related death rate was assumed to be reduced by the factor 𝜃 to 𝜎/𝜃. Their results implied that higher treatment rate could allow for coexistence between sensitive and resistant parasite strains. In [12], the authors formulated a deterministic model with multiple strains of schistosomes in order to explore the role of drug treatment in the maintenance of a polymorphism of parasite strains that differed in their resistance levels. And snails infected by parasite strains were divided into multistrain subclasses according to the different level 𝜃_𝑖 of parasite strains. Analysis of the model showed that the likelihood that resistant strains would increase in frequency depended on the interplay between their relative fitness, the cost of resistance, and the degree of selection pressure exerted by drug treatments.

Motivated by [11–13], we establish a new model considering hosts with sensitive and resistant strains in this paper. Our purpose is mainly to study the impact of drug treatment and genetic resistance on the transmission of schistosomiasis.

Our paper is organized as follows. In Section 2, we establish a mathematical model with praziquantel resistance and obtain basic reproductive number and existence of equilibria. And then the stability of disease-free equilibrium is obtained in Section 3. Section 4 devotes to stability analysis of boundary equilibria. In Section 5, stability analysis of endemic equilibrium is performed.

2. Mathematical Model

According to different level of parasite strains, we divide infected hosts into sensitive and resistant strains. Considering resistance and inheritance of resistance, new resistant strains are composed of two parts. We classify definitive and intermediate hosts as susceptible, sensitive and resistant in the following:(i)𝑋_𝑠(𝑡), the population of susceptible human; (ii)𝑋_𝑖(𝑡), the population of human infected with sensitive parasite strain; (iii)𝑋_𝑟(𝑡), the population of human infected with resistant parasite strain; (iv)𝑌_𝑠(𝑡), the population of susceptible snail host; (v)𝑌_𝑖(𝑡), the population of snail host carrying sensitive parasites; (vi)𝑌_𝑟(𝑡), the population of snail host carrying resistant parasites.

We follow some of the available models for schistosomiasis [17–19] and assume that the reproduction rate of hosts is constant, and we ignore the recovery class of host since the life span of infected is short in comparison to that of human [18, 20]. All parameters in the model are assumed to be nonnegative constants:(i)Λ₁, recruitment rate of human; (ii)𝜇₁, per capita natural death rate of human; (iii)𝛼₁, per capita disease-induced death rate of human; (iv)𝛽₁, per capita contact transmission rate from infected snails to susceptible human; (v)Λ₂, recruitment rate of snails; (vi)𝜇₂, per capita natural death rate of snails; (vii)𝛼₂, per capita disease-induced death rate of snails; (viii)𝛽₂, per capita contact transmission rate from infected human to susceptible snails; (ix)𝜃, the rate of treatment for infected human; (x)𝜌, the proportion of human with drug-resistant strain produced by treatment; (xi)𝑘, represents the impact of inheritance and the cost of resistance on transmission rate, we assume that 𝑘≤1 since the transmission rate is reduced due to resistance [12]; (xii)𝑘₁, represents the impact of resistance on disease-induced death rate of human; (xiii)𝑘₂, represents the impact of resistance on disease-induced death rate of snails.

Then, we have a model with the form: 𝑑𝑋_𝑠𝑑𝑡=Λ₁−𝛽₁𝑋_𝑠𝑌_𝑖−𝑘𝛽₁𝑋_𝑠𝑌_𝑟−𝜇₁𝑋_𝑠+1−𝜌𝜎𝑋_𝑖,𝑑𝑋_𝑖𝑑𝑡=𝛽₁𝑋_𝑠𝑌_𝑖−𝜇₁𝑋_𝑖−𝛼₁𝑋_𝑖−𝜎𝑋_𝑖,𝑑𝑋_𝑟𝑑𝑡=𝑘𝛽₁𝑋_𝑠𝑌_𝑟+𝜌𝜎𝑋_𝑖−𝜇₁𝑋_𝑟−𝑘₁𝛼₁𝑋_𝑟,𝑑𝑌_𝑠𝑑𝑡=Λ₂−𝛽₂𝑌_𝑠𝑋_𝑖−𝑘𝛽₂𝑌_𝑠𝑋_𝑟−𝜇₂𝑌_𝑠,𝑑𝑌_𝑖𝑑𝑡=𝛽₂𝑌_𝑠𝑋_𝑖−𝜇₂𝑌_𝑖−𝛼₂𝑌_𝑖,𝑑𝑌_𝑟𝑑𝑡=𝑘𝛽₂𝑌_𝑠𝑋_𝑟−𝜇₂𝑌_𝑟−𝑘₂𝛼₂𝑌_𝑟. (1) Using standard methods, it is easy to see that disease free equilibrium 𝐸₀=(Λ₁/𝜇₁,0,0,Λ₂/𝜇₂,0,0) always exists. Let 𝑋_𝑠0=Λ₁/𝜇₁, 𝑌_𝑠0=Λ₂/𝜇₂, 𝑚₁=𝜇₁+𝛼₁+𝜎, 𝑚₂=𝜇₁+𝑘₁𝛼₁, 𝑚₃=𝜇₂+𝛼₂, and 𝑚₄=𝜇₂+𝑘₂𝛼₂. According to the concept of next generation matrix [21] and the formula of the basic reproductive number for ODE compartmental models [22], if we let 𝐹=00𝛽₁𝑋_𝑠00000𝑘𝛽₁𝑋_𝑠0𝛽₂𝑌_𝑠00000𝑘𝛽₂𝑌_𝑠0,𝑚00𝑉=¹0000𝑚₂0000𝑚₃0000𝑚₄, (2) one can calculate that the eigenvalues of the next generation matrix 𝐹𝑉⁻¹ are given by 𝑅_0𝑖=𝛽₁𝛽₂𝑋_𝑠0𝑌_𝑠0𝑚₁𝑚₃,𝑅_0𝑟=𝑘²𝛽₁𝛽₂𝑋_𝑠0𝑌_𝑠0𝑚₂𝑚₄. (3) Then, it follows that the basic reproductive number for the system (1) is given by 𝑅₀=max𝑅_0𝑖,𝑅_0𝑟. (4)

To obtain other equilibria, we let the right-hand side of (1) equal to zero and obtain Λ₁−𝛽₁𝑋_𝑠𝑌_𝑖−𝑘𝛽₁𝑋_𝑠𝑌_𝑟−𝜇₁𝑋_𝑠+1−𝜌𝜎𝑋_𝑖𝛽=0,₁𝑋_𝑠𝑌_𝑖−𝜇₁𝑋_𝑖−𝛼₁𝑋_𝑖−𝜎𝑋_𝑖=0,𝑘𝛽₁𝑋_𝑠𝑌_𝑟+𝜌𝜎𝑋_𝑖−𝜇₁𝑋_𝑟−𝑘₁𝛼₁𝑋_𝑟Λ=0,₂−𝛽₂𝑌_𝑠𝑋_𝑖−𝑘𝛽₂𝑌_𝑠𝑋_𝑟−𝜇₂𝑌_𝑠𝛽=0,₂𝑌_𝑠𝑋_𝑖−𝜇₂𝑌_𝑖−𝛼₂𝑌_𝑖=0,𝑘𝛽₂𝑌_𝑠𝑋_𝑟−𝜇₂𝑌_𝑟−𝑘₂𝛼₂𝑌_𝑟=0. (5)(6)(7)(8)(9)(10) If 𝑋_𝑖=0, we have 𝑌_𝑖=0 and a formula of 𝑋_𝑟 as follows: 𝑘²𝛽₁𝛽₂Λ₂𝑚₂+𝜇₁𝑚₂𝑚₄𝑘𝛽₂𝑋_𝑟=𝑘²𝛽₁𝛽₂Λ₁Λ₂−𝜇₁𝜇₂𝑚₂𝑚₄. (11)

Note that 𝑅_0𝑟>1 equals to 𝑘²𝛽₁𝛽₂Λ₁Λ₂>𝜇₁𝜇₂𝑚₂𝑚₄. Hence, if 𝑅_0𝑟>1, there exists a boundary equilibrium with only resistant type, given as 𝐸_0𝑟=(𝑋_𝑠𝑟,0,𝑋_𝑟𝑟,𝑌_𝑠𝑟,0,𝑌_𝑟𝑟), where 𝑋_𝑠𝑟=Λ₁𝜇₁𝑚₄𝑅²_0𝑟+𝑘𝛽₁Λ₂𝜇²₁𝑚₄+𝑘𝛽₁𝜇₁Λ₂𝑅²_0𝑟,𝑋_𝑟𝑟=𝜇₁𝜇₂𝑚₄𝑅²_0𝑟−1𝜇₁𝑚₄+𝑘𝛽₁Λ₂𝑘𝛽₂,𝑌_𝑠𝑟=Λ₂𝜇₁𝑚₄+𝑘𝛽₁Λ₂𝜇₁𝜇₂𝑚₄𝑅²_0𝑟+𝑘𝛽₁Λ₂𝜇₂,𝑌_𝑟𝑟=Λ₂𝜇₁𝑅²_0𝑟−1𝜇₁𝑚₄𝑅²_0𝑟+𝑘𝛽₁Λ₂. (12)

If drug treatment does not lead to drug resistance, that is, 𝜌=0, then when 𝑅_0𝑖>1, we can obtain the other boundary equilibrium with only sensitive type, given as 𝐸_0𝑖=(𝑋_𝑠𝑖,𝑋_𝑖𝑖,0,𝑌_𝑠𝑖,𝑌_𝑖𝑖,0), where 𝑋_𝑠𝑖=Λ₁𝑚₁𝑚₃+𝜇₁+𝛼₁𝜇₂𝑚₁𝑚₃𝛽₁𝛽₂Λ₂𝜇₁+𝛼₁+𝜇₁𝑚₁𝑚₃,𝑋_𝑖𝑖=𝜇₁𝜇₂𝑚₁𝑚₃𝑅²_0𝑖−1𝛽₁𝛽₂Λ₂𝜇₁+𝛼₁+𝜇₁𝑚₁𝑚₃,𝑌_𝑠𝑖=𝛽₁𝛽₂Λ²₂𝜇₁+𝛼₁+Λ₂𝜇₁𝑚₁𝑚₃𝛽₂𝜇₁𝜇₂𝑚₁𝑚₃𝑅²_0𝑖−1+𝛽₁𝛽₂Λ₂𝜇₂𝜇₁+𝛼₁+𝜇₁𝜇₂𝑚₁𝑚₃,𝑌_𝑖𝑖=Λ₂𝛽₂𝜇₁𝜇₂𝑚₁𝑅²_0𝑖−1𝛽₂𝜇₁𝜇₂𝑚₁𝑚₃𝑅²_0𝑖−1+𝛽₁𝛽₂Λ₂𝜇₂𝜇₁+𝛼₁+𝜇₁𝜇₂𝑚₁𝑚₃. (13)

Now, we study existence of coexistence equilibrium for the system (1). From (5)–(10), we obtain 𝛽₁𝛽₂Λ₂𝜇₁+𝛼₁+𝜇₁𝑚₁𝑚₃𝛽₂𝑋_𝑖+𝜇₁𝑚₁𝑚₃𝑘𝛽₂+𝛽₁𝛽₂Λ₂𝑚₂𝑋_𝑟=𝛽₁𝛽₂Λ₁Λ₂−𝜇₁𝜇₂𝑚₁𝑚₃. (14) Following (6), (7), (9), and (10), we have 𝑌_𝑖𝑌_𝑟=𝑘𝑚₁𝑋_𝑖𝑚₂𝑋_𝑟−𝜌𝜎𝑋_𝑖=𝑚₄𝑋_𝑖𝑘𝑚₃𝑋_𝑟. (15)

In the case that 𝜌=0, (15) leads to 𝑘𝑚₁𝑋_𝑖𝑚₂𝑋_𝑟=𝑚₄𝑋_𝑖𝑘𝑚₃𝑋_𝑟. (16) Note that 𝑅_0𝑖=𝑅_0𝑟 equals to 𝑘𝑚₁/𝑚₂=𝑚₄/𝑘𝑚₃. It is easy to see that if 𝑅_0𝑖≠𝑅_0𝑟, (16) cannot hold, which implies that there is not coexistence equilibrium. If 𝑅_0𝑖=𝑅_0𝑟, (16) always holds. Hence, if 𝑅_0𝑟=𝑅_0𝑖>1, that is, 𝛽₁𝛽₂Λ₁Λ₂>𝜇₁𝜇₂𝑚₁𝑚₃, (14) represents the existence of coexistence equilibrium in the form of a line.

In the case that 𝜌≠0, (15) leads to 𝑚₄𝜌𝜎𝑋_𝑖=𝑚₂𝑚₄−𝑘²𝑚₁𝑚₃𝑋_𝑟. (17) Note that 𝑅_0𝑖>𝑅_0𝑟 equals to 𝑚₂𝑚₄−𝑘²𝑚₁𝑚₃>0. Following (14) and (17), we can obtain that if 𝑅_0𝑖>𝑅_0𝑟 and 𝑅_0𝑖>1, then the unique coexistence equilibrium is given by 𝐸_𝑐=(𝑋_𝑠𝑐,𝑋_𝑖𝑐,𝑋_𝑟𝑐,𝑌_𝑠𝑐,𝑌_𝑖𝑐,𝑌_𝑟𝑐), where 𝑋_𝑟𝑐=𝑚₄𝜌𝜎𝛽₁𝛽₂Λ₁Λ₂−𝜇₁𝜇₂𝑚₁𝑚₃×𝛽₁𝛽₂Λ₂𝜇₁+𝛼₁+𝜇₁𝑚₁𝑚₃𝛽₂×𝑚₂𝑚₄−𝑘²𝑚₁𝑚₃+𝑚₄𝜌𝜎𝜇₁𝑚₁𝑚₃𝑘𝛽₂+𝛽₁𝛽₂Λ₂𝑚₂⁻¹,𝑋_𝑖𝑐=𝑚₂𝑚₄−𝑘²𝑚₁𝑚₃𝑚₄𝑋𝜌𝜎^𝑟𝑐,𝑌_𝑠𝑐=Λ₂𝛽₂𝑋_𝑖𝑐+𝑘𝛽₂𝑋_𝑟𝑐+𝜇₂,𝑌_𝑟𝑐=𝑘Λ₂𝛽₂𝑋_𝑟𝑐𝑚₄𝛽₂𝑋_𝑖𝑐+𝑘𝛽₂𝑋_𝑟𝑐+𝜇₂,𝑌_𝑖𝑐=𝛽₂𝑌_𝑠𝑐𝑋_𝑖𝑐𝜇₂+𝛼₂,𝑋_𝑠𝑐=Λ₁+1−𝜌𝜎𝑋_𝑖𝑐𝛽₁𝑌_𝑖𝑐+𝑘𝛽₁𝑌_𝑟𝑐+𝜇₁. (18)

Summarizing above analyses, we have the following result.

Lemma 1. The existence of equilibria for the system (1) is as follows:(1)the disease free equilibrium 𝐸₀ always exists;(2)if 𝑅_0𝑟>1, there exists a boundary equilibrium with only resistant type 𝐸_0𝑟;(3)if 𝜌=0 and 𝑅_0𝑖>1, there exists a boundary equilibrium with only sensitive type 𝐸_0𝑖;(4)if 𝜌=0 and 𝑅_0𝑖=𝑅_0𝑟>1, there exists coexistence equilibrium in the form of a line (14);(5)if 𝜌≠0 and 𝑅_0𝑖>𝑅_0𝑟 and 𝑅_0𝑖>1, there exists a unique coexistence equilibrium 𝐸_𝑐.

The following section shows that the basic reproductive number 𝑅₀ provides a threshold condition for schistosoma extinction in (1).

3. Stability Analysis of the Disease Free Equilibrium

In this section, we will analyze stability of the disease free equilibrium of the model (1). The stability of the disease free equilibrium determines whether schistosomiasis will be permanent in an uninfected population. The following result shows that schistosome will go extinct if 𝑅₀<1.

Theorem 2. The disease free equilibrium 𝐸₀ of the system (1) is locally asymptotically stable if 𝑅₀<1 and unstable if 𝑅₀>1.

Proof. The Jacobian matrix for the system (1) is given by𝐽=−𝛽₁𝑌_𝑖−𝑘𝛽₁𝑌_𝑟−𝜇₁1−𝜌𝜎00−𝛽₁𝑋_𝑠−𝑘𝛽₁𝑋_𝑠𝛽₁𝑌_𝑖−𝑚₁00𝛽₁𝑋_𝑠0𝑘𝛽₁𝑌_𝑟𝜌𝜎−𝑚₂00𝑘𝛽₁𝑋_𝑠0−𝛽₂𝑌_𝑠−𝑘𝛽₂𝑌_𝑠−𝛽₂𝑋_𝑖−𝑘𝛽₂𝑋_𝑟−𝜇₂000𝛽₂𝑌_𝑠0𝛽₂𝑋_𝑖−𝑚₃000𝑘𝛽₂𝑌_𝑠𝑘𝛽₂𝑋_𝑟0−𝑚₄. (19)Then, the eigenvalues of 𝐸₀ are −𝜇₁, −𝜇₂ and roots of the following equations: 𝜆²+𝑚₁+𝑚₃𝜆+𝑚₁𝑚₃−𝛽₁𝛽₂𝑋_𝑠0𝑌_𝑠0𝜆=0,₂+𝑚₂+𝑚₄𝜆+𝑚₂𝑚₄−𝑘²𝛽₁𝛽₂𝑋_𝑠0𝑌_𝑠0=0. (20) Note that 𝑅₀<1 equals to 𝑅_0𝑖<1 and 𝑅_0𝑟<1, which leads to 𝑚₁𝑚₃−𝛽₁𝛽₂𝑋_𝑠0𝑌_𝑠0>0 and 𝑚₂𝑚₄−𝑘²𝛽₁𝛽₂𝑋_𝑠0𝑌_𝑠0>0. Hence, if 𝑅₀<1, then all the roots of (20) have negative real parts. Hence, using the Routh-Hurwitz criterion, we can obtain that the disease free equilibrium 𝐸₀ of the system (1) is locally asymptotically stable if 𝑅₀<1 and unstable if 𝑅₀>1.

Now, we turn to the study of the global stability of the disease free equilibrium of the model (1) by using Metzler matrix theory and the technique of Kamgang and Sallet [23].

Consider systems of the following form: 𝑑𝑥₁𝑥𝑑𝑡=𝑓₁,𝑥₂,𝑑𝑥₂𝑥𝑑𝑡=𝑔₁,𝑥₂, (21) where 𝑥₁∈𝑅^𝑛1₊, 𝑥₂∈𝑅^𝑛2₊, and 𝑓 and 𝑔 are 𝐶¹. We denote by 𝑥=(𝑥₁,𝑥₂) the state of the system and (𝑥^∗₁,0) is a disease free equilibrium on a positively invariant set Ω⊂𝑅^𝑛1+𝑛2₊. Now rewrite (21) as 𝑑𝑥₁𝑑𝑡=𝐴₁(𝑥)⋅𝑥₁−𝑥^∗₁+𝐴₁₂(𝑥)⋅𝑥₂,𝑑𝑥₂𝑑𝑡=𝐴₂(𝑥)𝑥₂. (22) For the system (22), we make the following assumptions.(h₁) The system is defined on the positively invariant set Ω of the nonnegative orthant. The system is dissipative on Ω.(h₂) The subsystem 𝑑𝑥₁/𝑑𝑡=𝐴₁(𝑥)⋅(𝑥₁−𝑥^∗₁) is globally asymptotically stable at the equilibrium 𝑥^∗₁ on the canonical projection of Ω on 𝑅^𝑛1₊.(h₃) The matrix 𝐴₂(𝑥) is Metzler and irreducible for any given 𝑥∈Ω.(h₄) There exists an maximum matrix 𝐴₂, then for any 𝑥∈Ω such that 𝐴₂=𝐴₂(𝑥), 𝑥∈𝑅^𝑛1₊×{0}.(h₅)𝛼(𝐴₂)≤0, that is, the greatest real part of eigenvalues of 𝐴₂ is nonnegative.For convenience, we state two lemmas due to Kamgang and Sallet [23].

Lemma 3. If the above hypotheses, h₁–h₅, are satisfied, then the disease free equilibrium is globally asymptotically stable in Ω.

Lemma 4. If the same notations and hypotheses in Lemma 3 hold and if, furthermore, we have 𝐴₂=𝐴₂(𝑥^∗₁,0), the disease free equilibrium is globally asymptotically stable if and only if 𝛼(𝐴₂)≤0.

Next, we discuss the global stability of the disease free equilibrium 𝐸₀ of the system (1) using the above two Lemmas.

From the system (1), we know 𝑑𝑋_𝑠+𝑋_𝑖+𝑋_𝑟𝑑𝑡=Λ₁−𝜇₁𝑋_𝑠+𝑋_𝑖+𝑋_𝑟−𝛼₁𝑋_𝑖−𝑘₂𝛼₁𝑋_𝑟≤Λ₁−𝜇₁𝑋_𝑠+𝑋_𝑖+𝑋_𝑟,𝑑𝑌_𝑠+𝑌_𝑖+𝑌_𝑟𝑑𝑡=Λ₂−𝜇₂𝑌_𝑠+𝑌_𝑖+𝑌_𝑟−𝛼₂𝑌_𝑖−𝑘₄𝛼₂𝑌_𝑟≤Λ₂−𝜇₂𝑌_𝑠+𝑌_𝑖+𝑌_𝑟. (23) This proves that the set Ω=𝑋_𝑠,𝑋_𝑖,𝑋_𝑟,𝑌_𝑠,𝑌_𝑖,𝑌_𝑟∈ℝ⁵₊∣𝑋_𝑠+𝑋_𝑖+𝑋_𝑟≤Λ₁𝜇₁,𝑌_𝑠+𝑌_𝑖+𝑌_𝑟≤Λ₂𝜇₂ (24) is a compact positively invariant absorbing set contained in the nonnegative orthant. Thus, the system (1) is dissipative on Ω because the trajectories of (1) are forward bounded. Now, we will study the system (1) on Ω.

We set for system (1) 𝑥₁=(𝑋_𝑠,𝑌_𝑠), 𝑥₂=(𝑋_𝑖,𝑋_𝑟,𝑌_𝑖,𝑌_𝑟), and 𝑥^∗₁=(Λ₁/𝜇₁,Λ₂/𝜇₂). As in [23], we express the subsystem as 𝑑𝑥₁/𝑑𝑡=𝐴₁(𝑥₁,0)(𝑥₁−𝑥^∗₁) and 𝑑𝑋_𝑠𝑑𝑡=Λ₁−𝜇₁𝑋_𝑠,𝑑𝑌_𝑠𝑑𝑡=Λ₂−𝜇₂𝑌_𝑠. (25) This is a linear system, and its unique equilibrium (𝐴₁/𝜇₁,𝐴₂/𝜇₂) (corresponding to the disease free equilibrium of (1)) is globally asymptotically stable, hence the assumptions (h₁) and (h₂) are satisfied.

The matrix 𝐴₂(𝑥) is given by 𝐴₂(𝑥)=−𝑚₁0𝛽₁𝑋_𝑠0𝜌𝜎−𝑚₂0𝑘𝛽₁𝑋_𝑠𝛽₂𝑌_𝑠0−𝑚₃00𝑘𝛽₂𝑌_𝑠0−𝑚₄. (26) As required by hypothesis h₃, for any 𝑥∈Ω, the matrix 𝐴₂(𝑥) is irreducible.

Now, let us check (h₄). There is a maximum which is uniquely realized in Ω if 𝑋_𝑠=Λ₁/𝜇₁ and 𝑌_𝑠=Λ₂/𝜇₂, which corresponds to the disease free equilibrium. This maximum matrix is then 𝐽₂, the subblock of the Jacobian matrix at the disease free equilibrium, corresponding to the matrix 𝐴₂(𝑥). The matrix 𝐽₂ is given by 𝐽₂=−𝑚₁0𝛽₁𝑋_𝑠00𝜌𝜎−𝑚₂0𝑘𝛽₁𝑋_𝑠0𝛽₂𝑌_𝑠00−𝑚₃00𝑘𝛽₂𝑌_𝑠00−𝑚₄. (27) Therefore, we are in the situation of Lemma 4, where the maximum is attained at the disease free equilibrium.

The hypothesis (h₅) requires that 𝛼(𝐽₂)≤0. Writing 𝐽₂ as a block matrix 𝐽₂=𝐴_2×2𝐵_2×2𝐶_2×2𝐷_2×2. Since 𝐴 is already a Metzler stable matrix, the condition 𝛼(𝐽₂)≤0 is equivalent to the condition 𝛼(𝐷−𝐶𝐴⁻¹𝐵)≤0 [23], where 𝐷−𝐶𝐴⁻¹𝛽𝐵=¹𝛽₂𝑋_𝑠0𝑌_𝑠0𝑚₁−𝑚₃0𝜌𝜎𝑘𝛽₁𝛽₂𝑋_𝑠0𝑌_𝑠0𝑚₁𝑚₂𝑘²𝛽₁𝛽₂𝑋_𝑠0𝑌_𝑠0𝑚₂−𝑚₄. (28) Then, 𝛼𝐷−𝐶𝐴⁻¹𝐵𝛽=max¹𝛽₂𝑋_𝑠0𝑌_𝑠0−𝑚₁𝑚₃𝑚₁,𝑘²𝛽₁𝛽₂𝑋_𝑠0𝑌_𝑠0−𝑚₂𝑚₄𝑚₂. (29) Hence, the condition 𝛼(𝐷−𝐶𝐴⁻¹𝐵)≤0 is equivalent to 𝑅₀≤1. We have seen that the hypotheses (h₁), (h₂), (h₃), (h₄), and (h₅) are satisfied. Then, by Lemma 4, we have the following result.

Theorem 5. The disease free equilibrium 𝐸₀ of the system (1) is globally asymptotically stable if 𝑅₀≤1.

4. Stability Analysis of the Boundary Equilibria

In this section, we turn to study stability of the two boundary equilibria. From Lemma 1, we know that if 𝑅_0𝑟>1, there exists a boundary equilibrium with only resistant type 𝐸_0𝑟. Through calculations, we can obtain the characteristic equation as following: 𝜆+𝑚₁𝜆+𝑚₃−𝛽₁𝛽₂𝑋_𝑠𝑟𝑌_𝑠𝑟×𝜆+𝑚₂𝜆+𝑚₄𝜆+𝜇₁+𝑘𝛽₁𝑌_𝑟𝑟×𝜆+𝜇₂+𝑘𝛽₂𝑋_𝑟𝑟−𝑘²𝛽₁𝛽₂𝑋_𝑠𝑟𝑌_𝑠𝑟𝜆+𝜇₁𝜆+𝜇₂=0. (30) Hence, the eigenvalue of 𝐸_0𝑟 are roots of the following equations: 𝜆²+𝑚₁+𝑚₃𝜆+𝑚₁𝑚₃−𝛽₁𝛽₂𝑋_𝑠𝑟𝑌_𝑠𝑟𝜆=0,₄+𝑎₁𝜆³+𝑎₂𝜆²+𝑎₃𝜆+𝑎₄=0. (31)(32) Here, 𝑎₁=𝑚₂+𝑚₄+Λ₁𝑋_𝑠𝑟+Λ₂𝑌_𝑠𝑟𝑎>0,₂=𝑚₂𝑚₄+𝑚₂+𝑚₄Λ₁𝑋_𝑠𝑟+Λ₂𝑌_𝑠𝑟+Λ₁𝑋_𝑠𝑟Λ₂𝑌_𝑠𝑟−𝑘²𝛽₁𝛽₂𝑋_𝑠𝑟𝑌_𝑠𝑟,𝑎₃=𝑚₂𝑚₄Λ₁𝑋_𝑠𝑟+Λ₂𝑌_𝑠𝑟+𝑚₂+𝑚₄Λ₁𝑋_𝑠𝑟Λ₂𝑌_𝑠𝑟−𝑘²𝛽₁𝛽₂𝑋_𝑠𝑟𝑌_𝑠𝑟𝜇₁+𝜇₂,𝑎₄=𝑚₂𝑚₄Λ₁𝑋_𝑠𝑟Λ₂𝑌_𝑠𝑟−𝑘²𝛽₁𝛽₂𝑋_𝑠𝑟𝑌_𝑠𝑟𝜇₁𝜇₂. (33)

From (12), we can obtain 𝑋_𝑠𝑟𝑌_𝑠𝑟=𝑋_𝑠0𝑌_𝑠0/𝑅²_0𝑟=𝑚₂𝑚₄/𝑘²𝛽₁𝛽₂. Then, 𝑚₁𝑚₃−𝛽₁𝛽₂𝑋_𝑠𝑟𝑌_𝑠𝑟=𝑚₁𝑚₃−𝑚₂𝑚₄𝑘². (34) It is easy to see that 𝑅²_0𝑖/𝑅²_0𝑟=(𝑚₂𝑚₄/𝑘²)/𝑚₁𝑚₃. Then, 𝑅_0𝑖<𝑅_0𝑟 equals to 𝑚₁𝑚₃>𝛽₁𝛽₂𝑋_𝑠𝑟𝑌_𝑠𝑟, which implies that the roots of (31) have negative real parts if 𝑅_0𝑖<𝑅_0𝑟.

Based on (Λ₁/𝑋_𝑠𝑟)(Λ₂/𝑌_𝑠𝑟)=𝜇₁𝜇₂𝑅²_0𝑟 and 𝑘²𝛽₁𝛽₂𝑋_𝑠𝑟𝑌_𝑠𝑟=𝑚₂𝑚₄, we can obtain 𝑎₂=𝑚₂+𝑚₄Λ₁𝑋_𝑠𝑟+Λ₂𝑌_𝑠𝑟+Λ₁𝑋_𝑠𝑟Λ₂𝑌_𝑠𝑟𝑎>0,₃=𝑚₂𝑚₄𝑘𝛽₁𝜇₁Λ₂𝑅²_0𝑟−1𝜇₁𝑚₄𝑅²_0𝑟+𝑘𝛽₁Λ₂+𝜇₁𝜇₂𝑚₄𝑅²_0𝑟−1𝜇₁𝑚₄+𝑘𝛽₁Λ₂+𝑚₂+𝑚₄𝜇₁𝜇₂𝑅²_0𝑟,𝑎₄=𝑚₂𝑚₄𝜇₁𝜇₂𝑅²_0𝑟−1. (35) Here, 𝑎₃>0 and 𝑎₄>0 if 𝑅_0𝑟>1.

If 𝑅_0𝑟>1, let 𝑘𝛽₁𝜇₁Λ₂(𝑅²_0𝑟−1)/(𝜇₁𝑚₄𝑅²_0𝑟+𝑘𝛽₁Λ₂)=𝑝₁>0 and 𝜇₁𝜇₂𝑚₄(𝑅²_0𝑟−1)/(𝜇₁𝑚₄+𝑘𝛽₁Λ₂)=𝑝₂>0, then Λ₁/𝑋_𝑠𝑟=𝑝₁+𝜇₁>𝜇₁ and Λ₂/𝑌_𝑠𝑟=𝑝₂+𝜇₂>𝜇₂. Note that 𝑚₂=𝜇₁+𝑘₁𝛼₁>𝜇₁ and 𝑚₄=𝜇₂+𝑘₂𝛼₂>𝜇₂, then 𝑚₂+𝑚₄>𝜇₁+𝜇₂ and 𝑚₂𝑚₄>𝜇₁𝜇₂. Hence, we can obtain 𝐻₁=𝑎₁𝐻>0,₂=𝑎₁𝑎₃1𝑎₂=𝑎₁𝑎₂−𝑎₃=𝑚²₂+𝑚²₄+𝑚₂𝑚₄+𝜇₁𝜇₂𝑅²_0𝑟Λ₁𝑋_𝑠𝑟+Λ₂𝑌_𝑠𝑟+𝑚₂+𝑚₄Λ₁𝑋_𝑠𝑟+Λ₂𝑌_𝑠𝑟²+𝑚₂𝑚₄𝜇₁+𝜇₂𝐻>0,₃=𝑎₁𝑎₃01𝑎₂00𝑎₁𝑎₃=𝑎₁𝑎₂𝑎₃−𝑎²₃−𝑎²₁𝑎₄=𝑚₂+𝑚₄²𝑚₂𝑚₄Λ₁𝑋_𝑠𝑟+Λ₂𝑌_𝑠𝑟²+𝑚₂+𝑚₄³Λ₁𝑋_𝑠𝑟+Λ₂𝑌_𝑠𝑟𝜇₁𝜇₂𝑅²_0𝑟+𝜇₁𝜇₂𝑅²_0𝑟𝑚₂+𝑚₄𝑚₂𝑚₄Λ₁𝑋_𝑠𝑟+Λ₂𝑌_𝑠𝑟+𝑚₂+𝑚₄𝑚₂𝑚₄Λ₁𝑋_𝑠𝑟+Λ₂𝑌_𝑠𝑟³+𝑚₂+𝑚₄²Λ₁𝑋_𝑠𝑟+Λ₂𝑌_s𝑟²𝜇₁𝜇₂𝑅²_0𝑟+𝑚₂+𝑚₄Λ₁𝑋_𝑠𝑟+Λ₂𝑌_𝑠𝑟𝜇²₁𝜇²₂𝑅⁴_0𝑟+2𝜇₁+𝜇₂𝑚²₂𝑚²₄Λ₁𝑋_𝑠𝑟+Λ₂𝑌_𝑠𝑟+𝑚₂+𝑚₄𝜇₁𝜇₂𝑅²_0𝑟𝑚₂𝑚₄𝜇₁+𝜇₂+𝑚₂+𝑚₄²𝜇₁𝜇₂𝑚₂𝑚₄+2𝑚₂+𝑚₄Λ₁𝑋_𝑠𝑟+Λ₂𝑌_𝑠𝑟𝜇₁𝜇₂𝑚₂𝑚₄+Λ₁𝑋_𝑠𝑟+Λ₂𝑌_𝑠𝑟²𝜇₁𝜇₂𝑚₂𝑚₄−𝑚₂+𝑚₄²𝑚₂𝑚₄Λ₁𝑋_𝑠𝑟+Λ₂𝑌_𝑠𝑟𝜇₁+𝜇₂−𝑚₂+𝑚₄𝑚₂𝑚₄Λ₁𝑋_𝑠𝑟+Λ₂𝑌_𝑠𝑟²𝜇₁+𝜇₂−𝜇₁𝜇₂𝑅²_0𝑟𝜇₁+𝜇₂𝑚₂𝑚₄Λ₁𝑋_𝑠𝑟+Λ₂𝑌_𝑠𝑟−Λ₁𝑋_𝑠𝑟+Λ₂𝑌_𝑠𝑟²𝑚²₂𝑚²₄−𝑚²₂𝑚²₄𝜇₁+𝜇₂²−2𝑚₂𝑚₄𝑚₂+𝑚₄Λ₁𝑋_𝑠𝑟+Λ₂𝑌_𝑠𝑟𝜇₁𝜇₂𝑅²_0𝑟−𝑚₂𝑚₄(𝑚₂+𝑚₄)²𝜇₁𝜇₂𝑅²_0𝑟>𝑚²₂+𝑚²₄²𝑚₂𝑚₄𝑝²₁+𝑝²₂+𝑝₁𝑝₂+𝑝₁𝜇₁+𝑝₂𝜇₂+𝑚²₂𝑚²₄𝑝²₁+𝑝²₂+2𝑝₁𝜇₁+2𝑝₂𝜇₂𝐻>0,₄=𝑎₁𝑎₃001𝑎₂𝑎₄00𝑎₁𝑎₃001𝑎₂𝑎₄=𝑎₄𝐻₃>0. (36) It follows from Routh-Hurwitz criterion that all roots of (32) have negative real parts if 𝑅_0𝑟>1. Summering above analyses, we have the following result.

Theorem 6. The boundary equilibrium 𝐸_0𝑟 of the system (1) is locally asymptotically stable if 𝑅_0𝑟>1 and 𝑅_0𝑖<𝑅_0𝑟.

Now, we study the global stability of the boundary equilibrium 𝐸_0𝑟. Consider the Lyapunov function 𝑉₁=𝑑₁𝑋_𝑖+𝑑₂𝑌_𝑖, (37) where 𝑑₁=(𝑚₃+𝛽₂(Λ₂/𝜇₂))/𝑚₁𝑚₃(𝑅_0𝑖+1) and 𝑑₂=(𝑚₁+𝛽₁(Λ₁/𝜇₁))/𝑚₁𝑚₃(𝑅_0𝑖+1). The Lyapunov derivative is 𝑑𝑉₁𝑑𝑡=𝑑₁𝛽₁𝑋_𝑠𝑌_𝑖−𝑚₁𝑋_𝑖+𝑑₂𝛽₂𝑌_𝑠𝑋_𝑖−𝑚₃𝑌_𝑖≤𝑑₁𝛽₁Λ₁𝜇₁−𝑑₂𝑚₃𝑌_𝑖+𝑑₂𝛽₂Λ₂𝜇₂−𝑑₁𝑚₁𝑋_𝑖≤𝑅_0𝑖−1𝑌_𝑖+𝑋_𝑖≤0for𝑅_0𝑖≤1. (38) This implies that the sensitive type dies out if 𝑅_0𝑖≤1. Then the largest compact invariant set of the system (1) in the set {𝑑𝑉₁/𝑑𝑡=0} is Ω₁=(𝑋_𝑠,𝑋_𝑖,𝑋_𝑟,𝑌_𝑠,𝑌_𝑖,𝑌_𝑟)∈Ω∣𝑋_𝑖=𝑌_𝑖=0. Using the LaSalle-Lyapunov theorem, we know that all trajectories in Ω eventually tend to Ω₁ as 𝑡→∞. Then, we only need to study the dynamical behavior of (1) in Ω₁. At this time, (1) reduces to the following system 𝑑𝑋_𝑠𝑑𝑡=Λ₁−𝑘𝛽₁𝑋_𝑠𝑌_𝑟−𝜇₁𝑋_𝑠,𝑑𝑋_𝑟𝑑𝑡=𝑘𝛽₁𝑋_𝑠𝑌_𝑟−𝑚₂𝑋_𝑟,𝑑𝑌_𝑠𝑑𝑡=Λ₂−𝑘𝛽₂𝑌_𝑠𝑋_𝑟−𝜇₂𝑌_𝑠,𝑑𝑌_𝑟𝑑𝑡=𝑘𝛽₂𝑌_𝑠𝑋_𝑟−𝑚₄𝑌_𝑟. (39) To show that all trajectories of (39) in the interior of Ω₁ approach the point (𝑋_𝑠𝑟,𝑋_𝑟𝑟,𝑌_𝑠𝑟,𝑌_𝑟𝑟) corresponding to the boundary equilibrium 𝐸_0𝑟, consider the Lyapunov function 𝑉₂=𝑑₁𝑋_𝑠−𝑋_𝑠𝑟−𝑋_𝑠𝑟𝑋ln^𝑠𝑋_𝑠𝑟+𝑋_𝑟−𝑋_𝑟𝑟−𝑋_𝑟𝑟𝑋ln^𝑟𝑋_𝑟𝑟+𝑑₂𝑌_𝑠−𝑌_𝑠𝑟−𝑌_𝑠𝑟𝑌ln^𝑠𝑌_𝑠𝑟+𝑌_𝑟−𝑌_𝑟𝑟−𝑌_𝑟𝑟𝑌ln^𝑟𝑌_𝑟𝑟, (40) where positive constants 𝑑₁ and 𝑑₂ are defined in the following. It is easy to see that 𝑉₂≥0 for (𝑋_𝑠,𝑋_𝑟,𝑌_𝑠,𝑌_𝑟)∈Ω₁, and 𝑉₂=0⇔(𝑋_𝑠,𝑋_𝑟,𝑌_𝑠,𝑌_𝑟)=(𝑋_𝑠𝑟,𝑋_𝑟𝑟,𝑌_𝑠𝑟,𝑌_𝑟𝑟). Hence, the function 𝑉₂ is positive definite with respect to the point (𝑋_𝑠𝑟,𝑋_𝑟𝑟,𝑌_𝑠𝑟,𝑌_𝑟𝑟).

Computing the derivative of 𝑉₂ along solutions of system (39), we have 𝑑𝑉₂=𝑑𝑡𝑑₁𝑑𝑋_𝑠−𝑑𝑡𝑑𝑋_𝑠𝑋𝑑𝑡𝑋^𝑠𝑟_𝑠+𝑑𝑋_𝑟−𝑑𝑡𝑑𝑋_𝑟𝑋𝑑𝑡𝑋^𝑟𝑟_𝑟+𝑑₂𝑑𝑌_𝑠−𝑑𝑡𝑑𝑌_𝑠𝑌𝑑𝑡𝑌^𝑠𝑟_𝑠+𝑑𝑌_𝑟−𝑑𝑡𝑑𝑌_𝑟𝑌𝑑𝑡𝑌^𝑟𝑟_𝑟=𝑑₁Λ₁−𝜇₁𝑋_𝑠−𝑋_𝑠𝑟𝑋_𝑠Λ₁+𝑘𝛽₁𝑋_𝑠𝑟𝑌_𝑟+𝜇₁𝑋_𝑠𝑟−𝑚₂𝑋_𝑟−𝑘𝛽₁𝑋_𝑠𝑌_𝑟𝑋_𝑟𝑟𝑋_𝑟+𝑚₂𝑋_𝑟𝑟+𝑑₂Λ₂−𝜇₂𝑌_𝑠−𝑌_𝑠𝑟𝑌_𝑠Λ₂+𝑘𝛽₂𝑌_𝑠𝑟𝑋_𝑟+𝜇₂𝑌_𝑠𝑟−𝑚₄𝑌_𝑟−𝑘𝛽₂𝑌_𝑠𝑋_𝑟𝑌_𝑟𝑟𝑌_𝑟+𝑚₄𝑌_𝑟𝑟. (41) Substituting Λ₁=𝑘𝛽₁𝑋_𝑠𝑟𝑌_𝑟𝑟+𝜇₁𝑋_𝑠𝑟, Λ₂=𝑘𝛽₂𝑌_𝑠𝑟𝑋_𝑟𝑟+𝜇₂𝑌_𝑠𝑟, 𝑚₂𝑋_𝑟𝑟=𝑘𝛽₁𝑋_𝑠𝑟𝑌_𝑟𝑟 and 𝑚₄𝑌_𝑟𝑟=𝑘𝛽₂𝑌_𝑠𝑟𝑋_𝑟𝑟 into 𝑑𝑉₂/𝑑𝑡, we obtain 𝑑𝑉₂=𝑑𝑡𝑑₁𝜇₁𝑋_𝑠𝑟𝑋2−^𝑠𝑋_𝑠𝑟−𝑋_𝑠𝑟𝑋_𝑠+𝑘𝛽₁𝑋_𝑠𝑟𝑌_𝑟−𝑚₂𝑋_𝑟+2𝑘𝛽₁𝑋_𝑠𝑟𝑌_𝑟𝑟−𝑘𝛽₁𝑌_𝑟𝑟𝑋_𝑠𝑟²𝑋_𝑠−𝑘𝛽₁𝑋_𝑠𝑌_𝑟𝑋_𝑟𝑟𝑋_𝑟+𝑑₂𝜇₂𝑌_𝑠𝑟𝑌2−^𝑠𝑌_𝑠𝑟−𝑌_𝑠𝑟𝑌_𝑠+𝑘𝛽₂𝑌_𝑠𝑟𝑋_𝑟−𝑚₄𝑌_𝑟+2𝑘𝛽₂𝑌_𝑠𝑟𝑋_𝑟𝑟−𝑘𝛽₂𝑋_𝑟𝑟𝑌_𝑠𝑟²𝑌_𝑠−𝑘𝛽₂𝑌_𝑠𝑋_𝑟𝑌_𝑟𝑟𝑌_𝑟. (42) Let 𝑑₁=𝑘𝛽₂𝑌_𝑠𝑟𝑋_𝑟𝑟 and 𝑑₂=𝑘𝛽₁𝑋_𝑠𝑟𝑌_𝑟𝑟. Following from 𝑚₂𝑋_𝑟𝑟=𝑘𝛽₁𝑋_𝑠𝑟𝑌_𝑟𝑟 and 𝑚₄𝑌_𝑟𝑟=𝑘𝛽₂𝑌_𝑠r𝑋_𝑟𝑟, we have 𝑑₁𝑘𝛽₁𝑋_𝑠𝑟𝑌_𝑟−𝑚₂𝑋_𝑟+𝑑₂𝑘𝛽₂𝑌_𝑠𝑟𝑋_𝑟−𝑚₄𝑌_𝑟=0. (43) Hence, for all (𝑋_𝑠,𝑋_𝑟,𝑌_𝑠,𝑌_𝑟)∈Ω₁, 𝑑𝑉₂=𝑑𝑡𝑑₁𝜇₁𝑋_𝑠𝑟𝑋2−^𝑠𝑋_𝑠𝑟−𝑋_𝑠𝑟𝑋_𝑠+𝑘𝛽₁𝑋_𝑠𝑟𝑌_𝑟𝑟𝑋2−^𝑠𝑟𝑋_𝑠−𝑋_𝑠𝑌_𝑟𝑋_𝑟𝑟𝑋_𝑠𝑟𝑌_𝑟𝑟𝑋_𝑟+𝑑₂𝜇₂𝑌_𝑠𝑟𝑌2−^𝑠𝑌_𝑠𝑟−𝑌_𝑠𝑟𝑌_𝑠+𝑘𝛽₂𝑌_𝑠𝑟𝑋_𝑟𝑟𝑌2−^𝑠𝑟𝑌_𝑠−𝑌_𝑠𝑋_𝑟𝑌_𝑟𝑟𝑌_𝑠𝑟𝑋_𝑟𝑟𝑌_𝑟=𝑑₁𝜇₁𝑋_𝑠𝑟𝑋2−^𝑠𝑋_𝑠𝑟−𝑋_𝑠𝑟𝑋_𝑠+𝑑₂𝜇₂𝑌_𝑠𝑟𝑌2−^𝑠𝑌_𝑠𝑟−𝑌_𝑠𝑟𝑌_𝑠+𝑑₁𝑑₂𝑋4−^𝑠𝑟𝑋_𝑠−𝑋_𝑠𝑌_𝑟𝑋_𝑟𝑟𝑋_𝑠𝑟𝑌_𝑟𝑟𝑋_𝑟−𝑌_𝑠𝑟𝑌_𝑠−𝑌_𝑠𝑋_𝑟𝑌_𝑟𝑟𝑌_𝑠𝑟𝑋_𝑟𝑟𝑌_𝑟≤0. (44) It is easy to see that 𝑑𝑉₂𝑑𝑡=0⟺𝑋_𝑠=𝑋_𝑠𝑟,𝑋_𝑟=𝑎𝑋_𝑟𝑟,𝑌_𝑠=𝑌_𝑠𝑟,𝑌_𝑟=𝑎𝑌_𝑟𝑟, (45) where 𝑎 is an arbitrary positive number. Substituting 𝑋_𝑠=𝑋_𝑠𝑟, 𝑋_𝑟=𝑎𝑋_𝑟𝑟, 𝑌_𝑠=𝑌_𝑠𝑟, and 𝑌_𝑟=𝑎𝑌_𝑟𝑟 into the first equation of system (39), we obtain 0=Λ₁−𝑎𝑘𝛽₁𝑋_𝑠𝑟𝑌_𝑟𝑟−𝜇₁𝑋_𝑠𝑟, (46) and then the above formula holds if and only if 𝑎=1. Therefore, the only compact invariant subset of the set where 𝑑𝑉₂/𝑑𝑡=0 is the point (𝑋_𝑠𝑟,𝑋_𝑟𝑟,𝑌_𝑠𝑟,𝑌_𝑟𝑟), corresponding to the boundary equilibrium 𝐸_0𝑟. By LaSalle’s Invariance Principle, 𝐸_0𝑟 is globally asymptotically stable if 𝑅_0𝑟>1 and 𝑅_0𝑖≤1. Summering above analyses, we have the following result.

Theorem 7. The boundary equilibrium 𝐸_0𝑟 of the system (1) is globally asymptotically stable if 𝑅_0𝑟>1≥𝑅_0𝑖.

Now, we turn to the other boundary equilibrium. From Lemma 1, we know that if 𝜌=0 and 𝑅_0𝑖>1, the boundary equilibrium with only sensitive type 𝐸_0𝑖 exists. Through calculations, we can obtain the characteristic equation as following: 𝜆+𝑚₂𝜆+𝑚₄−𝑘²𝛽₁𝛽₂𝑋_𝑠𝑖𝑌_𝑠𝑖×𝜆+𝑚₁𝜆+𝑚₃𝜆+𝜇₁𝜆+𝜇₂+𝛽₂𝑋_𝑖𝑖+𝛽₁𝑌_𝑖𝑖𝜆+𝑚₃𝜆+𝜇₂+𝛽₂𝑋_𝑖𝑖𝜆+𝑚₁−𝜎−𝛽₁𝛽₂𝑋_𝑠𝑖𝑌_𝑠𝑖𝜆+𝜇₁𝜆+𝜇₂=0. (47) Hence, the eigenvalue of 𝐸_0𝑖 are roots of the following equations: 𝜆²+𝑚₂+𝑚₄𝜆+𝑚₂𝑚₄−𝑘²𝛽₁𝛽₂𝑋_𝑠𝑖𝑌_𝑠𝑖𝜆=0,₄+𝑏₁𝜆³+𝑏₂𝜆²+𝑏₃𝜆+𝑏₄=0. (48)(49) Here, 𝑏₁=𝑚₁+𝑚₃+𝜇₁+𝜇₂+𝛽₂𝑋_𝑖𝑖+𝛽₁𝑌_𝑖𝑖𝑏>0,₂=𝜇₁𝑚₁+𝑚₃+𝜇₂+𝛽₂𝑋_𝑖𝑖+𝑚₁𝑚₃+𝜇₂+𝛽₂𝑋_𝑖𝑖+𝑚₃𝜇₂+𝛽₂𝑋_𝑖𝑖+𝛽₁𝑌_𝑖𝑖𝑚₁−𝜎+𝑚₃+𝜇₂+𝛽₂𝑋_𝑖𝑖−𝑚₁𝑚₃=𝜇₁𝑚₁+𝑚₃+𝜇₂+𝛽₂𝑋_𝑖𝑖+𝑚₁𝜇₂+𝛽₂𝑋_𝑖𝑖+𝑚₃𝜇₂+𝛽₂𝑋_𝑖𝑖+𝛽₁𝑌_𝑖𝑖𝜇₁+𝛼+𝑚₃+𝜇₂+𝛽₂𝑋_𝑖𝑖𝑏>0,₃=𝜇₁𝑚₁𝑚₃+𝜇₂+𝛽₂𝑋_𝑖𝑖+𝑚₁𝑚₃𝜇₂+𝛽₂𝑋_𝑖𝑖+𝛽₁𝑌_𝑖𝑖𝑚₃𝜇₂+𝛽₂𝑋_𝑖𝑖+𝑚₁−𝜎𝑚₃+𝜇₂+𝛽₂𝑋_𝑖𝑖−𝑚₁𝑚₃𝜇₁+𝜇₂=𝜇₁𝑚₁𝜇₂+𝛽₂𝑋_𝑖𝑖+𝑚₁𝑚₃𝛽₂𝑋_𝑖𝑖+𝛽₁𝑌_𝑖𝑖𝑚₃𝜇₂+𝛽₂𝑋_𝑖𝑖+𝜇₁+𝛼𝑚₃+𝜇₂+𝛽₂𝑋_𝑖𝑖𝑏>0,₄=𝜇₁𝑚₁𝑚₃𝜇₂+𝛽₂𝑋_𝑖𝑖+𝛽₁𝑌_𝑖𝑖𝑚₃𝜇₂+𝛽₂𝑋_𝑖𝑖𝑚₁−𝜎−𝑚₁𝑚₃𝜇₁𝜇₂=𝜇₁𝑚₁𝑚₃𝛽₂𝑋_𝑖𝑖+𝛽₁𝑌_𝑖𝑖𝑚₃𝜇₂+𝛽₂𝑋_𝑖𝑖𝜇₁+𝛼>0. (50)

From (13), we can obtain that 𝛽₁𝛽₂𝑋_𝑠𝑖𝑌_𝑠𝑖=𝑚₂𝑚₄. Then, it is easy to see that the roots of (48) have negative real parts if 𝑅_0𝑖>𝑅_0𝑟. Similarly to the case of 𝐸_0𝑟, using Routh-Hurtwitz criterion, we can obtain that all roots of (49) have negative real parts if 𝑅_0𝑖>1. Summering above analysis, we have the following result.

Theorem 8. When 𝜌=0, the boundary equilibrium 𝐸_0𝑖 of the system (1) is locally asymptotically stable if 𝑅_0𝑖>1 and 𝑅_0𝑖>𝑅_0𝑟.

Now, we study the global stability of the boundary equilibrium 𝐸_0𝑖. Consider the Lyapunov function 𝐿₁=𝑙₁𝑋_𝑟+𝑙₂𝑌_𝑟, (51) where 𝑙₁=(𝑚₄+𝑘𝛽₂(Λ₂/𝜇₂))/𝑚₂𝑚₄(𝑅_0𝑟+1) and 𝑙₂=(𝑚₂+𝑘𝛽₁(Λ₁/𝜇₁))/𝑚₂𝑚₄(𝑅_0𝑟+1). The Lyapunov derivative is 𝑑𝐿₁𝑑𝑡=𝑙₁𝑘𝛽₁𝑋_𝑠𝑌_𝑟−𝑚₂𝑋_𝑟+𝑙₂𝑘𝛽₂𝑌_𝑠𝑋_𝑟−𝑚₃Y_𝑟≤𝑙₁𝑘𝛽₁Λ₁𝜇₁−𝑙₂𝑚₄𝑌_𝑟+𝑙₂𝑘𝛽₂Λ₂𝜇₂−𝑙₁𝑚₂𝑋_𝑟≤𝑅_0𝑟−1𝑌_𝑟+𝑋_𝑟≤0for𝑅_0𝑟≤1. (52) This implies that the resistant type dies out if 𝑅_0𝑟≤1. Then, the largest compact invariant set of the system (1) in the set {𝑑𝐿₁/𝑑𝑡=0} is Ω₂={(𝑋_𝑠,𝑋_𝑖,𝑋_𝑟,𝑌_𝑠,𝑌_𝑖,𝑌_𝑟)∈Ω∣𝑋_𝑟=𝑌_𝑟=0}. Using the LaSalle-Lyapunov theorem, we know that all trajectories in Ω eventually tend to Ω₂ as 𝑡→∞. Then, we only need to study the dynamical behavior of (1) in Ω₂. At this time, (1) reduces to the following system: 𝑑𝑋_𝑠𝑑𝑡=Λ₁−𝛽₁𝑋_𝑠𝑌_𝑖−𝜇₁𝑋_𝑠+𝜎𝑋_𝑖,𝑑𝑋_𝑖𝑑𝑡=𝛽₁𝑋_𝑠𝑌_𝑖−𝑚₁𝑋_𝑖,𝑑𝑌_𝑠𝑑𝑡=Λ₂−𝛽₂𝑌_𝑠𝑋_𝑖−𝜇₂𝑌_𝑠,𝑑𝑌_𝑖𝑑𝑡=𝛽₂𝑌_𝑠𝑋_𝑖−𝑚₃𝑌_𝑖. (53) To show that all trajectories of (53) in the interior of Ω₂ approach the point (𝑋_𝑠𝑖,𝑋_𝑖𝑖,𝑌_𝑠𝑖,𝑌_𝑖𝑖) corresponding to the boundary equilibrium 𝐸_0𝑖, consider the Lyapunov function 𝐿₂=𝑙₁𝑋_𝑠−𝑋_𝑠𝑖−𝑋_𝑠𝑖𝑋ln^𝑠𝑋_𝑠𝑖+𝑋_𝑖−𝑋_𝑖𝑖−𝑋_𝑖𝑖𝑋ln^𝑖𝑋_𝑖𝑖+𝑙₂𝑌_𝑠−𝑌_𝑠𝑖−𝑌_𝑠𝑖𝑌ln^𝑠𝑌_𝑠𝑖+𝑌_𝑖−𝑌_𝑖𝑖−𝑌_𝑖𝑖𝑌ln^𝑖𝑌_𝑖𝑖, (54) where positive constants 𝑙₁ and 𝑙₂ are defined in the following. It is easy to see that 𝐿₂≥0 for (𝑋_𝑠,𝑋_𝑖,𝑌_𝑠,𝑌_𝑖)∈Ω₁ and 𝐿₂=0⇔(𝑋_𝑠,𝑋_𝑖,𝑌_𝑠,𝑌_𝑖)=(𝑋_𝑠𝑖,𝑋_𝑖𝑖,𝑌_𝑠𝑖,𝑌_𝑖𝑖). Hence, the function 𝐿₂ is positive definite with respect to the point (𝑋_𝑠𝑖,𝑋_𝑖𝑖,𝑌_𝑠𝑖,𝑌_𝑖𝑖).

Computing the derivative of 𝐿₂ along solutions of system (39), we have 𝑑𝐿₂=𝑑𝑡𝑙₁𝑑𝑋_𝑠−𝑑𝑡𝑑𝑋_𝑠𝑋𝑑𝑡𝑋^𝑠𝑖_𝑠+𝑑𝑋_𝑖−𝑑𝑡𝑑𝑋_𝑖𝑋𝑑𝑡𝑋^𝑖𝑖_𝑖+𝑙₂𝑑𝑌_𝑠−𝑑𝑡𝑑𝑌_𝑠𝑌𝑑𝑡𝑌^𝑠𝑖_𝑠+𝑑𝑌_𝑖−𝑑𝑡𝑑𝑌_𝑖𝑌𝑑𝑡𝑌^𝑖𝑖_𝑖=𝑙₁Λ₁−𝜇₁𝑋_𝑠−𝑋_𝑠𝑖𝑋_𝑠Λ₁+𝛽₁𝑋_𝑠𝑖𝑌_𝑖+𝜇₁𝑋_𝑠𝑖−𝑚₁𝑋_𝑖−𝛽₁𝑋_𝑠𝑌_𝑖𝑋_𝑖𝑖𝑋_𝑟+𝑚₁𝑋_𝑖𝑖+𝜎𝑋_𝑖𝑋1−^𝑠𝑖𝑋_𝑠+𝑙₂Λ₂−𝜇₂𝑌_𝑠−𝑌_𝑠𝑖𝑌_𝑠Λ₂+𝛽₂𝑌_𝑠𝑖𝑋_𝑖+𝜇₂𝑌_𝑠𝑖−𝑚₃𝑌_𝑖−𝛽₂𝑌_𝑠𝑋_𝑖𝑌_𝑖𝑖𝑌_𝑖+𝑚₃𝑌_𝑖𝑖. (55) Substituting Λ₁=𝛽₁𝑋_𝑠𝑖𝑌_𝑖𝑖+𝜇₁𝑋_𝑠𝑖−𝜎𝑋_𝑖𝑖, Λ₂=𝛽₂𝑌_𝑠𝑖𝑋_𝑖𝑖+𝜇₂𝑌_𝑠𝑖, 𝑚₁𝑋_𝑖𝑖=𝛽₁𝑋_𝑠𝑖𝑌_𝑖𝑖, and 𝑚₃𝑌_𝑖𝑖=𝛽₂𝑌_𝑠𝑖𝑋_𝑖𝑖 into 𝑑𝐿₂/𝑑𝑡, we obtain 𝑑𝐿₂=𝑑𝑡𝑙₁𝜇₁𝑋_𝑠𝑖𝑋2−^𝑠𝑋_𝑠𝑖−𝑋_𝑠𝑖𝑋_𝑠+𝛽₁𝑋_𝑠𝑖𝑌_𝑖−𝑚₁𝑋_𝑖+𝜎𝑋_𝑖−𝑋_𝑖𝑖𝑋1−^𝑠𝑖𝑋_𝑠+2𝛽₁𝑋_𝑠𝑖𝑌_𝑖𝑖−𝛽₁𝑌_𝑖𝑖𝑋_𝑠𝑖²𝑋_𝑠−𝛽₁𝑋_𝑠𝑌_𝑖𝑋_𝑖𝑖𝑋_𝑖+𝑙₂𝜇₂𝑌_𝑠𝑖𝑌2−^𝑠𝑌_𝑠𝑖−𝑌_𝑠𝑖𝑌_𝑠+𝛽₂𝑌_𝑠𝑖𝑋_𝑖−𝑚₃𝑌_𝑖+2𝛽₂𝑌_𝑠𝑖𝑋_𝑖𝑖−𝛽₂𝑋_𝑖𝑖𝑌_𝑠𝑖²𝑌_𝑠−𝛽₂𝑌_𝑠𝑋_𝑖𝑌_𝑖𝑖𝑌_𝑖. (56) Note that for the limiting system 𝑋_𝑠+𝑋_𝑖=𝑋_𝑠𝑖+𝑋_𝑖𝑖=Λ₁/𝜇₁, then (𝑋_𝑖−𝑋_𝑖𝑖)(𝑋_𝑠−𝑋_𝑠𝑖)≤0. Let 𝑙₁=𝛽₂𝑌_𝑠𝑖𝑋_𝑖𝑖 and 𝑙₂=𝛽₁𝑋_𝑠𝑖𝑌_𝑖𝑖. Following from 𝑚₁𝑋_𝑖𝑖=𝛽₁𝑋_𝑠𝑖𝑌_𝑖𝑖 and 𝑚₃𝑌_𝑖𝑖=𝛽₂𝑌_𝑠𝑖𝑋_𝑖𝑖, we have 𝑙₁𝛽₁𝑋_𝑠𝑖𝑌_𝑖−𝑚₁𝑋_𝑖+𝑙₂𝛽₂𝑌_𝑠𝑖𝑋_𝑖−𝑚₃𝑌_𝑖=0. (57) Hence, for all (𝑋_𝑠,𝑋_𝑖,𝑌_𝑠,𝑌_𝑖)∈Ω₂, 𝑑𝐿₂≤𝑑𝑡𝑙₁𝜇₁𝑋_𝑠𝑖𝑋2−^𝑠𝑋_𝑠𝑖−𝑋_𝑠𝑖𝑋_𝑠+𝛽₁𝑋_𝑠𝑖𝑌_𝑖𝑖𝑋2−^𝑠𝑖𝑋_𝑠−𝑋_𝑠𝑌_𝑖𝑋_𝑖𝑖𝑋_𝑠𝑖𝑌_𝑖𝑖𝑋_𝑖+𝑙₂𝜇₂𝑌_𝑠𝑖𝑌2−^𝑠𝑌_𝑠𝑖−𝑌_𝑠𝑖𝑌_𝑠+𝛽₂𝑌_𝑠𝑖𝑋_𝑖𝑖𝑌2−^𝑠𝑖𝑌_𝑠−𝑌_𝑠𝑋_𝑖𝑌_𝑖𝑖𝑌_𝑠𝑖𝑋_𝑖𝑖𝑌_𝑖≤𝑙₁𝜇₁𝑋_𝑠𝑖𝑋2−^𝑠𝑋_𝑠𝑖−𝑋_𝑠𝑖𝑋_𝑠+𝑙₂𝜇₂𝑌_𝑠𝑖𝑌2−^𝑠𝑌_𝑠𝑖−𝑌_𝑠𝑖𝑌_𝑠+𝑙₁𝑙₂𝑋4−^𝑠𝑖𝑋_𝑠−𝑋_𝑠𝑌_𝑖𝑋_𝑖𝑖𝑋_𝑠𝑖𝑌_𝑖𝑖𝑋_𝑖−𝑌_𝑠𝑖𝑌_𝑠−𝑌_𝑠𝑋_𝑖𝑌_𝑖𝑖𝑌_𝑠𝑖𝑋_𝑖𝑖𝑌_𝑖≤0. (58) It is easy to see that 𝑑𝐿₂𝑑𝑡=0⟺𝑋_𝑠=𝑋_𝑠𝑖,𝑋_𝑖=𝑋_𝑖𝑖,𝑌_𝑠=Y_𝑠𝑖,𝑌_𝑖=𝑌_𝑖𝑖. (59) Therefore, the only compact invariant subset of the set where 𝑑𝐿₂/𝑑𝑡=0 is the point (𝑋_𝑠𝑖,𝑋_𝑖𝑖,𝑌_𝑠𝑖,𝑌_𝑖𝑖), corresponding to the boundary equilibrium 𝐸_0𝑖. By LaSalle’s Invariance Principle, 𝐸_0𝑖 is globally asymptotically stable if 𝑅_0𝑖>1 and 𝑅_0𝑟≤1. Summering above analysis, we have the following result.

Theorem 9. When 𝜌=0, the boundary equilibrium 𝐸_0𝑖 of the system (1) is globally asymptotically stable if 𝑅_0𝑖>1≥𝑅_0𝑟.

5. Stability Analysis of the Coexistence Equilibrium

In this section, we turn to study the local stability of the coexistence equilibrium 𝐸_𝑐 in the limiting system of (1) by using Krasnoselskii sublinearity trick [24], as in [25, 26]. In detail, if 𝑥^′=𝑓(𝑥) is a system of differential equations and 𝑥^∗ is an equilibrium point, then to prove the local asymptotical stability of 𝑥^∗ is to prove that the linearized equation 𝑍^′=𝑓^′(𝑥^∗)𝑍 has no solutions of the form 𝑍(𝑡)=𝑍₀𝑒^𝜔𝑡, (60) with 𝑍₀∈𝐶^𝑛, 𝜔∈𝐶 and, Re𝜔≥0. This implies that the eigenvalues of the characteristic polynomial associated with the linearized equations have negative real part, that is, Re𝜔<0. Then, the coexistence equilibrium 𝐸_𝑐 is locally asymptotically stable.

Considering the limiting system 𝑑𝑋_𝑖𝑑𝑡=𝛽₁𝑁₁−𝑋_𝑖−𝑋_𝑟𝑌_𝑖−𝑚₁𝑋_𝑖,𝑑𝑋_𝑟𝑑𝑡=𝑘𝛽₁𝑁₁−𝑋_𝑖−𝑋_𝑟𝑌_𝑟+𝜌𝜎𝑋_𝑖−𝑚₂𝑋_𝑟,𝑑𝑌_𝑖𝑑𝑡=𝛽₂𝑁₂−𝑌_𝑖−𝑌_𝑟𝑋_𝑖−𝑚₃𝑌_𝑖,𝑑𝑌_𝑟𝑑𝑡=𝑘𝛽₂𝑁₂−𝑌_𝑖−𝑌_𝑟𝑋_𝑟−𝑚₄𝑌_𝑟, (61) where 𝑁₁=Λ₁/𝜇₁ and 𝑁₂=Λ₂/𝜇₂. In this way, let 𝑍₀=(𝑍₁,𝑍₂,𝑍₃,𝑍₄), 𝑍_𝑖∈𝐶. Substituting a solution of the form (60) into the linearized system (61) of the coexistence equilibrium 𝐸_𝑐, we obtain the following linear equations: 𝜔𝑍₁=−𝛽₁𝑌_𝑖𝑐+𝑚₁𝑍₁−𝛽₁𝑌_𝑖𝑐𝑍₂+𝛽₁𝑁₁−𝑋_𝑖𝑐−𝑋_𝑟𝑐𝑍₃,𝜔𝑍₂=𝜌𝜎−𝑘𝛽₁𝑌_𝑟𝑐𝑍₁−𝑘𝛽₁𝑌_𝑟𝑐+𝑚₂𝑍₂+𝑘𝛽₁𝑁₁−𝑋_𝑖𝑐−𝑋_𝑟𝑐𝑍₄,𝜔𝑍₃=𝛽₂𝑁₂−𝑌_𝑖𝑐−𝑌_𝑟𝑐𝑍₁−𝛽₂𝑋_𝑖𝑐+𝑚₃−𝛽₂𝑋_𝑖𝑐𝑍₄,𝜔𝑍₄=𝑘𝛽₂𝑁₂−𝑌_𝑖𝑐−𝑌_𝑟𝑐𝑍₂−𝑘𝛽₂𝑋_𝑟𝑐𝑍₃−𝑘𝛽₂𝑋_𝑟𝑐+𝑚₄𝑍₄, (62) which is equivalent to the system 11+𝑚₁𝜔+𝛽₁𝑌_𝑖𝑐𝑍₁𝛽=−¹𝑌_𝑖𝑐𝑚₁𝑍₂+𝛽₁𝑁₁−𝑋_𝑖𝑐−𝑋_𝑟𝑐𝑚₁𝑍₃,11+𝑚₂𝜔+𝑘𝛽₁𝑌_𝑟𝑐𝑍₂=𝜌𝜎−𝑘𝛽₁𝑌_𝑟𝑐𝑚₂𝑍₁+𝑘𝛽₁𝑁₁−𝑋_𝑖𝑐−𝑋_𝑟𝑐𝑚₂𝑍₄,11+𝑚₃𝜔+𝛽₂𝑋_𝑖𝑐𝑍₃=𝛽₂𝑁₂−𝑌_𝑖𝑐−𝑌_𝑟𝑐𝑚₃𝑍₁−𝛽₂𝑋_𝑖𝑐𝑚₃𝑍₄,11+𝑚₄𝜔+𝑘𝛽₂𝑋_𝑟𝑐𝑍₄=𝑘𝛽₂𝑁₂−𝑌_𝑖𝑐−𝑌_𝑟𝑐𝑚₄𝑍₂−𝑘𝛽₂𝑋_𝑟𝑐𝑚₄𝑍₃. (63) Moving all the negative terms to the left-hand side, after some manipulations we obtain the system 1+𝐹₁(𝜔)𝑍₁+1+𝐹₂(𝜔)𝑍₂=(𝐻𝑍)₁+(𝐻𝑍)₂,1+𝐹₃(𝜔)𝑍₃+1+𝐹₄(𝜔)𝑍₄=(𝐻𝑍)₃+(𝐻𝑍)₄, (64) where 𝐹₁(𝜔)=𝜔+𝛽₁𝑌_𝑖𝑐𝑚₁+𝑘𝛽₁𝑌_𝑟𝑐𝑚₂,𝐹₂(𝜔)=𝜔+𝑘𝛽₁𝑌_𝑟𝑐𝑚₂+𝛽₁𝑌_𝑖𝑐𝑚₁,𝐹₃(𝜔)=𝜔+𝛽₂𝑋_𝑖𝑐𝑚₃+𝑘𝛽₂𝑋_𝑟𝑐𝑚₄,𝐹₄(𝜔)=𝜔+𝑘𝛽₂𝑋_𝑟𝑐𝑚₄+𝛽₂𝑋_𝑖𝑐𝑚₃,𝛽𝐻=00¹𝑁₁−𝑋_𝑖𝑐−𝑋_𝑟𝑐𝑚₁0𝜌𝜎𝑚₂00𝑘𝛽₁𝑁₁−𝑋_𝑖𝑐−𝑋_𝑟𝑐𝛽₂𝑁₂−𝑌_𝑖𝑐−𝑌_𝑟𝑐𝑚₃000000𝑘𝛽₂𝑁₂−𝑌_𝑖𝑐−𝑌_𝑟𝑐𝑚₄. (65) Note that the matrix 𝐻 has nonnegative entries, and 𝐸_𝑐=(𝑋_𝑖𝑐,𝑋_𝑟𝑐,𝑌_𝑖𝑐,𝑌_𝑟𝑐) satisfies 𝐸_𝑐=𝐻𝐸_𝑐. (66)

To show that Re𝜔<0, we distinguish two cases: 𝜔=0 and 𝜔≠0. In the first case, (62) is a homogeneous linear system. Through calculations, we have the determinant of (62) is −𝛽Δ=¹𝑌_𝑖𝑐+𝑚₁−𝛽₁𝑌_𝑖𝑐𝛽₁𝑋_𝑠𝑐0𝜌𝜎−𝑘𝛽₁𝑌_𝑟𝑐−𝑘𝛽₁𝑌_𝑟𝑐+𝑚₂0𝑘𝛽₁𝑋_𝑠𝑐𝛽₂𝑌_𝑠𝑐0−𝛽₂𝑋_𝑖𝑐+𝑚₃−𝛽₂𝑋_𝑖𝑐0𝑘𝛽₂𝑌_𝑠𝑐−𝑘𝛽₂𝑋_𝑟𝑐−𝑘𝛽₂𝑋_𝑟𝑐+𝑚₄=−𝛽₁𝑌_𝑖𝑐+𝑚₁𝐴₁+𝛽₁𝑌_𝑖𝑐A₂+𝛽₁𝑋_𝑠𝑐𝐴₃, (67) where 𝐴₁=𝑚₁𝑚₃𝑘²𝛽₂𝑋_𝑖𝑐+𝑚₁𝑚²₃𝑘²−𝑚₄𝑘𝛽₁𝛽₂𝑋_𝑖𝑐𝑌_𝑟𝑐−𝑚₃𝑘²𝛽₁𝛽₂𝑋_𝑟𝑐𝑌_𝑟𝑐−𝑚₃𝑚₄𝑘𝛽₁𝑌_𝑟𝑐−𝑚₂𝑚₄𝛽₂𝑋_𝑖𝑐−𝑚₂𝑚₃𝑘𝛽₂𝑋_𝑟𝑐−𝑚₂𝑚₃𝑚₄,𝐴₂=𝜌𝜎𝑚₃𝑘𝛽₂𝑋_𝑟𝑐+𝑚₄𝛽₂𝑋_𝑖𝑐+𝑚₃𝑚₄−𝑚₄𝑘𝛽₁𝛽₂𝑋_𝑖𝑐𝑌_𝑟𝑐−𝑚₃𝑘²𝛽₁𝛽₂𝑋_𝑟𝑐𝑌_𝑟𝑐−𝑚₃𝑚₄𝑘𝛽₁𝑌_𝑟𝑐−𝑚₁𝑚₃𝑘²𝛽₂𝑋_𝑟𝑐,𝐴₃=𝑘²𝛽₁𝛽²₂𝑋_𝑠𝑐𝑌²_𝑠𝑐+𝜌𝜎𝑘𝛽²₂𝑋_𝑖𝑐𝑌_𝑠𝑐−𝑘²𝛽₁𝛽²₂𝑋_𝑖𝑐𝑌_𝑟𝑐𝑌_𝑠𝑐−𝛽₂𝑌_𝑠𝑐𝑘²𝛽₁𝛽₂𝑋_𝑟𝑐𝑌_𝑟𝑐+𝑚₄𝑘𝛽₁𝑌_𝑟𝑐+𝑚₂𝑘𝛽₂𝑋_𝑟𝑐+𝑚₂𝑚₄. (68) Then, we can obtain Δ=𝜌𝜎𝑚₃𝑘𝛽₁𝛽₂𝑋_𝑟𝑐𝑌_𝑖𝑐+𝑚₄𝛽₁𝛽₂𝑋_𝑖𝑐𝑌_𝑖𝑐+𝑚₃𝑚₄𝛽₁𝑌_𝑖𝑐+𝑚₁𝑚₃𝑘𝛽₂𝑋_𝑖𝑐+𝛽₁𝛽₂𝑋_𝑖𝑐𝑌_𝑖𝑐+𝑚₁𝛽₂𝑋_𝑖𝑐+𝑚₃𝛽₂𝑌_𝑖𝑐×𝑚₂𝑚₄−𝑚₁𝑚₃𝑘²+𝛽₁𝛽₂𝑋_𝑟𝑐𝑌_𝑖𝑐𝑚₃𝑘𝑚₂−𝑚₁𝑘+𝛽₁𝛽₂𝑋_𝑖𝑐𝑌_𝑟𝑐𝑚₁𝑘𝑚₄−𝑚₃𝑘. (69) Note that the coexistence equilibrium 𝐸_𝑐 exists under the condition 𝑅_0𝑖>𝑅_0𝑟, that is, 𝑚₂𝑚₄>𝑚₁𝑚₃𝑘². Then Δ>0 if 𝑚₂−𝑚₁𝑘≥0 and 𝑚₄−𝑚₃𝑘≥0. Hence, for the case that 𝜔=0, (62) has only the trivial solution if 𝑚₂−𝑚₁𝑘≥0 and 𝑚₄−𝑚₃𝑘≥0. This implies that 𝜔≠0.

Now assume that 𝜔≠0 and Re𝜔≥0. Let 𝐹(𝜔)=min{|1+𝐹_𝑖(𝜔)|,𝑖=1,2,3,4}. It is easy to prove that |1+𝐹_𝑖(𝜔)|>1 for all 𝑖, and therefore 𝐹(𝜔)>1. In the other hand, since the coordinates of 𝐸_𝑐 are positive, if 𝑍 is any solution of (64), then there exists a minimal positive real number 𝑠 such that 𝑍≤𝑠𝐸_𝑐, (70) where |𝑍|=(|𝑍₁|,|𝑍₂|,|𝑍₃|,|𝑍₄|), and |⋅| is the norm in 𝐶. Note that 𝑠 is also the minimal positive real number such that 𝑍₁+𝑍₂≤𝑠𝑋_𝑖𝑐+𝑋_𝑟𝑐,𝑍₃+𝑍₄≤𝑠𝑌_𝑖𝑐+𝑌_𝑟𝑐. (71) Taking norms on both side of (64) and using (66), (71), and the fact that 𝐻 is nonnegative, we obtain 𝐹(𝜔)𝑍₁+𝑍₂≤𝐻(|𝑍|)₁+𝐻(|𝑍|)₂≤𝑠𝑋_𝑖𝑐+𝑋_𝑟𝑐, (72) that is, 𝑍₁+𝑍₂≤𝑠𝐹(𝜔)𝑋_𝑖𝑐+𝑋_𝑟𝑐<𝑠𝑋_𝑖𝑐+𝑋_𝑟𝑐, (73) which contradicts the minimality of 𝑠. Hence, Re𝜔<0.

Summering above analysis, we have the following result.

Theorem 10. When 𝜌≠0, the coexistence equilibrium 𝐸_𝑐 of the system (1) is locally asymptotically stable if 𝑘≤min{𝑚₂/𝑚₁,𝑚₄/𝑚₃} and 𝑅_0𝑖>1.

The existence and stability of equilibria can be summered in Table 1.

Table 1 

Existence and stability of equilibria. The sign “∃  !” means the existence and uniqueness. LAS means locally asymptotically stable, and GAS means globally asymptotically stable.

6. Discussion

In this paper we established a new schistosomiasis model. In contrast to previous schistosomiasis models with drug resistance, the model established in this study consider many aspects. First, snail is considered as a variable in the model since the resistance of schistosoma to PZQ can be expressed in snails [14–16]. Second, previous models considered that resistance was caused by drug treatment, while a large number of the literature show that some cases are due to inheritance [1, 3, 4]. Therefore, in this study, we consider the reasons for resistance are drug treatment and genetic. Last, in previous models a resistance level (𝜃) was used to discount the treatment rate. But the resistance level could not be measured, and its value was assumed. This model consider what percentage of infected human after treatment will recover and what percentage will emerge drug resistance. In medicine, from the occurrence of cases the value of this ratio can be identified. For example, this ratio was given in [27, 28] (1−𝜌=0.28∼0.609). Therefore, it is easy to operate. In addition, we separate the discussion of the case that treatment will cause resistance (𝜌≠0) and treatment will not cause resistance (𝜌=0).

The reproductive number 𝑅_0𝑖 and 𝑅_0𝑟 of the sensitive and resistant strains are given, respectively. It is easy to see that 𝑅_0𝑖 is a decreasing function of 𝜎, and 𝑅_0𝑟 is an increasing function of 𝑘. If the basic reproductive number 𝑅₀=max{𝑅_0𝑖,𝑅_0𝑟} of the model (1) is less than 1, one can prove the stability of the disease free equilibrium. This means that the spread of schistosomiasis can be effectively controlled. When the basic reproductive number 𝑅₀ is greater than 1, we first consider the case that 𝜌=0. If 𝑅_0𝑖=𝑅_0𝑟>1, the two strains can coexist. There is a line of coexistence equilibria in this case. The infected human will evolve to one of them with higher reproductive number (see Table 1).

When drug treatment can not cause resistance (𝜌=0), that is, the new resistant strain is due to the inheritance of resistance, the sensitive strain will dominate if the treatment rate is smaller and the impact of the inheritance of resistance is smaller such that 𝑅_0𝑖>1≥𝑅_0𝑟. This result accords to the results of previous models. On the other hand, if the treatment rate and the impact of the inheritance of resistance are both larger such that 𝑅_0𝑟>1≥𝑅_0𝑖, the resistant strain will dominate. This shows that although the assumption is that drug treatment does not result in the emergence of drug resistance, once the treatment rate is greater than a value, and the impact of genetic resistance is larger, there will still be the emergence of resistant strain, and the resistant strain is dominant. This further shows that genetic resistance has a great impact on the system.

When drug treatment can cause resistance (𝜌≠0), we can show that the sensitive strain either does not appear or coexist with the resistant strain under certain condition. If 𝑅_0𝑟>1, there is only resistant strain. It can be seen that, regardless of whether drug treatment leads to the emergence of resistance, once the impact of genetic resistance is larger, resistant strain will be dominant, which is detrimental to the control of schistosomiasis.

Finally, from the formula of the coexistence equilibrium 𝐸_𝑐, it is easy to see that the value of the resistant strain is increased with the value of 𝜌. This means once the proportion of human with drug-resistant strain produced by drug treatment is larger, the number of human and snails with resistant strain is larger. This is not a good result for drug treatment with praziquantel.

Hence, for poor treatment, there are two possible reasons: drug therapy and genetic. An important priority in developing new control strategies is to search new drug targets, in combination with selection of genetic methods such as that viable vaccine candidates. And there is already a need for alternative drugs to treat PZQ-resistant schistosomiasis, such as already exists in northern Senegal [10].

Acknowledgments

This paper is supported by National Natural Science Foundation of China (11126177 and 11071011), Natural Science Foundation of Anhui Province (1208085QA15), and by the Foundation for Young Talents in College of Anhui Province (2012SQRL021) and also supported by Doctoral Fund of Ministry of Education of China (20113401110001 and 20103401120002), the Key Natural Science Foundation of the Anhui Higher Education Institutions of China (KJ2009A49), and Funding Project for Academic Human Resources Development in Institutions of Higher Learning Under the Jurisdiction of Beijing Municipality (PHR201107123). The authors would like to thank anonymous reviewers for very helpful suggestions which improved greatly this paper.