Optimal control strategies for a two-group epidemic model with vaccination-resource constraints

Applied Mathematics and Computation 371 (2019) 124956

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Optimal control strategies for a two-group epidemic model with vaccination-resource constraints Heting Zhang a, Zhanwen Yang a, Kasia A Pawelek b, Shengqiang Liu a,c,∗ a

School of Mathematics, Harbin Institute of Technology, Harbin, 150001, China Department of Mathematics, University of South Carolina Beaufort, Bluffton, 29909, USA c School of Mathematical Sciences, Tiangong University, Tianjin, 300387, China b

a r t i c l e

i n f o

Article history: Received 29 April 2019 Revised 25 September 2019 Accepted 1 December 2019

MSC: 34C25 92C45 92C60 Keywords: Epidemic model Vaccination Basic reproduction number Global stability Optimal control Pontryagin’s maximum principle

a b s t r a c t In a resource-limited environment, we propose and study an epidemic model consisting of two groups of hosts to reflect the different risk levels of becoming infected and transmitting infection. To further investigate vaccine implementation policies, we study the global dynamics of the model with group-dependent constant coverage rates by utilizing Lyapunov functions. We find that when the vaccination-dependent basic reproduction number is less than one (0 < 1), the disease-free equilibrium is globally asymptotically stable. Otherwise, if 0 > 1, the endemic equilibrium is globally attractive. For the case when the disease is endemic, we apply the optimal control theory to the model under the assumption of a resource-constrained environment and we derive optimal coverage strategies. We also perform numerical simulations for the optimal control problem to verify and further support the theoretical results. In limited control resource cases with both constant and variable coverage rate strategies, we show that when the total control resources are rather scarce, the resources should be biased towards at the group with the high level of infection risk; and that the control resources should preferably be directed at the lower risk group if the total resources are relatively abundant. © 2019 Elsevier Inc. All rights reserved.

1. Introduction A relatively infectious disease, known as an epidemic disease, can spread very rapidly once it erupts in a society. An outbreak, may lead to a serious prevalence of the disease in a population and, in some cases, result in deaths. Among various epidemic diseases, influenza is typical in significant ways, and its prevention has drawn great attention from both academic studies and the medical industry. Since Edward Jenner discovered the smallpox vaccine in the eighteenth century, vaccine research has made great progress and has led to the development of medical science in the past two centuries [2,7–9]. Nevertheless, in practice the situation can be highly complex. For instance, school children and the elderly usually become infected more easily than other age group [10,12,13,43]. In [15,16], certain groups of individuals were shown to spread the measles or HIV more rapidly than other groups. Hence, allocation of an efficient vaccination program among these groups characterized as risky and critical is of great importance. During recent years a number of studies on the spread of infectious diseases among individuals of heterogeneous groups have been carried out [1,21,48,49]. Moreover, mathematical models ∗

Corresponding author. E-mail addresses: [email protected] (H. Zhang), [email protected] (Z. Yang), [email protected] (K.A. Pawelek), [email protected] (S. Liu). https://doi.org/10.1016/j.amc.2019.124956 0 096-30 03/© 2019 Elsevier Inc. All rights reserved.

2

H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956

have been developed to predict the course of an infectious disease [1,5,6,22–24] and to investigate the effects of vaccination, quarantine-isolation, or other control and prevention strategies [3,11,14,20]. It is of great importance to categorize the host into different groups to measure the outcome of the vaccination, as the infection may vary widely among the individuals. Motivated by these points, in [1], Ashrafur and Zou proposed a mathematical model which grouped individuals with a high risk of getting infected in a specific susceptible group. They proposed a mathematical model with group-dependent incidence and vaccination rates to evaluate this measure. Their results offered guidelines on vaccine implementations focusing on effectiveness strategies and cost efficacy. Although much progress has been made in [1] and the above references in investigating effective vaccination strategies, some assumptions of earlier models should be relaxed so that more realistic situations can be studied. In particular, most known results are based on the assumption that the vaccine is fully protective and vaccine rates are constant. However, such assumptions are not realistic for many infectious diseases. For instance, the minimum recommended interval between the hepatitis B vaccine doses is four weeks (World Health Organization [42]). Halota et al. [4] showed that after the second dose, 12.7%, and after the third dose, 1.3% of patients did not achieve an anti-HBs protective titer ( > 410 mIU/ml), which is considered a reliable marker of immediate and long-term protection against an infection. Furthermore, there is evidence that the vaccine supply is usually limited or not readily available for emerging diseases; as shown for the H1N1 influenza A vaccine during 2009–2010 pandemic season. In addition, public health faced challenges from the differences in group-specific susceptibility and group-fatality during this influenza A pandemic. Therefore, it is important for public health purposes to determine an optimal vaccination policy incorporating the cost as well as the benefits of vaccination [25,26]. The optimal control vaccine method has been characterized as an efficient treatment approach among the various prevention strategies and has been well investigated [27,28]. To our knowledge, the optimal vaccination policy problem for a multi-group epidemic model in a resource-limited situation with non-fully immune vaccinated compartment has not been investigated thoroughly, and has motivated this study. In this paper, we will incorporate vaccination control into a two-group compartmental deterministic epidemic model with non-fully immune vaccinated compartment, and our main purpose is to investigate the optimal vaccination policy problem in the resource-limited environment for both the situations of constant and variant coverage rates. The rest of this paper is organized as follows. In Section 2, a two-group SIVI epidemic model with constant groupdependent coverage is given and some preliminary results are obtained. The main results of this model are presented in Section 3. Section 4 extends the model of Section 2 into the variable coverage and investigates the optimal coverage strategies. Numerical simulations are shown in Section 5 to verify the effectiveness of the proposed approach. Finally, the conclusions are given in Section 6. 2. Model formulation The total population is divided into the following two groups: Group R: the risky group with higher infection rates than others in the group; Group C: the critical group which contains the individuals who change their behavior and act with caution during an outbreak, and individuals who avoid contact while infected; hence reducing the rate at which they can spread the infection. Assumption 2.1. For susceptible, vaccinated individuals, infected, and recovered individuals, we make the following assumptions: (1) Susceptible individuals can become infected due to contact with infected individuals. These infected individuals remain the infected state until they recover; (2) Susceptible individuals can be vaccinated and move to the vaccinated state however they may still become infected; (3) We do not consider the vertical infection. The symbols and notation used in the model are explained in Table 1. Schematic representation of the model (2.1) can be shown in Fig. 6.1. The infection mechanism is considered to be Ir Ic followed by saturating incidence [17–19] defined by fr = , fc = where α ≥ 0 determines the saturation 1 + αr I r 1 + αc I c level when the infectious population is large. When α = 0,this reduces to the mass action incidence rate. The infection rate increases with the number of infected individuals when this number is small. As the infected number increases, the infection rate becomes plateaus. This phenomenon reflects the saturation of infected numbers also known as ‘crowding effect’ [1]. For convenience, we set θ i = θi ηi , i = r, c. With this assumption, the dynamics of the population are governed by the following equations:





dSr (t ) Ir Ic = r − βrr + βcr S r − ( μ + θ r )S r , r dt 1 + αr I 1 + αc Ic   dSc (t ) Ir Ic = c − βrc + βcc S c − ( μ + θ c )S c , r c dt 1 + αr I 1 + αc I dIr (t ) = dt

 βrr

Ir Ic + βcr 1 + αr Ir 1 + αc Ic



 S + κrr r



Ir Ic + κcr V r − ( μ + vr + γ r )I r , 1 + αr Ir 1 + αc Ic

H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956

3

Table 1 Definitions of the parameters and variables. Symbol r

S (t) Sc (t) Ir (t) Ic (t) Vr (t) Vc (t) Rr (t) Rc (t)

βij κij vi

μ θi γi ηi

Meaning number of risky susceptible individuals in t time number of critical susceptible individuals in t time number of risky infected individuals in t time number of critical infected individuals in t time number of risky vaccinated individuals in t time number of critical vaccinated individuals in t time number of risky recovered individuals in t time number of critical recovered individuals in t time contact rate between susceptible individuals and infected individuals (i, j = r, c ) contact rate between vaccinated individuals and infected individuals (i, j = r, c ) disease-induced death rate (i, j = r, c ) natural death rate rate at which susceptible individuals become successfully vaccinated (i = r, c ) recovery rate (i = r, c ) coverage rate (i = r, c )

Fig. 6.1. Schematic representation of the model (2.1).

  Ir Ic c + κ V c − ( μ + vc + γ c )I c , κrc c 1 + αr Ir 1 + αc Ic   dV r (t ) Ir Ic r = θ r Sr − κrr + κ V r − μV r , c dt 1 + αr Ir 1 + αc Ic   dV c (t ) Ir Ic c = θ c Sc − κrc + κ V c − μV c , c dt 1 + αr Ir 1 + αc Ic dRr (t ) = γ r I r − μR r , dIc (t ) = dt

dt dRc (t ) = dt

 βrc

Ir Ic + βcc r 1 + αr I 1 + αc Ic



Sc +

γ c I c − μR c .

(2.1)

According to the various contact ways of individuals in this paper we take the following:

βrr  βcc , βrr  κrr , βcc  κcc

(2.2)

Thereinto, • βrr  βcc indicates that the contact rate between risky susceptible individuals and risky infected individuals is far more than that between critical susceptible individuals and critical infected individuals; • βrr  κrr indicates that the contact rate between risky susceptible individuals and risky infected individuals is far more than that between risky vaccinated individuals and risky infected individuals; • βcc  κcc indicates that the contact rate between critical susceptible individuals and critical infected individuals is far more than that between critical vaccinated individuals and critical infected individuals;

4

H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956

All of the equations in model (2.1) except consider the following reduced model of (2.1):

dRr dRc and being independent of Rr (t) and Rc (t), then it suffices for us to dt dt

 Ir Ic r βrr + β S r − ( μ + θ r )S r , c 1 + αr Ir 1 + αc Ic   dSc (t ) Ir Ic c = c − βrc + β S c − ( μ + θ c )S c , c dt 1 + αr Ir 1 + αc Ic     dIr (t ) Ir Ic Ir Ic r r r r κ = βrr + β S + + κ V r − ( μ + vr )I r , c r c dt 1 + αr Ir 1 + αc Ic 1 + αr Ir 1 + αc Ic     dIc (t ) Ir Ic Ir Ic c c c c c = βr + βc S + κr + κc V c − ( μ + vc )I c , dt 1 + αr Ir 1 + αc Ic 1 + αr Ir 1 + αc Ic   dV r (t ) Ir Ic r r r r = θ S − κr + κc V r − μV r , dt 1 + αr Ir 1 + αc Ic   dV c (t ) Ir Ic c c c c = θ S − κr + κc V c − μV c . dt 1 + αr Ir 1 + αc Ic dSr (t ) = r − dt



(2.3)

2.1. Positivity and boundary of solutions In this subsection, we will investigate the positivity and boundedness of the solutions of model (2.3) based on epidemic reasons. Theorem 2.1. For t > 0, the solutions of the model (2.3) which initial values are Sr (0) ≥ 0, Sc (0) ≥ 0, Ir (0) > 0, Ic (0) > 0, Vr (0) ≥ 0, Vc (0) ≥ 0 are positive in R6 . Proof. The first equation can be written as

dSr (t ) = r − φ (t )Sr (t ) dt and

φ (t ) = βrr

Ir Ic + βcr + μ + θ r. r 1 + αr I 1 + αc Ic

As long as S0r ≥ 0, on the basis of the method of variation of constant, we can get 



t t Sr (t ) = S0r e− 0 φ (s )ds + r e− 0 φ (s )ds



t

0



τ e 0 φ (s )ds dτ ,

thus Sr (t) > 0, ∀t > 0. Similarly, one can show that Sc (t) > 0, ∀t > 0. Next we will prove Vr (t) > 0, ∀t > 0. In fact, according to the fifth equation in model (2.3), we have 



t t V r (t ) = V0r e− 0 ψ (s )ds + e− 0 ψ (s )ds

r



t 0

τ

θ r Sr (τ )e 0 ψ (s)ds dτ ,

(2.4)

c

where ψ (t ) = κrr 1+Iαr Ir + κcr 1+Iαc Ic + μ. Then as long as V0r is non-negative, we have Vr (t), ∀t > 0. Similarly arguments can be done for Vc (t). To show that Ir (t) and Ic (t) are always positive for any positive t, we choose to prove by contradiction. If not, there must exist some t0 > 0, s.t. min{Ir (t0 ), Ic (t0 )} = 0 while min {Ir (t), Ic (t)} > 0, ∀t ∈ [0, t0 ]. However, from the third and r forth equations of model (2.3) and the positiveness of Sr (t), Vr (t), we have dIdt(t ) ≥ −(μ + vr )Ir , ∀t ∈ [0, t0 ], thus we have Ir (t ) ≥ Ir (0 ) · e−(μ+v )t > 0, ∀t ∈ [0, t0 ]; Similarly, we get Ic (t ) ≥ Ic (0 ) · e−(μ+v )t > 0, ∀t ∈ [0, t0 ]. A contradiction to min{Ir (t0 ), Ic (t0 )} = 0, proving the positiveness of Ir (t), Ic (t). This completes the proof of Theorem 2.1.  r

c

Now, we are concerned with the boundary of the model (2.3). Theorem 2.2. All positive solutions of the system (2.3) have the ultimate upper boundary in R6 . Proof. By summing all of the equations in (2.3), it can be seen that the total number of individuals satisfies:

dSr (t ) dSc (t ) dIr (t ) dIc (t ) dV r (t ) dV c (t ) + + + + + dt dt dt dt dt dt = r + c − μSr (t ) − μSc (t ) − (μ + vr )Ir (t ) − (μ + vc )Ic (t ) − μV r (t ) − μV c (t ) therefore we obtain:

dSr (t ) dSc (t ) dIr (t ) dIc (t ) dV r (t ) + + + + ≤ r + c − μ(Sr (t ) + Sc (t ) + Ir (t ) + Ic (t ) + V r (t ) + V c (t ) ) dt dt dt dt dt

H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956

5

Method of constant variation provides:

lim sup(Sr (t ) + Sc (t ) + Ir (t ) + Ic (t )V r (t ) + V c (t )) ≤ t→∞

r + c μ 

So the biologically feasible region of the model (2.3) is



r + c = (Sr (t ), Sc (t ), Ir (t ), Ic (t ), V r (t ), V c (t )) ∈ R6+ : Sr (t ) + Sc (t ) + Ir (t ) + Ic (t ) + V r (t ) + V c (t ) ≤ μ

(2.5)

3. Main results The model (2.3) has a disease-free equilibrium (DFE), denoted by:



r c θ r r θ c c , , 0, 0, 2 , 2 r c r μ+θ μ+θ μ + μθ μ + μθ c

E0 =



The stability of E0 is closely related to the concept of the model basic reproduction number. In this paper, we use 0 to denote the basic reproduction number and we will prove that it plays an important role in determining the disease persistence. The number 0 is defined as ‘the expected number of secondary cases produced, in a completely susceptible population, by a typically infected individual’ (see, e.g. [39]). Next we use next-generation matrix [40] to compute 0 . To this end, we define the non-negative matrix F and the non-singular M-matrix M, represented as new-infection and transition matrices severally, for the system (2.3), are provided by

⎛

r θ r r ⎞ r + κ c μ + θr μ2 + μθ r ⎟ ⎠ c  θ c c c βcc + κ c μ + θc μ2 + μθ c

r θ r r + κrr 2 r μ + θ μ + μθ r ⎜ F =⎝ c  θ c c βrc + κrc 2 c μ+θ μ + μθ c   μ + vr 0 M= 0 μ + vc

βrr

βcr

(3.1)

(3.2)

The definition of the basic reproduction number is



0 = ρ ( F M

−1

)=ρ

F J

H K



=

F +K+



(F − K )2 + 4HJ

(3.3)

2

where

r 1 θ r r + κr , (μ + θ r )(μ + vr ) μ + vr r μ2 + μθ r

H = βcr

r 1 θ r r + κr , (μ + θ r )(μ + vc ) μ + vc c μ2 + μθ r

c 1 θ c c + κrc 2 , r r (μ + θ )(μ + v ) μ + v μ + μθ c

K = βcc

c 1 θ c c + κcc 2 . c c (μ + θ )(μ + v ) μ + v μ + μθ c

F = βrr J = βrc

c

c

By Theorem 2 in [40], we can gain the following result on the stability of E0 . Theorem 3.1. For system (2.3), if 0 < 1, the DFE E0 is locally asymptotically stable in ; it changes into unstable if 0 > 1. 3.1. Global stability of E0 In this section, we investigate the global stability of the DFE E0 for the model (2.3). Theorem 3.2. For system (2.3), if 0 < 1, E0 is globally asymptotically stable in . Proof. In the model (2.3), we can obtain

lim sup Sr (t ) ≤ t→∞

dSr (t ) ≤ r − (μ + θ r )Sr (t ) from the equation about Sr (t) . This indicates dt

r μ + θr

(3.4)

In the same way

lim sup Sc (t ) ≤ t→∞

c , μ + θc

lim sup V r (t ) ≤ t→∞

θ r r , μ + μθ r 2

lim sup V c (t ) ≤ + t→∞

θ c c μ + μθ c 2

(3.5)

Thus, for any ε > 0, there exists T1 > 0, when t ≥ T1 , such that

Sr (t ) ≤

r + ε μ + θr

(3.6)

6

H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956

and Sc (t) also satisfies the conditions, we obtain

c + ε μ + θc

Sc (t ) ≤

(3.7)

Analogously,

V r (t ) ≤

θ r r + ε , μ2 + μθ r

V c (t ) ≤

θ c c + ε μ2 + μθ c

(3.8)

Applying the estimation in Eqs. (3.6)–(3.8), we can receive

⎞ θ r r + ε θ c c + ε r c r r ⎜ β (t ) + β (t ) μ + θ r + κ (t ) μ2 + μθ r + κc I (t ) μ2 + μθ c − (μ + v )I (t )⎟ ⎟ ⎝ c ⎠≤ ⎜  c + ε ⎝ ⎠ r r c c θ  + ε θ  + ε dI (t ) c r c c c c βrc Ir (t ) + βcc Ic (t ) + κr I (t ) 2 + κc I (t ) 2 − (μ + v )I (t ) c r c μ+θ μ + μθ μ + μθ dt    r P W I (t ) ≤ G Q Ic (t ) ⎛

⎞

dIr (t ) ⎜ dt ⎟

⎛

r r rI

r c cI

 r + ε

r r rI

where

P = βrr

r + ε θ r r + ε + κrr 2 − ( μ + vr ), μ + θr μ + μθ r

G = βrc

c + ε θ r r + ε + κrc 2 , c μ+θ μ + μθ r

Q = βcc

W = βcr

r + ε θ c c + ε + κcr 2 μ + θr μ + μθ c

c + ε θ c c + ε + κcc 2 − ( μ + vc ) c μ+θ μ + μθ c

Thus, the sub-system (2.3) has a linear and cooperative comparison system which has upper bound with following coefficient matrix:



P G

A (ε ) =

W Q



(3.9)

As t → ∞, we can choose such sufficiently small ε that s(A(ε )) < 0, where s(A(ε )) denotes the maximum real part of all the eigenvalues of the matrix A(ε ) (the spectral abscissa of A(ε ) [37]). Then, all solutions of this comparing linear system tend to (0, 0)T . Its parameters Ir and Ic also approach 0. When t → ∞, we conclude that Ir (t) → 0 and Ic (t) → 0, this indicates that the sub-system of the Eqs. (2.3) have the following limit system:

⎧ dSr (t ) ⎪ = r − (μ + θ r )Sr (t ), ⎪ ⎪ dt ⎪ ⎪ ⎪ ⎪ dSc (t ) ⎪ ⎪ = c − (μ + θ c )Sc (t ), ⎨ dt

(3.10)

⎪ dV r (t ) ⎪ ⎪ = θ r Sr − μV r , ⎪ ⎪ dt ⎪ ⎪ ⎪ ⎪ ⎩ dV c (t ) = θ c Sc − μV c . dt

r +ε c +ε , V r (t ) → θ r r +ε , V c (t ) → θ c c +ε . Though any soluc When t → ∞, we obtain that Sr (t ) →  μ+θ r , S (t ) → μ+θ c  μ2 +μθ r μ2 +μθ c  r  , c , 0 , 0 , θ r r , θ c c tion of Eq. (2.3) tends to , b y the theory of asymptotically autonomous sysμ+ θ r μ+ θ c μ2 +μθ r μ2 +μθ c



tems, the



r

c

Sr (t ), Sc (t ), Ir (t ), Ic (t ), V r (t ), V c (t ) θ r r

θ c c





part of any non-negative solution of Eq. (2.3) also tends to

μ+θ r , μ+θ c , 0, 0, μ2 +μθ r , μ2 +μθ c . Hence every non-negative solution of Eq. (2.3) converges to the (DFE) E0 . The global

attractiveness of E0 and the local stability result in the global asymptotical stability of E0 , accomplishing the proof of the theorem.  3.2. Persistence of the disease When 0 > 1, DFE becomes unstable, therefore it is natural to think that under this conditions the infectious populations Ir (t) and Ic (t) will keep continuous. In this subsection, we confirm this expectation. Indeed, we will prove the following theorem. Theorem 3.3. Hypothesize that the system (2.3)is uniformly persistent if 0 > 1; that is, there exists a positive constant η > 0 such that all solutions of (2.3) satisfy:

lim inf Ir (t ) > η, t→∞

lim inf Ic (t ) > η t→∞

H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956

7

Furthermore, there exists an endemic equilibrium on this occasion. Proof. We shall apply Theorem 3.4 in [36] to prove the uniform persistence. To this end, we set

X = {(Sr (t ), Sc (t ), Ir (t ), Ic (t ), V r (t ), V c (t )) ∈ R6+ : Ir (0 ) ≥ 0, Ic (0 ) ≥ 0} X0 = {(Sr (t ), Sc (t ), Ir (t ), Ic (t ), V r (t ), V c (t )) ∈ X : Ir (0 ) > 0, Ic (0 ) > 0}

∂ X0 = \ X0 = {(Sr (t ), Sc (t ), 0, 0, V r (t ), V c (t )) ∈ X : Ir (0 ) = 0 or Ic (0 ) = 0} which is relatively closed in X. Let (t) : X → X be the solution flow associated with system (2.3), that is, (t )(A0 ) = A(t ). By the positivity, it is easy to know that (t) is positively invariant about X when Sr (0) ≥ 0, Sc (0) ≥ 0, Ir (0) > 0, Ic (0) > 0, Vr (0) ≥ 0, Vc (0) ≥ 0, then Sr (t) > 0, Sc (t) > 0, Ir (t) > 0, Ic (t) > 0, Vr (t) > 0, Vc (t) > 0 for ∀t > 0. Now we prove X0 is positively invariant for (t). By the third and forth equations of (1) we have,

dIr (t ) ≥ −(μ + vr )Ir (t ), dt

∀t > 0

(3.11)

dIc (t ) ≥ −(μ + vr )Ic (t ), dt

∀t > 0

(3.12)

If Ir (0) ≥ 0, Ic (0) ≥ 0, it follows from (3.11) and (3.12) that r Ir (t ) ≥ Ir (0 ) · e−(μ+v )t ,

c Ic (t ) ≥ Ic (0 ) · e−(μ+v )t ,

∀t > 0

(3.13)

Thus X0 is positively invariant for (t). Moreover, by Theorem 2.2, system (2.3) has the ultimate upper boundary, thus we have X is point dissipative for (t). We set

M∂ = {A0 ∈ ∂ X0 , ∀t ≥ 0}

(3.14)

We claim that

M∂ = { ( S r , S c , 0 , 0 , V r , V c )}

(3.15)

Assuming A(t) ∈ M∂ , ∀t ≥ 0, it suffices to show that = that either (a) Ir (t0 ) > 0, Ic (t0 ) = 0; or (b) Ir (t0 ) = 0, Ic (t0 ) > 0. For case (a), from the forth equation of system (2.3) we have Ir (t )

Ic (t )

= 0, ∀t ≥ 0. If it is not true, then there exists t0 > 0 such

dIc (t ) Ir (t0 ) Ir (t0 ) |t=t0 = βrc Sc (t0 ) + κrc V c (t0 ) > 0 r dt 1 + αr I (t0 ) 1 + αr Ir (t0 )

(3.16)

Hence there is an ε 0 > 0 such that Ic (t) > 0, ∀t ∈ (t0 , t0 + ε0 ). On the other hand, from Ir (t0 ) > 0 there exists an ε 1 (0 < ε 1 < ε 0 ) such that Ir (t) > 0, ∀t ∈ (t0 , t0 + ε1 ). Thus, we have Ir (t) > 0, Ic (t) > 0, ∀t ∈ (t0 , t0 + ε1 ), which contradicts the assumption that {Sr (t), Sc (t), Ir (t), Ic (t), Vr (t), Vc (t) ∈ M∂ , ∀t ≥ 0}. Similarly, we can obtain a contradiction for case (b). This proves the claim (3.15).  Let A = x∈A ω (x ), where Ab is the global attractor of system (2.3) restricted to ∂ X0 . We show that A = E0 . In fact, from b A ⊆ M∂ and the first and second equations of system (2.3), we have

lim Sr (t ) =

t→∞

r + ε c + ε θ r r + ε θ c c + ε , lim Sc (t ) = , lim V r (t ) = 2 , lim V c (t ) = 2 . r c r t→∞ t→∞ μ+θ μ+θ μ + μθ t→∞ μ + μθ c

Thus, E0 is the isolated invariant set in X. Next we show that W s (E0 ) ∩ X0 = ∅. If this is not true, then there exists a solution

(Sr (t ), Sc (t ), Ir (t ), Ic (t ), V r (t ), V c (t )) ∈ X0 such that



lim (S (t ), S (t ), I (t ), I (t ), V (t ), V (t )) = r

t→∞

c

r

c

r

c

 r c θ r r θ c c , , 0, 0, 2 , . μ + θr μ + θc μ + μθ r μ2 + μθ c

Hence for any sufficiently small constant ξ > 0, there exists a positive constant T1 = T1 (ξ ) such that we have the following inequalities for all t ≥ T1 :

r − ξ r + ξ c − ξ c + ξ ≤ Sr (t ) ≤ , ≤ Sc (t ) ≤ , r r c μ+θ μ+θ μ+θ μ + θc 0 ≤ Ir (t ) ≤ ξ , 0 ≤ Ic (t ) ≤ ξ ,

8

H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956

θ r r − ξ θ r r + ξ ≤ V r (t ) ≤ 2 , 2 r μ + μθ μ + μθ r

θ c c − ξ θ c c + ξ ≤ V c (t ) ≤ 2 . 2 c μ + μθ μ + μθ c

(3.17)

Consider the following auxiliary system with the constant ξ given above:

    Ir (t ) Ic (t ) Ir (t ) Ic (t ) r r r r + β S ( t ) + + κ V r (t ) − (μ + vr )Ir (t ) βrr κ c r c 1 + αr Ir (t ) 1 + αc Ic (t ) 1 + αr Ir (t ) 1 + αc Ic (t )   r − ξ θ r r − ξ r r ≥ βrr + κ − ( μ + v ) Ir (t ) r (1 + αr ξ )(μ + θ r ) (1 + αr ξ )(μ2 + μθ r )   c − ξ θ r r − ξ r + βcr + κ Ic (t ); c (1 + αc ξ )(μ + θ c ) (1 + αc ξ )(μ2 + μθ r )

dIr (t ) = dt

    Ir (t ) Ic (t ) Ir (t ) Ic (t ) c c c c + β S ( t ) + + κ V c (t ) − (μ + vc )Ic (t ) βrc κ c r c 1 + αr Ir (t ) 1 + αc Ic (t ) 1 + αr Ir (t ) 1 + αc Ic (t )   c − ξ θ c c − ξ c ≥ βrc + κ Ir (t ) r (1 + αr ξ )(μ + θ c ) (1 + αr ξ )(μ2 + μθ c )   c − ξ θ c c − ξ c c + βcc + κ − ( μ + v ) Ic (t ). c (1 + αc ξ )(μ + θ c ) (1 + αc ξ )(μ2 + μθ c )

dIc (t ) = dt

Considering the following comparing system, where I˜r (0 ) = Ir (0 ), I˜c (0 ) = Ic (0 ):

⎧   d Ir (t ) r − ξ θ r r − ξ ⎪ r r r ⎪ = + κ − ( μ + v ) Ir (t ) β ⎪ r r ⎪ dt (1 + αr ξ )(μ + θ r ) (1 + αr ξ )(μ2 + μθ r ) ⎪ ⎪   ⎪ ⎪ ⎪ c − ξ θ r r − ξ ⎪ ⎪ + βcr + κcr Ic (t ) ⎨ (1 + αc ξ )(μ + θ c ) (1 + αc ξ )(μ2 + μθ r )   c c c ⎪ d Ic (t )  −ξ θ  −ξ ⎪ c c ⎪ = βr + κr Ir (t ) ⎪ ⎪ dt (1 + αr ξ )(μ + θ c ) (1 + αr ξ )(μ2 + μθ c ) ⎪ ⎪   ⎪ ⎪ c − ξ θ c c − ξ ⎪ c c ⎪ + βcc + κ − ( μ + v ) Ic (t ) ⎩ c (1 + αc ξ )(μ + θ c ) (1 + αc ξ )(μ2 + μθ c )

(3.18)

, represented as new-infection and transition matrices Then, the non-negative matrix F and the non-singular M-matrix M severally, for the system (3.18), are provided by

⎛

r − ξ θ r r − ξ + κrr r (1 + αr ξ )(μ2 + μθ r ) ⎜ (1 + αr ξ )(μ + θ ) F = ⎝ c − ξ θ c c − ξ βrc + κrc (1 + αr ξ )(μ + θ c ) (1 + αr ξ )(μ2 + μθ c )

βrr

⎞ c − ξ θ r r − ξ + κcr c 2 r (1 + αc ξ )(μ + θ ) (1 + αc ξ )(μ + μθ ) ⎟ ⎠ c  −ξ θ c c − ξ c βcc + κ c (1 + αc ξ )(μ + θ c ) (1 + αc ξ )(μ2 + μθ c ) βcr

(3.19) = M

 μ + vr 0

0 μ + vc

−1 = Denote  J (ξ ) = FM −ξ κcr (1+α ξ )(θμ2+μθ r )(μ+vc ) , r

c

r

−ξ κcc (1+α ξ )(θμ2+μθ c )(μ+vc ) . c

c



A C

 (3.20)



r r r −ξ B c −ξ r , where A = βrr (1+α ξ )(μ+−θξr )(μ+vr ) + κrr (1+α ξ )(θμ2+μθ r )(μ+vr ) , B = βc (1+αc ξ )(μ+θ c )(μ+vc ) + r r D

−ξ C = βrc (1+α ξ )(μ+−θξc )(μ+vr ) + κrc (1+α ξ )(θμ2+μθ c )(μ+vr ) , r r c

c

c

D = βcc (1+α ξ )(μ+−θξc )(μ+vc ) + c c

c

By Lemma 2.1 in [32], as long as 0 > 1, when ξ and s(J˜(0 )) > 0, s(J˜(ξ )) is continuous. We can select such sufficiently small ξ > 0 that obtain s(J˜(ξ )) > 0 and this can explain that the positive solutions of the lower comparing system increase exponentially. By contrasting to the standard comparing system, as t → ∞, the solutions Ir (t) and Ic (t) of Eq. (2.3) are infinite which are opposite to the truth that the solutions of the system (2.3) are bounded. Therefore, W s (E0 ) ∩ X0 = ∅. Clearly, every orbit in M∂ converges to E0 . Thus by Theorem 3 in [35], we have

lim inf(Ir (t ), Ic (t )) > (η, η ), t→∞

f or some constant

η > 0.

(3.21)

By Theorem 4.3 and Remark 4.3 in [33], we conclude that system (2.3) is uniformly persistent with respect to (X0 , ∂ X0 ). And from Theorem 2.4 in [34], system (2.3) has at least one equilibrium (S∗r , S∗c , I∗r , I∗c , V∗r , V∗c ) ∈ X0 , with I∗r  0 and I∗c  0. As a result, there exists the endemic equilibrium in the system (2.3). The proof of the theorem is complete. 

H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956

9

3.3. Global stability of E∗ In this subsection, we research the global stability of the endemic equilibrium E∗ under the condition which is 0 > 1. For convenience, we simplify • Ir = Ir (t ) • Ic = Ic (t ) r r r • fr (Ir ) = 1+Iαr Ir , then it obviously meets ( ffr ((IIr )) − IIr )( fr (Ir ) − fr (I∗r )) ≤ 0 r ∗ ∗ c c c • The same to fc (Ic ) = 1+Iαc Ic , it also satisfies ( ffc ((IIc )) − IIc )( fc (Ic ) − fc (I∗c )) ≤ 0 c

∗

∗

The following conditions are given to illustrate the attractiveness of endemic equilibrium. (H1 )

S∗r Sc = ∗c r V∗ V∗

(3.22)

(H2 ) For the identical individual, the vaccines take the same effect of inhibition.

βcr βrc = κcr κrc

(3.23)

Therefore, there exist appropriate k1 , k2 > 0, then the followings result hold:

k1 βcr fc (I∗c )S∗r = k2 βrc fr (I∗r )S∗c

(3.24)

k1 κcr fc (I∗c )V∗r = k2 κrc fr (I∗r )V∗c

(3.25)

For this purpose, we use to the Lyapunov function to solve as this type of functions have the following advantages:

g(x ) = x − 1 − ln(x )

(3.26)

which is positive in (0, ∞) except at x = 1 which is 0. Theorem 3.4. For system (2.3), if 0 > 1 and conditions (H1 ) − (H2 ) hold, then the endemic equilibrium E∗ exists and it is globally attractive. Proof. We establish the Lyapunov function as follows:



 r S S∗r

N = k1 S∗r g

 r I I∗r

+ I∗r g



+ V∗r g

Vr V∗r





 c

+ k2 S∗c g

S S∗c

 c

+ I∗c g

I I∗c

 + V∗c g

Vc V∗c



It is clear that N(t) is non-negative in the area and arrive 0 at E∗ . For the purpose of demonstrating that negative definite, we investigate N(t) along the trajectories of Eq. (2.3) and obtain



Sr (t ) dN (t ) = k1 1 − ∗r dt S (t )



+ k2 1 −





Ir (t ) dSr (t ) + k1 1 − ∗r dt I (t )

S∗c (t ) Sc (t )











V r (t ) dIr (t ) + k1 1 − ∗r dt V (t )

Ic (t ) dSc (t ) + k2 1 − ∗c dt I (t )







dV r (t ) dt

V c (t ) dIc (t ) + k2 1 − ∗c dt V (t )



dV c (t ) dt

   Ir (t ) Ic (t ) r r r r r − βrr + β S ( t ) − ( μ + θ ) S ( t ) c 1 + αr Ir (t ) 1 + αc Ic (t )       c Sc (t ) Ir (t ) c I (t ) c c c + k2 1 − ∗c c − βrc + β S ( t ) − ( μ + θ ) S ( t ) c S (t ) 1 + αr Ir (t ) 1 + αc Ic     I∗r (t ) Ir Ic (t ) r r βr + k1 1 − r + βc Sr (t ) I (t ) 1 + αr Ir (t ) 1 + αc Ic (t )    Ir (t ) Ic (t ) r r r r r + κr + κc V (t ) − (μ + v )I (t ) 1 + αr Ir (t ) 1 + αc Ic (t )    I∗c (t ) Ir (t ) Ic (t ) c c + k2 1 − c + βc Sc (t ) βr I (t ) 1 + αr Ir (t ) 1 + αc Ic (t )    Ir (t ) Ic (t ) c c c c c + κr + κc V (t ) − (μ + v )I (t ) 1 + αr Ir (t ) 1 + αc Ic (t ) 

= k1 1 −

S∗r (t ) Sr (t )

dN (t ) is dt

10

H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956

    Ir (t ) Ic (t ) r r r r r r θ S (t ) − κr + k1 + κc V (t ) − μV (t ) 1 + αr Ir (t ) 1 + αc Ic (t )      V c (t ) Ir (t ) Ic (t ) c c c θ c Sc (t ) − κrc + k2 1 − ∗c + κ V ( t ) − μ V ( t ) c V (t ) 1 + αr Ir (t ) 1 + αc Ic (t ) 

V r (t ) 1 − ∗r V (t )

Now we simplify the equilibrium equations with E∗ and obtain:

dN (t ) dSr (t ) dSc (t ) dIr (t ) dIc (t ) dV r (t ) dV c (t ) = k1 + k2 + k1 + k2 + k1 + k2 dt dt dt dt dt dt dt S∗r dSr (t ) S∗c dSc (t ) I∗r dIc (t ) I∗c dIc (t ) V∗r dV (t ) V c dV (t ) − k1 r − k2 c − k1 r − k2 c − k1 − k2 ∗ S dt S dt I dt I dt V dt V dt = k1 r + k2 c − k1 μSr − k2 μSc − k1 (μ + vr )Ir − k2 (μ + vc )Ic − k1 μV r − k2 μV c S∗r dSr (t ) Sc dSc (t ) Ir dIr (t ) Ic dIc (t ) V r dV r (t ) V c dV c (t ) − k2 ∗c − k1 ∗r − k2 ∗c − k1 ∗c − k2 ∗c Sr dt S dt I dt I dt V dt V dt   r r r r r r S I f f ( I ) ( I ) S I r r ∗ = k1 βrr fr (I∗r )S∗r 2 + − ∗r − r − fr (I∗r ) S I∗ fr (I∗r )S∗r Ir − k1



+ k1 βcr fc (I∗c )S∗r 2 +

 + k2 βrc fr (I∗r )S∗c 2 +



Sr Ir fc (Ic )Sr I∗r f c (I c ) − ∗r − r − c fc (I∗ ) S I∗ fc (I∗c )S∗r Ir Sc Ic fr (Ir )Sc I∗c f r (I r ) − ∗c − c − r fc (I∗ ) S I∗ fr (I∗r )S∗c Ic

+ k2 βcc fc (I∗c )S∗c 2 + k2

 + k1 κrr fr (I∗r )V∗r 3 +

 + k1 κcr fc (I∗c )V∗r 3 +

 + k2 κrc fr (I∗r )V∗c 3 +

 +

κ

( ) 

k2 cc fc

+ k1 μ

S∗r

I∗c

+ k1 μ + k2 μ

Sc Ic fc (Ic )Sc I∗c f c (I c ) − ∗c − c − c fc (I∗ ) S I∗ fc (I∗c )S∗c Ic

 

fc (Ic )V r I∗r Sr Ir fc (Ic ) V∗r Sr − r r − − ∗r − r c c r r fc (I∗ ) V S∗ fc (I∗ )V∗ I S I∗ fr (Ir )V c I∗c Sc Ic fr (Ir ) V∗c Sc − c c − − ∗c − c r r c c fr (I∗ ) V S∗ fr (I∗ )V∗ I S I∗

fc (Ic )V c I∗c Sc Ic fc (Ic ) V∗c Sc 3+ − c c − − ∗c − c c c c c fc (I∗ ) V S∗ fc (I∗ )V∗ I S I∗





+ k2 μ

Vr V r Sr Sr 3 − r − ∗r r − ∗r V∗ V S∗ S

 V∗c



fr (Ir )V r I∗r Sr Ir fr (Ir ) V∗r Sr − r r − − ∗r − r r r r r fr (I∗ ) V S∗ fr (I∗ )V∗ I S I∗

Sr Sr 2 − ∗r − r S S∗

 V∗r

V∗c



Vr V c Sc Sc 3 − r − ∗c c − ∗c V∗ V S∗ S

S∗c



Sc Sc 2 − ∗c − c S S∗

  





In the above equation, the last four terms can be proven apparently to be non-positive. We only need to prove that the front square hence similar terms are non-positive as well. By adding g function, we have:

2+

Sr Ir fr (Ir )Sr I∗r f r (I r ) − ∗r − r − r fr (I∗ ) S I∗ fr (I∗r )S∗r Ir

In order to show that the above expressions are non-positive, we set:

2+

Sr Ir fr (Ir )Sr I∗r f r (I r ) − ∗r − r − fr (I∗r ) S I∗ fr (I∗r )S∗r Ir

=3−



S∗r fr (I∗r )Ir fr (Ir )Sr I∗r Ir f r (I r ) fr (I∗r )Ir − − + − r + −1 r r r r r r r S fr (I )I∗ fr (I∗ )S∗ I fr (I∗ ) I∗ fr (Ir )I∗r

Sr fr (I∗r )Ir fr (Ir )Sr I∗r = 3 − ∗r − − S fr (Ir )I∗r fr (I∗r )S∗r Ir





 r

Sr fr (I∗r )Ir fr (Ir )Sr I∗ = 3 − ∗r − − S fr (Ir )I∗r fr (I∗r )S∗r Ir



+

 +

Ir f r (I r ) − fr (I∗r ) I∗r

Ir f r (I r ) − r fr (I∗r ) I∗





fr (I∗r ) + f r (I r )



fr (I∗r ) 1− f r (I r )

f r (I r ) Ir − I∗r fr (I∗r )





H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956

since

11

S∗r fr (I∗r )Ir fr (Ir )Sr I∗r + + ≥ 3 with equality holding if and only if Sr = S∗r , fr (Ir ) = fr (I∗r ), Ir = I∗r . And then r r r S fr (I )I∗ fr (I∗r )S∗r Ir



Ir f r (I r ) − r fr (I∗r ) I∗



fr (I∗r ) 1− f r (I r )



1 = f r (I r )



Ir f r (I r ) − r fr (I∗r ) I∗



( fr (Ir ) − fr (I∗r ) ) ≤ 0

(3.27)

Therefore, it can be maintained

Sr Ir fr (Ir )Sr I∗r f r (I r ) − ∗r − r − ≤0 r fr (I∗ ) S I∗ fr (I∗r )S∗r Ir

2+

(3.28)

Besides, due to (H1 ) − (H2 ),

k1 βcr fc (I∗c )S∗r = k2 βrc fr (I∗r )S∗c

(3.29)

then, these coefficients can be combined to calculate which means



k1 β

r c fc

( ) I∗c

S∗r

Sr Ir fc (Ic )Sr I∗r f c (I c ) 2+ − ∗r − r − fc (I∗c ) S I∗ fc (I∗c )S∗r Ir



= k1 β

r c fc

( ) I∗c

S∗r



+ k2 β

c r fr

( ) I∗r

 S∗c

Sc Ic fr (Ir )Sc I∗c f r (I r ) 2+ − ∗c − c − fc (I∗r ) S I∗ fr (I∗r )S∗c Ic

Sr Ir fc (Ic )Sr I∗r Sc Ic fr (Ir )Sc I∗c f c (I c ) f r (I r ) 2+ − ∗r − r − +2+ − ∗c − c − fc (I∗c ) S I∗ fc (I∗c )S∗r Ir fc (I∗r ) S I∗ fr (I∗r )S∗c Ic





= k1 βcr fc (I∗c )S∗r  Hence, we should only prove the combined items are non-positive.



= 4+ 

f c (I c ) f r (I r ) + c fc (I∗ ) fc (I∗r )

f r (I r ) Ir − r r fr (I∗ ) I∗

Owing to

 



−

S∗r Ir fc (Ic )Sr I∗r Sc Ic fr (Ir )Sc I∗c − r − − ∗c − c − r c r r S I∗ fc (I∗ )S∗ I S I∗ fr (I∗r )S∗c Ic



fr (Ir ) − fr (I∗r ) ≤ 0, then it can be obtained

fr (I∗r )Ir f r (I r ) Ir ≤1+ r − r fr (I∗ ) I∗ fr (Ir )I∗r

(3.30)

f c (I c ) , so we can conclude fc (I∗c )

The same to

fr (I∗r )Ir fc (I∗c )Ic f r (I r ) f c (I c ) Ir Ic + ≤2+ r − + c − r c r r fr (I∗ ) fc (I∗ ) I∗ fr (I )I∗ I∗ fc (Ic )I∗c

(3.31)

Then,

≤ 4+2+ = 6− Since fr (I∗r ), fc

S∗r Sr (I c )

+

fr (I∗r )Ir fc (I∗c )Ic Sr Ir fc (Ic )Sr I∗r Sc Ic fr (Ir )Sc I∗c Ir Ic − + c − − ∗r − r − − ∗c − c − r r r c c c r r I∗ fr (I )I∗ I∗ fc (I )I∗ S I∗ fc (I∗ )S∗ I S I∗ fr (I∗r )S∗c Ic

S∗r S∗c fr (I∗r )Ir fc (I∗c )Ic fr (Ir )Sc I∗c fc (Ic )Sr I∗r − − − − − Sr Sc fr (Ir )I∗r fc (Ic )I∗c fr (I∗r )S∗c Ic fc (I∗c )S∗r Ir

fr (I∗r )Ir f ( I c )I c f ( I r )S c I c f ( I c )S r I r + fc (I∗c )Ic + f r (Ir )Sc I∗c + f c (Ic )Sr I∗r fr (Ir )I∗r c r ∗ ∗ c ∗ ∗ ∗ fc (I∗c ), Ir = I∗r , Ic = I∗c . Therefore, the first

S∗c Sc

+

≥ 6 with equality holding if and only if Sr = S∗r , Sc = S∗c , fr (Ir ) =

= four terms can be proven apparently to be non-positive. The same theory can also be applied to prove that the middle four terms are non-positive.

3+

fr (Ir )V r I∗r Sr Ir fr (Ir ) V∗r Sr − r r − − ∗r − r r r r r fr (I∗ ) V S∗ fr (I∗ )V∗ I S I∗

=4−



S∗r V r Sr fr (I∗r Ir ) fr (Ir )V r I∗r Ir fr (I∗r Ir ) fr (I∗r Ir ) − ∗r r − − + − r + −1 r r r r r r r r S V S∗ fr (I )I∗ fr (I∗ )V∗ I fr (I )I∗ I∗ fr (Ir )I∗r

= 4−

S∗r V r Sr fr (I∗r Ir ) fr (Ir )V r I∗r − ∗r r − − r r r S V S∗ fr (I )I∗ fr (I∗r )V∗r Ir



+

1 f r (I r )



Ir f r (I r ) − r r fr (I∗ ) I∗



( fr (Ir ) − fr (I∗r ) )

≤0 Thus, we only need to prove the combined items are non-positive. Owing to (H1 ) − (H2 ), then we have

k1 κcr fc (I∗c )V∗r = k2 κrc fr (I∗r )V∗c

(3.32)

12

H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956

And similar to the above calculation, the combined items are also non-positive.



fc (Ic )V r I∗r Sr Ir fc (Ic ) V∗r Sr 3+ − − − ∗r − r fc (I∗c ) V r S∗r fc (I∗c )V∗r Ir S I∗ =6+ ≤8−





fr (Ir )V c I∗c Sc Ic fr (Ir ) V∗c Sc + 3+ − − − ∗c − c fr (I∗r ) V c S∗c fr (I∗r )V∗c Ic S I∗



fc (Ic )V r I∗r Sr Ir V c Sc fr (Ir )V c I∗c Sc Ic f c (I c ) fr (Ir ) V∗r Sr + − r r − − ∗r − r − ∗c c − − ∗c − c c r c r r r c c fc (I∗ ) fr (I∗ ) V S∗ fc (I∗ )V∗ I S I∗ V S∗ fr (I∗ )V∗ I S I∗

fr (I∗r Ir ) fc (I∗c Ic ) V∗c Sc fr (Ir )V c I∗c Sc V r Sr fc (Ic )V r I∗r Sr − − c c − − ∗c − ∗r r − − ∗r fr (Ir )I∗r fc (Ic )I∗c V S∗ fr (I∗r )V∗c Ic S V S∗ fc (I∗c )V∗r Ir S

≤0 dN (t ) ≤ 0 for all {Sr (t), Sc (t), Ir (t), Ic (t), Vr (t), Vc (t) ∈ }, where  = dt dN (t ) {S0r , S0c , I0r , I0c , V0r , V0c ∈ X : S0r  ≤ S∗r , S0c  ≤ S∗c , I0r  ≤ I∗r , I0c  ≤ I∗c , V0r  ≤ V∗r , V0c  ≤ V∗c }. Consequently, ≤ 0 with equaldt ity holding only at the equilibrium E∗ . By [41], all positive solutions approach the largest invariant subset of the set dN (t ) dN (t ) ≤ 0. Since is zero only at E∗ , the largest invariant subset is just E∗ which is a singleton set. By LaSalle indt dt variant principle, endemic equilibrium E∗ is globally asymptotically stable when 0 > 1.  For t > 0, the positive definiteness of g(x) implies

4. Characterization of optimal control In this section, we calculate an optimal control model for the disease transmission dynamics by broadening the system (2.3) to contain control functions under the limited immune resources, and our objective here is to research the optimal control strategies to impose restrictions on the epidemic and maintain the cost minimum. The optimal control system of differential equations under the condition αi = 0 is as follows:

 Ir Ic r βrr + β Sr − (μ + θr u1 )Sr , c 1 + αr Ir 1 + αc Ic   dSc (t ) Ir Ic = c − βrc + βcc Sc − (μ + θc u2 )Sc , r c dt 1 + αr I 1 + αc I     dIr (t ) Ir Ic Ir Ic r r r r = βrr + β S + + κ V r − ( μ + vr )I r , κ c r c dt 1 + αr Ir 1 + αc Ic 1 + αr Ir 1 + αc Ic     dIc (t ) Ir Ic Ir Ic c c c c = βrc + β S + + κ V c − ( μ + vc )I c , κ c r c dt 1 + αr Ir 1 + αc Ic 1 + αr Ir 1 + αc Ic   dV r (t ) Ir Ic r r r  r = θ u1 S − κr + κc V r − μV r , dt 1 + αr Ir 1 + αc Ic   dV c (t ) Ir Ic c c c  c = θ u2 S − κr + κc V c − μV c . dt 1 + αr Ir 1 + αc Ic dSr (t ) = r − dt



(4.1)

with initial conditions

Sr (0 ) = S0r ,

Sc (0 ) = S0c ,

Ir (0 ) = I0r ,

Ic (0 ) = I0c ,

V r (0 ) = V0r ,

V c (0 ) = V0c , and 0 ≤ t ≤ t f .

(4.2)

In system (4.1), u1 and u2 denote the coverage rate of the Sr and Sc individuals respectively who vaccinate. tf is the final time and control functions ui (t ), i = 1, 2 are defined on the closed time interval [0, tf ]. The control functions u1 , u2 are bounded and Lebesgue integrable. The coefficients C1 , C2 , ε i , (i = 1, 2 ) denote severally the corresponding weight constants which balance the cost elements on the basis of their size and importance in the parts of the objective functional. Our goal is to reduce the number of infectious population and the cost for implementing strategy. The objective functional J(u1 , u2 ) is defined as:

J ( u1 , u2 ) =



tf 0





C1 Ir + C2 Ic + ε1 u21 + ε2 u22 dt

(4.3)

The control goal is to find optimal controls u1 and u2 to

J (u∗1 , u∗2 ) =

min J (u1 , u2 ),

(u1 ,u2 )∈

(4.4)

where the constraint is  = {(u1 , u2 ) ∈ L1 [0, t f ]| 0 ≤ u1 ≤ 1, 0 ≤ u2 ≤ 1}. The existence, uniqueness and the characteristic will be discussed in the following subsection.

H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956

13

4.1. Existence of optimal solutions Theorem 4.1. There exist optimal controls (u∗1 , u∗2 ) ∈  such that

J (u∗1 , u∗2 ) =

min J (u1 , u2 )

(4.5)

(u1 ,u2 )∈

satisfying the control system (4.1) with the primary condition (4.2). Proof. This result will be yielded by the general discussions in [29] and the following Lemma 4.1.



Lemma 4.1. The control system (4.1) with the initial condition (4.2) satisfy the following conditions. 1.  is not empty and the interval [0,1] is closed and convex. − → − → − → 2. Let X = [Sr , Sc , Ir , Ic , V r , V c ]T , U = [u1 , u2 ]T . Then the right-hand function of the control system (4.1) has the form of ( X ) + − → − → ( X ) U , where

⎛

 r − βrr



Ir Ic + βcr S r − μS r 1 + αr Ir 1 + αc Ic

⎜ ⎜   ⎜ Ir Ic ⎜ c − βrc + βcc S c − μS c r c ⎜ 1 + αr I 1 + αc I ⎜    ⎜ Ic Ir Ic ⎜ β r Ir ⎜ r 1 + α Ir + βcr 1 + α Ic Sr + κrr 1 + α Ir + κcr 1 + α Ic V r − (μ + vr )Ir r c r c ⎜ − →   ( X ) = ⎜  c r c ⎜ c Ir I I I c c c c ⎜ βr + βc S + κr + κc V c − ( μ + vc )I c ⎜ 1 + αr Ir 1 + αc Ic 1 + αr Ir 1 + αc Ic ⎜   ⎜ Ir Ic ⎜ r r − κr + κc V r − μV r ⎜ 1 + αr Ir 1 + αc Ic ⎜   ⎜ ⎝ Ir Ic c c − κr + κc V c − μV c 1 + αr Ir 1 + αc Ic ⎛  r −θ r S ⎜ 0 ⎜ − → ⎜ 0 ( X ) = ⎜ ⎜ 0 ⎝ θr Sr

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

⎞

0 −θc Sc ⎟ 0 0 0 θc Sc

0

⎟ ⎟ ⎟ ⎟ ⎠

− → − → − → − → − → − → 3. L( X , U ) = C1 Ir + C2 Ic + ε1 u21 + ε2 u22 is a convex function of U and satisfies L( X , U ) ≥ min{ε1 , ε2 } U 22 . Proof. It is clear that 1 and 2 are trivial. The condition 3 comes from the positivity of the state variables.



4.2. Characterization of optimal control The control goal is to search the optimal control function (u∗1 , u∗2 ) such that

J (u∗1 , u∗2 ) =

min J (u1 , u2 )

(4.6)

(u1 ,u2 )∈

reaches its minimum. The necessary condition for the optimal solution is given by Pontryagin’s maximum principle, which converts the optimal control (4.1) into a problem of minimizing the Hamiltonian function. For convenience, we consider the x condition that αi = 0, (i = r, c ) and then simplify (i = c, r ) into x. Then the Hamiltonian function is given by 1 + αi x

H (X, U, λ, w ) = C1 Ir + C2 Ic + ε1 u21 + ε2 u22



  βrr Ir + βcr Ic Sr − (μ + θr u1 )Sr     + λ2 c − βrc Ir + βcc Ic Sc − (μ + θc u2 )Sc    + λ3 βrr Ir + βcr Ic Sr + (κrr Ir + κcr Ic )V r − (μ + vr )Ir    + λ4 βrc Ir + βcc Ic Sc + (κrc Ir + κcc Ic )V c − (μ + vc )Ic + λ1 r −



14

H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956









+ λ5 θr u1 Sr − (κrr Ir + κcr Ic )V r − μV r

+ λ6 θc u2 Sc − (κrc Ir + κcc Ic )V c − μV c + λ7 (u1 + w21 ) + λ8 (u2 + w22 )

where λ = (λ1 , · · · , λ8 )T is the adjoint variable and w = (w1 , w2 )T is the relaxing variable. By Pontryagin’s maximum principle in [30] and the existence result for the optimal control, the following theorem is formulated. Theorem 4.2. There exists an optimal control (u∗1 , u∗2 ) and corresponding solutions Sr , Sc , Ir , Ic , Vr , Vc which minimizes J(u1 , u2 ) over  . Furthermore, there exist adjoint functions λ1 (t), , λ8 (t) such that the adjoint equations:

dλ8 (t ) dλ1 (t ) ∂ H dλ2 (t ) ∂H ∂H =− r, = − c, ··· =− , dt ∂S dt ∂S dt ∂ w2

(4.7)

with the transversality conditions

λi (t f ) = 0 i = 1, · · · , 8.

(4.8)

Then, the characteristic values u∗1 (t ) and u∗2 (t ) are given by

⎧ ⎪ (λ1 − λ5 )θr Sr ⎪ or 0, ⎨u∗1 (t ) = 2ε1 c c ⎪ ⎪ ⎩u∗2 (t ) = (λ2 − λ6 )θ S or 0. 2ε2

(4.9)

Remark 4.1. For the reader’s convenience, the specific form of







dλ1 (t ) (λ1 − λ5 )θr Sr = λ1 (t ) μ + θr max 0, min 1, dt 2ε1 −λ3 (t )(βcr Ic + βrr Ir ),

dλ1 (t ) dλ8 (t ) ,··· , are given as follows: dt dt

 +β

r c cI +

β

r r rI



  (λ − λ )θc Sc dλ2 (t ) 2 6 = λ2 (t ) μ + θc max 0, min 1, dt 2ε2 −λ4 (t )(βcc Ic + βrc Ir ),

 + βcc Ic + βrc Ir









(λ1 − λ5 )θr Sr − θr λ5 (t ) max 0, min 1, 2ε1

− θc λ6 (t ) max 0, min 1,

(λ2 − λ6 )θc Sc 2ε2

dλ3 (t ) = λ3 (t )(μ + vr − βrr Sr − κrrV r ) − C1 − λ4 (t )(βrc Sc + κrcV c ) + βrc λ2 (t )Sc + βrr λ1 (t )Sr + κrc λ6 (t )Sc + κrr λ5 (t )V r , dt dλ4 (t ) = λ4 (t )(μ + vc − βcc Sc − κccV c ) − C2 − λ3 (t )(βcr Sr + κcrV r ) + βcc λ2 (t )Sc + βcr λ1 (t )Sr + κcc λ6 (t )V c + κcr λ5 (t )V r , dt dλ5 (t ) = λ5 (t )(μ + κcr Ic + κrr Ir ) − λ3 (t )(κcr Ic + κrr Ir ), dt dλ6 (t ) = λ6 (t )(μ + κcc Ic + κrc Ir ) − λ4 (t )(κcc Ic + κrc Ir ) dt dλ7 (t ) = −2λ7 (t )w1 , dt dλ8 (t ) = −2λ8 (t )w2 . dt (4.10)

H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956

15

Remark 4.2. For the reader’s convenience, the specific form of the optimality system will be given.

⎧ r    (λ − λ )θr Sr  dS (t ) 1 5 ⎪ ⎪ = r − (βrr Ir + βcr Ic )Sr − μ + θr max 0, min 1, Sr ⎪ dt 2 ε ⎪ 1 ⎪ ⎪ ⎪    (λ − λ )θc Sc  ⎪ ⎪ dSc (t ) 2 6 ⎪ ⎪ = c − (βrc Ir + βcc Ic )Sc − μ + θc max 0, min 1, Sc ⎪ ⎪ dt 2ε2 ⎪ ⎪ ⎪ ⎪ dIr (t ) ⎪ ⎪ ⎪ = (βrr Ir + βcr Ic )Sr + (κrr Ir + κcr Ic )V r − (μ + vr )Ir ⎪ ⎪ dt ⎪ ⎪ ⎪ ⎪ ⎪ dIc (t ) ⎪ = (βrc Ir + βcc Ic )Sc + (κrc Ir + κcc Ic )V c − (μ + vc )Ic ⎪ ⎪ dt ⎪ ⎪ ⎪ ⎪   (λ − λ )θr Sr ⎪ ⎪ dV r (t ) 1 5 ⎪ ⎪ = θr max 0, min 1, Sr − (κrr Ir + κcr Ic )V r − μV r ⎪ ⎪ dt 2ε1 ⎪ ⎪ ⎪ ⎪   c c c ⎪ ⎪ ⎪ dV (t ) = θc max 0, min 1, (λ2 − λ6 )θ S Sc − (κrc Ir + κcc Ic )V c − μV c ⎪ ⎪ dt 2 ε ⎪ 2 ⎪ ⎪   ⎪   ⎪ r r ⎪ dλ1 (t ) ⎪ r c r r r max 0, min 1, (λ1 − λ5 )θ S ⎪ = λ ( t ) + θ + β I + β I μ 1 c r ⎪ ⎪ dt 2ε1 ⎪ ⎪   ⎪ r  r ⎪ (λ1 − λ5 )θ S ⎨ −θr λ5 (t ) max 0, min 1, − λ3 (t )(βcr Ic + βrr Ir ), 2ε1 ⎪   ⎪   (λ − λ )θc Sc ⎪ dλ2 (t ) ⎪ 2 6 c c c r  c ⎪ = λ2 (t ) μ + θ max 0, min 1, + βc I + βr I ⎪ ⎪ dt 2ε2 ⎪ ⎪   ⎪ ⎪ c c ⎪ ⎪ −θc λ6 (t ) max 0, min 1, (λ2 −2λε62)θ S − λ4 (t )(βcc Ic + βrc Ir ), ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ dλ3 (t ) ⎪ ⎪ = λ3 (t )(μ + vr − βrr Sr − κrrV r ) −C1 − λ4 (t )(βrc Sc + κrcV c ) + βrc λ2 (t )Sc + βrr λ1 (t )Sr + κrc λ6 (t )Sc + κrr λ5 (t )V r , ⎪ ⎪ dt ⎪ ⎪ ⎪ ⎪ ⎪ dλ4 (t ) ⎪ ⎪ = λ4 (t )(μ + vc − βcc Sc − κccV c ) −C2 − λ3 (t )(βcr Sr + κcrV r ) + βcc λ2 (t )Sc + βcr λ1 (t )Sr + κcc λ6 (t )V c + κcr λ5 (t )V r , ⎪ ⎪ dt ⎪ ⎪ ⎪ ⎪ ⎪ dλ5 (t ) ⎪ ⎪ = λ5 (t )(μ + κcr Ic + κrr Ir ) − λ3 (t )(κcr Ic + κrr Ir ), ⎪ ⎪ dt ⎪ ⎪ ⎪ ⎪ ⎪ dλ6 (t ) ⎪ ⎪ = λ6 (t )(μ + κcc Ic + κrc Ir ) − λ4 (t )(κcc Ic + κrc Ir ) ⎪ ⎪ dt ⎪ ⎪ ⎪ ⎪ ⎪ dλ7 (t ) ⎪ = −2λ7 (t )w1 , ⎪ ⎪ dt ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ dλ8 (t ) = −2λ8 (t )w2 . dt

(4.11) Proof. In view of Pontryagin’s maximum principle, we only need to consider the extremum condition ui which implies that

⎧ ∂H   r r ⎪ ⎨u1 ∂ u = 2ε1 u1 − λ1 θr S + λ5 θr S u1 = 0, 1

⎪ ⎩u ∂ H = 2ε u − λ θc Sc + λ θc Sc u = 0. 2 2 2 2 6 2 ∂ u2

∂H = 0, i = 1, 2, ∂ ui

(4.12)

Namely

⎧ ⎪ (λ1 − λ5 )θr Sr ⎪ or 0, ⎨u∗1 (t ) = 2ε1 c c ⎪ ⎪ ⎩u∗2 (t ) = (λ2 − λ6 )θ S or 0. 2ε2

(4.13)

Through the range of (u∗1 (t ), u∗2 (t )), we can obtain the properties of (u∗1 (t ), u∗2 (t )). Furthermore, the second derivative of the optimal control ui (i = 1, 2 ) is positive, which means that the optimal problem arrives minimum at controls ui (i = 1, 2 ).

16

H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956

Then, combined with the adjoint equations, the state equations, the initial and transversality conditions, the optimal system can be formulated as Remark 4.2.  4.3. Uniqueness of the optimal control In this subsection, we consider the uniqueness of the optimal control problem on the state and adjoint variables. From the bounds of the state variables, the adjoint system has bounded coefficients and it is linear in every adjoint variable. Therefore, the solutions Sr , Sc , Ir , Ic , Vr , Vc are bounded and the adjoint Eq. (4.10) is a bounded adjoint system as well which is dependent on the final time tf . Hence, for all coefficients and bounds on all solution variables, there exists a positive constant D such that |λi | < Dtf on [0, tf ]. Theorem 4.3. For tf sufficiently small, the solution to the optimality is unique. Proof. We suppose that (Sr , Sc , Ir , Ic , Vr , Vc ) and (S¯r , S¯c , I¯r , I¯c , V¯r , V¯c ) are two distinct solutions to the optimality system (4.11). For m0 > 0, we denote

Sr = em0 t g, Sc = em0 t h, Ir = em0 t j, Ic = em0 t k, V r = em0 t p, V c = em0 t n λ1 = e−m0 t q, λ2 = e−m0 t s, λ3 = e−m0 t w, λ4 = e−m0 t y, λ5 = e−m0 t z, λ6 = e−m0 t v

(4.14)

Sr = em0 t g, S¯c = em0 t h¯ , I¯r = em0 t j¯, I¯c = em0 t k¯ , V r = em0 t p, V c = em0 t n λ1 = e−m0 t q, λ2 = e−m0 t s, λ3 = e−m0 t w, λ4 = e−m0 t y, λ5 = e−m0 t z, λ6 = e−m0 t v

(4.15)

and









u∗1 = max 0, min 1, u∗1 = max 0, min 1,









(λ1 − λ5 )θr Sr 2ε1

, u∗2 = max 0, min 1,

(λ1 − λ5 )θr Sr 2ε1

, u∗2 = max 0, min 1,

(λ2 − λ6 )θc Sc 2ε2

(4.16)

(λ2 − λ6 )θc Sc 2ε2

Substitute Sr = em0 t g and λ1 = e−m0 t q into the first and the sixth differential equation of the optimality system (4.11), we obtain that



m0 e

m0 t

g+e

m0 t 

g = − (β r

r m0 t re

j+β

r m0 t k ce

)e

m0 t

g−

  (λ − λ )θr Sr 5 μ + θr max 0, min 1, 1 2ε1

Multiply both sides by e−m0 t and simplify the above equation,

m0 g + g = where  =

d . dt

e

r −m0 t

− (β j + β r r

r ck

)e

m0 t

   (λ − λ )θr Sr 1 5 g − μ + θr max 0, min 1, 2ε1



)em0 t g

 )g

(4.17)

The same to the sixth equation of the optimality system(4.11) and we can derive that

−m0 q + q =

   (λ − λ )θr Sr  5 q βrr em0 t j + βcr em0 t k + μ + θr max 0, min 1, 1 2ε1   (λ − λ )θr Sr 1 5 −(βrr em0 t j + βcr em0 t k )w − z θr max 0, min 1, , 2ε1

(4.18)

Next we subtract the above two equations corresponding to g and g, q and q and multiply each equation by appropriate difference of functions and integrate from 0 to tf . It yields from the g − g equation after multiplying by g − g that



tf

m0 0

1 2

1 2

(g − g)2 dt − [g(0 ) − g(0 )]2 + [g(t f ) − g(t f )]2 

= −βrr em0 t −θr



tf



0

tf 0

( jg − jg)(g − g)dt − βcr em0 t





max 0, min 1,



tf

0

(kg − kg)(g − g)dt − μ



(λ1 − λ5 )θr Sr 2ε1





g − θr max 0, min 1,

tf

0

(g − g)2 dt

(λ1 − λ5 )θr Sr 2ε1

(4.19)



g (g − g)dt.

Then we estimate several terms for acquiring the uniqueness result. Since em0 t < em0 t f for 0 < t < tf , we have



m0 0

tf

1 2

1 2

(g − g)2 dt − [g(0 ) − g(0 )]2 + [g(t f ) − g(t f )]2

= −βrr em0 t f −θr

 0

tf





tf 0

( jg − jg)(g − g)dt − βcr em0 t f





max 0, min 1,

(λ1 − λ5 )θr Sr 2ε1



0

tf

(kg − kg)(g − g)dt − μ





g − θr max 0, min 1,

 0

tf

(g − g)2 dt

(λ1 − λ5 )θr Sr 2ε1

(4.20)



g (g − g)dt.

H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956

17

Table 2 Parameters and their estimated values. Parameter

 c μ r

395.44 247.15 0.00029 [0,0.2] [0,0.2] 8.44 × 10−6 7.61 × 10−6 4.62 × 10−8 2.1 × 10−7 7.17 × 10−6 6.47 × 10−6 3.93 × 10−8 1.785 × 10−7

vr vc

βrr βcr βrc βcc κrr κcr κrc κcc

 tf 0

value

unit 3

10 /week 103 /week −1 week

Reference [45,46] [45,46] [46] [38] [38] [45,47] [45,47] [45,47] [45,47] [45,47] [45,47] [45,47] [45,47]

( jg − jg)(g − g)dt should be particularly analyzed. To find the estimate, Cauchy’s inequality is employed to divide the

linear terms into quadratic terms. Also, we recognize that jg − jg = jg − jg + jg − jg. As a result,



tf

( jg − jg)(g − g)dt 

0

= =

tf

0 t f 0

M1 = 2

( jg − jg + jg − jg)(g − g)dt ( j − j ) g( g − g ) + 

tf 0



tf 0

j (g − g)2 dt

M + 2 M3 ( j − j ) dt + 1 2



2

0

tf

(4.21)

(g − g)2 dt

where M1 , M3 are the upper bounds for g and j respectively. To complete the uniqueness proof of the optimal control, we imply the method in [26,31] and obtain as follow

(m0 − D − Cem0 t )



0

tf

[ ( g − g )2 + ( h − h )2 + ( j − j )2 + ( k − k )2 + ( p − p )2 + ( n − n )2

(4.22)

+(q − q )2 + (s − s )2 + (w − w )2 + (y − y )2 + (z − z )2 + (v − v )2 ]dt ≤ 0 where D and C depend on all coefficients and bounds on all solution variables. According to Theorem 3.2, C and D are independent of (Sr , Sc , Ir , Ic , Vr , Vc ) and (Sr , Sc , Ir , Ic , V r , V c ). Therefore, a contradiction comes from m0 − D − Cem0 t f > 0 for 1 m0 = Ce + D and t f < . Based on the uniqueness of the state and the adjoint solutions the optimal control is thus Ce + D unique. The proof is completed.  5. Control under limited vaccination resources In order to verify the predictability and the control rate of the model, the numerical simulation will be applied to deliberate the inoculation efficiency. The situation will be considered in which the vaccinated individuals are infected after vaccination to determine the coverage rate of the relationship between the ultimate infected individuals and the susceptible individuals (ηi , i = r, c ) under the limited vaccination resources. We first draw the impact to the ultimate infected individuals which ηr , ηc change alone. Under the Assumptions (2.1) (where κ ij = 0, i, j = r, c), we consider the relationship between the ultimate infected individuals and the change of ηr , ηr which is bounded. We take parameter values as in Table 2. 5.1. Disease control with constant coverage strategies In this paper, we aim to investigate the vaccine implementation policy of an infectious disease in a resource constrained environment. When coverage rate ηc is constant and one susceptible person occupies an unit of immune resources, thus the immune resources (s) can be denoted by



s=

0

tf

(ηr Sr + ηc Sc )dt

(5.1)

Some numerical simulation results will be presented. Fig. 6.2 shows the ultimate infected individuals of model (2.3) with constant coverage rate ηc in different resources s. From Fig. 6.2, it can be expressed that when the total resource is less than 55, bigger ηc is not better. ηc has the best choice to minimize the ultimately infected individuals and occupy resource

18

H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956 100

90 X 0.3 Y 75.65

The ultimate infected individuals

80

X 0.35 Y 64.71

70

60 X 0.4 Y 45.61

50

X 0.55 Y 39.31 40 X 0.85 Y 26.08

resource 1 = 25 resource 2 = 27.4747 resource 3 = 35 resource 4 = 45 resource 5 = 55 resource 6= 65

30

20

X1 Y 16.01

10 0

0.1

0.2

0.3

0.4

0.5

Constant vaccination rate

0.6

0.7

0.8

0.9

1

c

Fig. 6.2. Under the different resources, the ultimate infected individuals based on different constant coverage rate ηc in model (2.3).

at the same time. However, this theory is not appropriate for the abundant resource. If the total resource is more than 55, such as 65, ηc = 1 makes the ultimate infected individuals least. At this time, ηc is as large as possible. These results agree with the discussion shown by Ashrafur and Zou in [1]. When the immune resource (s) is relative enough, it is surprising that vaccinating to risky individuals does not bring any benefit to reduce the ultimately infected individuals as the critical susceptible individuals in this group are more vulnerable to infection. In conclusion, the value of ηc is significant and it can not be too big or small for the control of the disease. 5.2. Disease control with optimal coverage strategies As a vital role in the modeling of actual situations, ordinary differential systems with state or control variables are applied in various areas. Many papers have been settled down to optimal control problems and the derivation of necessary optimality conditions. In papers by Thalya and Jon [44], there were systems of differential equations considered that consisted of a combination of factors which could be analyzed numerically to confirm the solution. In this section, some numerical results are presented with the purpose of recommending that adding optimal controls into a SVIR model can have an intensive influence on the solution. We use the gradient method to solve this problem. Gradient method is a direct method and can be applied in many ways. In this method, we first pre-propose a control function u(t), which may not satisfy the necessary condition that H tends to the minimum. Then, according to the direction which the gradient decreases of H, u(t) can be improved until it meets the necessary condition. There are two systems of differential equations, the first system being the state equations involving the control, and the second being the adjoint equations (λ). Calculation steps are listed below: • Step 1. An initial guess was made for the λ gives and an initial guess (u1 , u2 ) = (u10 , u20 ) for the control over the interval. • Step 2. From here, the state equations were solved using the initial conditions

(Sr , Sc , Ir , Ic , V r , V c ) = (S0r , S0c , I0r , I0c , V0r , V0c ) and the values for (u1 , u2 ), solve (Sr , Sc , Ir , Ic , Vr , Vc ) in system (4.1) forward in time according to its differential equation in the optimality system on interval [0, tf ]. • Step 3. Using the transversality conditions (λ1 (t f ), · · · , λ8 (t f )) = (0, 0, · · · , 0 ) and the values for (u1 , u2 ) as well as (Sr (t), Sc (t), Ir (t), Ic (t), Vr (t), Vc (t)), solve (Sr (t), Sc (t), Ir (t), Ic (t), Vr (t), Vc (t)) in (4.6) backward in time according to its differential equation in the optimality system on interval [0, tf ].Since the adjoint equations depend upon the state equations, the adjoint equations use the updated state variables to determine new solutions for the adjoints. • Step 4. Then, update the solutions to which a new control is formulated. Update (u1 , u2 ) in (4.13) by entering the new (Sr (t), Sc (t), Ir (t), Ic (t), Vr (t), Vc (t)) and (λ1 (tf ), , λ8 (tf )) values into the characterization of the optimal control. • Step 5. Check convergence. The process continues until the difference in the current and previous values for the states, adjoints, and control are within an acceptable error range. If values of the variables in this iteration and the last iteration are negligibly close, output the current values as solutions. If values are not close, return to Step 2.

H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956

70

19

40 30 days

30 days

35

60

30 50 25

Ir (t)

Ic(t)

40 20

30 15 20 10 10

5

0

0 0

10

20

30

0

10

Terminal time (days)

20

30

Terminal time (days)

Fig. 6.3. Under the optimal control, the change of the ultimate infected individuals Ir (t), Ic (t) about time during 30 days in model (4.1).

10 4

4

10 4

2.5 30 days

30 days

3.5 2 3

1.5

V c(t)

V r (t)

2.5

2

1

1.5

1 0.5 0.5

0

0 0

10

20

Terminal time (days)

30

0

10

20

30

Terminal time (days)

Fig. 6.4. Under the optimal control, the change of the vaccinated individuals Vr (t), Vc (t) about time during 30 days in model (4.1).

The time frame used for this algorithm is 8 and 30 days respectively. We consider several scenarios below. In our examples, the parameter values are used in Table 2. The bounds on the controls are derived from the estimation on the effectiveness of the vaccination strategy. For the control u1 (t) and u2 (t), which mean coverage rate, it is reasonable to assume that its upper bound is 1.0 and lower bound is 0. The initial values for the states are given by Sr (0 ) = 39544, Sc (0 ) = 24715, Ir (0 ) = 60, Ic (0 ) = 35, V r (0 ) = 5, V c (0 ) = 1, respectively. The weights in the objective function are C1 = 1, C2 = 1, ε1 = 1.6, ε2 = 1.6,

20

H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956 1 u1 u2

0.9 0.8

control input

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0

5

10

15

20

25

30

Terminal time (days) Fig. 6.5. The change of the optimal control u1 and u2 about time during 30 days in model (4.1).

70

38 8 days

65

36

60

34

55

32

Ic(t)

Ir (t)

8 days

50

30

45

28

40

26

35

24

30

22 0

2

4

6

Terminal time (days)

8

0

2

4

6

8

Terminal time (days)

Fig. 6.6. Under the optimal control, the change of the ultimate infected individuals Ir (t), Ic (t) about time during 8 days in model (4.1).

which means that the minimization of the number of infections and cost has different importance. Figs. 6.3–6.8 present the time-varying optimal coverage control strategy in order to reduce the cost of infected quantity and vaccinations. Figs. 6.3-6.5 show trajectories for infectious individuals, vaccinated individuals, and control variate over time. The ultimate infected individuals initially shows a short increase to be followed by a sharp reduction. Overall, the risky infected individuals Ir remain at a low level, which is much less than their initial values, and Ic stay in a controlled range, which can be eliminated by isolation. The vaccinated population steadily increased within the first 5 days followed by a flat level. The combination of optimal control is efficacious in reducing the infectious population and promoting the vaccinated population to increase. The results in Figs. 6.3-6.4 suggest that this control strategy is sufficient in controlling the spread of disease. So far, we have been considering time horizons of 8 days for the control. In the absence of control, the infectious population of Ir and Ic rise sharply at first and decrease afterward (Fig. 6.6). Fig. 6.7 suggests that the combination of optimal control is efficacious in vaccinating the individuals. Thus, for short time scale in Figs. 6.6-6.8, if the immune resource is limited, the strategy which uses the lowest funds to minimize the infectious individuals suffice as a best control strategy in controlling disease.

H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956

10 4

4

21

10 4

2.5 8 days

8 days

3.5 2 3

1.5

V c(t)

V r (t)

2.5

2

1

1.5

1 0.5 0.5

0

0 0

2

4

6

8

0

Terminal time (days)

2

4

6

8

Terminal time (days)

Fig. 6.7. Under the optimal control, the change of vaccinated individuals Vr (t), Vc (t) about time during 8 days in model (4.1).

1 u1 u2

0.9 0.8

control input

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0

1

2

3

4

5

6

7

8

Terminal time (days)

Fig. 6.8. The change of the optimal control u1 and u2 about time during 8 days in model (4.1).

5.3. Contrast of constant and optimal coverage strategies Fig. 6.9 emerges that the ultimate infected individuals under the optimal control and the constant controls respectively. From this figure, it can be concluded that the optimal control is effective. Under the optimal control, the ultimate number of infected individuals number approximately 34.87% of that with the constant control of the total resource s = 65. Further, the optimal control only needs the total resource of s = 27.4747 which comparatively speaking, saves least about 57.73% of consumption. On the other hand, if the immune resource is assured, the ultimate number of infected individuals based on different control strategies are also demonstrated in Fig. 6.9. Under the same resource, the ultimate infected individuals based on constant is about 64.71, while that is approximate to 5.582 with the optimal control. The optimal control strategy drastically reduces the infected population about 91.37% and plays a significant role in containing the disease outbreak Figure 6.10. Thus, it can be concluded that the optimal control strategy has a positive effect on the reduction of disease and seems more costly.

22

H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956 100

90 X 0.3 Y 75.65

The ultimate infected individuals

80

X 0.35 Y 64.71

70

60 X 0.4 Y 45.61

50

X 0.55 Y 39.31 40 X 0.85 Y 26.08

resource 1 = 25 resource 2 = 27.4747 resource 3 = 35 resource 4 = 45 resource 5 = 55 resource 6 = 65 under the optimal control

30

20

10

X1 Y 16.01 X 0.6 Y 5.582

0 0

0.1

0.2

0.3

0.4

0.5

Variant vaccination rate

0.6

0.7

0.8

0.9

1

c

Fig. 6.9. The ultimate number of infected individuals based on constant coverage rate ηc strategies under different immune resource in model (2.3) VS that based on the optimal control strategy with immune resource s = 27.4747 in model (4.1). (a) Under the optimal control, the number of the ultimate infected individuals is calculated as shown in the very bottom red line; (b) The corresponding constant controls are chosen as the other colors and notations; (c) The coordinates in different color and notation lines present the least number of the ultimate infected individuals under appropriate ηc . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6.10. Under the same resources, the ultimate infected individuals based on constant coverage rate ηc and the optimal control in model (4.1).

6. Conclusions and discussion In this paper, we focus on a modified epidemic disease SVIR model with different groups incorporating an ordered classification of vaccinations status. Then, a mathematical model of sequence is established, which contains partition vaccination strategies. For the underlying system, we gave the sections to give a rigorous series of analyses in the above sections, such as showing the positive invariance of solutions, the basic reproductive number, the persistence of disease, the local stability,

H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956

23

and the global stability of equilibrium. It is worth mentioning that, by constructing a Lyapunov function, we prove the global attractor to be a set of single points. The global dynamics are totally determined by the 0 (threshold).When 0 < 1 is the disease-free equilibrium, the disease-free equilibrium is globally asymptotically stable, and disease is tends toward extinction. When 0 > 1, the local disease equilibrium point is globally asymptotically stable, and the disease will persist at that point. To control the spread of disease, the optimal regulations of an extended model with vaccination has been proposed, and optimal control strategies for several scenarios have been identified. We concentrate on the application of optimal control theory to minimize the spread of infection. The optimality is measured by the minimality of the infectious burden and cost. Then, we change the original problem into the Hamiltonian function and find the extreme value. Finally, we derive the necessary conditions for the optimal control problem based on the Pontryagin’s maximum principle. The uniqueness and existence of the optimality system guarantee that the optimal control both exists and is unique in the depicted cases. Furthermore, numerical simulations of the results obtained from the optimal control technique have some interesting implications. With enough large immune resources, the disease in system (4.1) will be always restrained within limits, tending to elimination. The terminal time largely influenced the optimal control largely and the longer time disease sustained in the population, the better the effect of vaccination. In detail, for a disease with a shorter infectious period, there was a higher cost to controlling the disease in order to achieve the optimal control objective, because it was optimal to start treatment with maximal effort at the start of the epidemic and to continue treatment with maximal effort until the switch time had arrived. Otherwise, the lower total cost needs to be expended to control the disease with a longer infectious period. Next, in the situation with limited-resources, vaccinations to the critical group cannot have positive effect on the reduction of disease. Vaccinations of the ‘risky group only’ strategy will be better policy. Besides, numerical experiments also simulate the following disease control strategies to demonstrate values for the optimal control: (a)disease control with constant coverage strategies under limited resources; (b) disease control with optimal coverage strategies under limited resources. In the end, for our system (4.1), minimizing the cost of disease and the ultimate number of infected individuals can only be achieved by the optimal control strategy. Therefore, the optimal vaccination strategy is economical and effective in reducing disease transmission. Our results have potential applications in practical situations on the basis of more realistic data in [40,47-49]. In comparison with our approach, the viewpoint of Ashrafur and Zou [1] for controlling the disease is to reduce the basic reproduction number R0 and only implement the constant control. Thus, it follows that our work seems more precise and enriches further their results. The model can also be improved by adding several realistic aspects. For instance, to express the result of immunization, we can incorporate the delay and the impact of passive immunity. We leave these as possible future research projects. Contributions S.L. conceived and designed the study, H.Z. and S.L. wrote the first draft of the paper. K.P. developed the model framework; H.Z. and Z.Y. performed the experiments and analyzed the outcomes; K.P., H.Z. and S.L. contributed to the revisions and writing of the final draft, and all authors approved its content. Declaration of Competing Interest None. Acknowledgments The authors would like to thank Professors Donald L. DeAngelis, Yuming Chen, Jinliang Wang and Dr. Wenqiang Ji for their valuable suggestions and comments. H Zhang and S. Liu are supported by the Natural Science Foundation of China (NSFC) (No. 11871179, 11771374). Z Yang is supported by NSFC (No. 11501568, 11871179). References [1] S.M. Ashrafur, X. Zou, Modelling the impact of vaccination on infectious diseases dynamics, J. Biol. Dyn. 9 (2015) 307–320, doi:10.1080/17513758.2014. 986545. [2] E. Belongia, A.L. Naleway, Smallpox vaccine: the good, the bad, and the ugly, Clin. Med. Res. 1 (2) (2003) 87–92, doi:10.3121/cmr.1.2.87. [3] A. Gumel, S. Ruan, T. Day, J. Watmough, et al., Modeling strategies for controlling SARS outbreaks based on Toronto, Hong Kong, Singapore and Beijing experience, Proc. R. Soc. Lond. B 271 (2004) 2223–2232, doi:10.1098/rspb.2004.2800. ´ , M. Pawlowska, Hepatitis B virus serologic markers and anti-hepatitis B vaccination in patients with diabetes, Med. Sci. Monit. [4] W. Halota, M. Muszynka 8 (2002) 516–519. [5] H.W. Hethcote, The mathematics of infectious diseases, SIAM Rev. 42 (4) (20 0 0) 599–653, doi:10.1137/S0036144500371907. [6] W.O. Kermack, A.G. McKendrick, A contribution to the mathematical theory of epidemics, Proc. R. Soc. Lond. A 115 (772) (1927) 700–721. [7] S. Lakhani, Early clinical pathologists: edward jenner, J. Clin. Pathol. 45 (9) (1992) 756–758, doi:10.1136/jcp.44.8.621. [8] I.M. Longini, M.E. Halloran, Strategy for distribution of influenza vaccine to high-risk groups and children, Am. J. Epidemiol. 161 (4) (2005) 303–306, doi:10.1093/aje/kwi053. [9] S.A. Plotkin, Correlates of vaccine-induced immunity, Clin. Infect. Dis. 47 (3) (2008) 401–409, doi:10.1086/589862. [10] H.M. Foy, M.K. Cooney, I. Allan, Longitudinal studies of types A and B influenza among Seattle school children and families, J. Infect. Dis. 134 (4) (1976) 362–369, doi:10.1093/infdis/134.4.362.

24

H. Zhang, Z. Yang and K.A. Pawelek et al. / Applied Mathematics and Computation 371 (2019) 124956

[11] S. Gandon, M. Mackinnon, S. Nee, A. Read, Imperfect vaccination: some epidemiological and evolutionary consequences, Proc. R. Soc. Lond. B 270 (2003) 1129–1136, doi:10.1098/rspb.2003.2370. [12] I.M. Longini, M.E. Halloran, Strategy for distribution of influenza vaccine to high-risk groups and children, Am. J. Epidemiol. 161 (4) (2005) 303–306, doi:10.1093/aje/kwi053. [13] M. Loeb, M.L. Russell, L. Moss, K. Fonseca, et al., Effect of influenza vaccination of children on infection rates in Hutterite communities: a randomized trial, JAMA 303 (10) (2010) 943–950, doi:10.1001/jama.2010.250. [14] C. Suppo, J.M. Naulin, M. Langlais, M. Artois, A modelling approach to vaccination and contraception programmes for rabies control in fox populations, Proc. R. Soc. Lond. B 267 (1452) (20 0 0) 1575–1582, doi:10.1098/rspb.20 0 0.1180. [15] A.R. Tuite, D.N. Fisman, J.C. Kwong, A.L. Greer, Optimal pandemic influenza vaccine allocation strategies for the canadian population, PLoS One 5 (5) (2010) e10520, doi:10.1371/journal.pone.0010520. [16] E. Simons, M. Mort, A. Dabbagh, P. Strebel, L. Wolfson, Strategic planning for measles control: using data to inform optimal vaccination strategies, JID 204 (2011) 28–34, doi:10.1093/infdis/jir095. [17] L. Cai, X. Li, J. Yu, Analysis of a delayed HIV/AIDS epidemic model with saturation incidence, J. Appl. Math. Comput. 27 (2008) 365–377, doi:10.1007/ s12190- 008- 0070- 3. [18] V. Capasso, G. Serio, A generalization of the Kermack-McKendrick deterministic epidemic model, Math. Biosci. 42 (1978) 41–61, doi:10.1016/ 0 025-5564(78)90 0 06-8. [19] J.K. Hale, P. Waltman, Persistence in infinite-dimensional systems, SIAM J. Math. Anal. 20 (1989) 388–395, doi:10.1137/0520025. [20] K.A. Pawelek, S. Tobin, C. Griffin, et al., Impact of a waning vaccine and altered behavior on the spread of influenza, AIMS Med. Sci. 4 (2) (2017) 217–232, doi:10.3934/medsci.2017.2.217. [21] J. Medlock, A.P. Galvani, Optimizing influenza vaccine distribution, Science 325 (2009) 1705–1708, doi:10.1126/science.1175570. [22] R. Xu, Z.E. Ma, Global stability of a SIR epidemic model with nonlinear incidence rate and time delay, J. Nonlinear Anal. 10 (5) (2009) 3175–3189, doi:10.1016/j.nonrwa.2008.10.013. [23] L.M. Cai, X.Z. Li, M. Ghosh, Global stability of a stage-structured epidemic model with a nonlinear incidence, Discrete Dyn. Nat. Soc. 214 (1) (2009) 73–82, doi:10.1016/j.amc.2009.03.088. [24] D. Kirschner, S. Lenhart, S. Serbin, Optimal control of the chemotherapy of HIV, Math. Biol. 35 (7) (1997) 775–792, doi:10.10 07/s0 02850 050 076. [25] J. Muller, Optimal vaccination patterns in age-structured populations: endemic case, Math. Comput. Model. 31 (20 0 0) 149–160, doi:10.1016/ s0895-7177(0 0)0 0 033-9. [26] T. Burden, J. Ernstberger, K.R. Fister, Optimal control applied to immunotherapy, J. Discrete Contin. Dyn. Syst. Ser. B 4 (1) (2004) 135–146, doi:10.3934/ dcdsb.2004.4.135. [27] Y. Yang, S. Tang, X. Ren, H. Zhao, C. Guo, Global stability and optimal control for a tuberculosis model with vaccination and treatment, Discrete. Cont. Dyn. Syst. B 21 (3) (2016) 1009–1022, doi:10.3934/dcdsb.2016.21.1009. [28] T. Malik, M. Imran, R. Jayaraman, Optimal control with multiple human papillomavirus vaccines, J. Theor. Biol. 393 (2016) 179–193, doi:10.1016/j.jtbi. 2016.01.004. [29] W.H. Fleming, R. Rishel, Deterministic and Stochastic Optimal Control, Springer-Verlag, New York, 1975. [30] F.H. Richard, P.S. Suresh, G.V. Raymond, A survey of the maximum principles for optimal control problems with state constraints, SIAM Rev. 37 (2) (1995) 181–218, doi:10.1137/1037043. [31] Y.Z. Pei, M.M. Chen, X.Y. Liang, Z.M. Xia, Y.F. Lv, C.G. Li, Optimal control problem in an epidemic disease SIS model with stages and delays, Int. J. Biomath. 9 (5) (2016) 1650072, doi:10.1142/S1793524516500728. [32] W.D. Wang, X.Q. Zhao, An epidemic model in a patchy environment, Math. Biosci. 190 (1) (2004) 97–112, doi:10.1016/j.mbs.2002.11.001. [33] M.W. Hirsch, H. Smith, X.Q. Zhao, Chain transitivity, attractivity, and strong repellors for semidynamical systems, J. Dyn. Differ. Equations 13 (1) (2001) 107–131, doi:10.1023/A:1009044515567. [34] X.Q. Zhao, Uniform persistence and periodic coexistence states in infinite-dimensional periodic semiflows with applications, Can. Appl. Math. 3 (1995) 473–495. [35] H. Smith, X.Q. Zhao, Robust persistence for semidynamical systems, Nonlinear Anal. 46 (2001) 6169–6179, doi:10.1016/S0362- 546X(01)00678- 2. [36] H.R. Thieme, Persistence under relaxed point-dissipativity (with applications to an endemic model), SIAM J. Math. Anal. 24 (1993) 407–435, doi:10. 1137/0524026. [37] P. Dreessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci. 180 (2002) 29–48, doi:10.1016/S0025-5564(02)00108-6. [38] E. Simons, M. Ferrari, J. Fricks, et al., Assessment of the 2010 global measles mortality reduction goal: results from a model of surveillance data, Lancet 379 (2012) 2173–2178, doi:10.1016/s0140- 6736(12)60522- 4. [39] O. Diekmann, J.S.P. Heesterbeek, J.A.J. Metz, On the definition and the computation of the basic reproduction ratio r0 in models for infectious diseases in heterogeneous populations, J. Math. Biol. 28 (4) (1990) 365–382, doi:10.10 07/BF0 0178324. [40] P. van den Driessche, J. Watmough, Reproduction numbers and subthreshold endemic equilibria for compartmental models of disease transmission, Math. Biosci. 180 (2002) 29–48, doi:10.1016/S0025-5564(02)00108-6. [41] J. Hale, S.V. Lunel, Introduction to Functional Differential Equations, Springer, New York, 1993. [42] World Health Organization, Hepatitis B vaccines, in: Weekly Epidemiological Record, vol. 79, 2004, p. 255. [43] Y. Muroya, Y. Enatsub, T. Kuniya, Global stability for a multi-group SIRS epidemic model with varying population sizes, Nonlinear Anal 14 (3) (2013) 1693–1704, doi:10.1016/j.nonrwa.2012.11.005. [44] T. Burden, J. Ernstberger, K.R. Fister, Optimal control applied to immunotherapy, Discrete Contin. Dyn. Syst.Ser. B 41 (12) (2004) 135–146, doi:10.1023/B: DESI.0 0 0 0 012445.67356.7a. [45] Z. Linhua, W. Yan, X. Yanyu, L. Michael, Global dynamics of a discrete age-structured SIR epidemic model with applications to measles vaccination strategies, Math. Biosci. 308 (2019) 27–37, doi:10.1016/j.mbs.2018.12.003. [46] W.H. Organization, Global health observatory (GHO) data: India, 2017, http://www.who.int/gho/countries/ind/country_profiles/en(accessed on April 16). [47] J. Mossong, N. Hens, M. Jit, et al., Social contacts and mixing patterns relevant to the spread of infectious diseases, PLoS Med. 5 (2008) E74, doi:10. 1371/journal.pmed.0 050 074. [48] S. Ruoyan, S. Junping, Global stability of multigroup epidemic model with group mixing and nonlinear incidence rates, Appl. Math. Comput. 218 (2011) 280–286, doi:10.1016/j.amc.2011.05.056. [49] W. Jinliang, W. Jing, Toshikazu Kuniya, analysis of an age-structured multi-group heroin epidemic model, Appl. Math. Comput. 347 (2019) 78–100, doi:10.1016/j.amc.2018.11.012.