Model 1: Zika Virus SIR Transmission Model

Bonyah and Okosun [17] formulated a mathematical model for transmission of Zika virus (ZIKV), which incorporated some optimal control strategies for prevention, treatment, and use of an insecticide. Pontryagin's maximum principle was used to characterize the necessary conditions for optimal control of ZIKV. The authors explored an SIR model to examine the dynamics of transmission of ZIKV to humans [40]. The use of preventive and insecticide techniques to minimize the spread of the virus showed a reduction in both infected humans and mosquitoes. The human population at time t is given by Nh (t) = SH(t) + IH(t) + Rn(t), where SH(t), IH{t), and RH(f) respectively denote the number of individuals with Zika symptoms, infectious individuals, and individuals recovered from Zika. The overall vector (mosquito) population at time t is given by Nv(t) = Sv(t) + Iv(t),

where Sv(t) and Iv(t) respectively denote the susceptible and infectious mosquitoes. The transmission is governed by the following set of ODEs [17]:

All the model parameters are positive and the initial conditions are given by

The parameters in the model are defined as follows:

/3H,Pv : Transmission rates of ZIKV from humans to mosquitoes and from the vector (mosquitoes) to humans respectively.

Hh,Hv ■ Natural death rates of host and vector respectively.

AH : Recruitment rate into susceptible population.

Т)и : Recovery rate from treatment.

Av : Recruitment rate into susceptible mosquito population. r)v : Vector death rate from the use of an insecticide.

1/у: Average infectious period for humans.

(Control functions) Prevention, treatment, and insecticide controls respectively. These control functions are bounded and Lebesgue integrable.

о : Effective contact rate between infected and susceptible humans.

Analysis of equilibrium points: The DFE point is given by E0 = (AH//UH, 0,0,Ау//Лу, 0). Let us consider the case when the control measures are not employed. In that case, the EE point is given by Ег =(S'H,Ih,Rh,Sv,Iv), where

where Гн is a positive root of the equation + ЬГН + с = 0, with

The solution given here differs from the solution given by the authors [17]. We analyzed the model using the above solutions. (In a private communication, the authors are informed of the results). We have a > 0. If c < 0, then the equation has always a positive root irrespective of the sign of b. (The signs of the coefficients are then +, +, - or +, - -.)

Calculation of the basic reproduction number 7Z0: The next-generation operator approach of Okosun and Makinde [73,85] was employed to determine the reproduction number. The matrices describing the new infections generated and the transition terms, F and V, are given by

The basic reproduction number R0 (the spectral radius of FV ') is obtained as We can write t, = + = K=^,

Цн О pv 2

'R0 depends on the infection rate of human and vector populations and is driven by the respective recruitment rates. The reduction of ZIKV is driven by the recovery rate and natural mortality of human and vector populations.

Stability analysis: The Jacobian matrix of system (5.5) at the DFE E0 is obtained as

Three of the eigenvalues of JEo are Н/- цН/- цу, and the remaining two eigenvalues are the solutions of the polynomial P(x) = x2 + r^x + = 0, where

A sufficient condition that > 0 and r0 > 0, is (2 + R3) < 1. In this case, the equation has negative roots or a complex pair with negative real parts. Then, the DFE is asymptotically stable.

Consider the Lyapunov function

Taking the time derivative of (5.7) along the solutions of the system (5.5), we obtain

Substituting the DFE solution, H = Анн and S“ = Av jnv in equation (5.8) and rearranging, we get

since 'R0 = '/?.| + 'Ro. L'(t)<0, if 'R0 < 1. L'(t) = 0, if SH = Sh,Sv = S°,IH = Iv = 0. Hence, the largest compact invariant set (SH,IHlRH,SV,IV) e £2: L'(t) = 0 is the singleton set {E0}. Thus, E0 is globally asymptotically stable under the above conditions, when all solutions of system (5.5) in R5 are bounded [65].

Stability of the endemic equilibrium: Denote t4 = S'Hl5H,t5 = Svpv, t6 = /3h( +/v), f8 = t6 + ци, t4 = nv + 1‘npv, f10 = t3 - crf4, and t7 = t9 + tw. The Jacobian matrix of the system (5.5) at the EE point E, is given by

The characteristic equation of the Jacobian matrix J£l is given by where ?7i2 = f8 + tq + ho = t7 + f8, T)u = hh + hho + — hh/

(The solution given here differs from the solution given by the authors [17]). Two of the eigenvalues of J£] are -Цн,-Цу, and the remaining three eigenvalues are the solutions of the polynomial q(x) = 0. The roots of q(x) = 0 are negative or have negative real parts when the following Routh-Hurwitz criteria are satisfied: pw >0, r/u >0, rp2 > 0,ПчЛп >0. From these conditions, we obtain the following sufficient conditions:

If these conditions are satisfied, then the EE point E] is asymptotically stable.

Global stability of the endemic equilibrium: Consider the following non-linear Lyapunov function to establish the global stability of the EE

Differentiating the above equation, simplifying, and setting a, = a2 = a3 = a4 = 1, we obtain

The equation can be written as

Sufficient conditions for the quadratic form V", to be negative definite are the following:

Thus, the EE point whenever it exists is globally asymptotically stable under the conditions (5.10).

Numerical simulations: To study the stability, numerical computations are performed using the following parameter values: Лн = 100, pH = 0.0001, (7 = 0.05, цн = 1/(365x60), y = l,000, Av = 1,000, pv = 0.00002, цУ = 0.07142. The initial conditions are taken as (100,000, 2, 2, 2,000, 10). The DFE is obtained as E0(2.19 x 106,0,0,1.4x 104,0) with 7?o = 0.932033. Local asymptotic stability of DFE is displayed in Figure 5.1. Now, if we take the values /),, = 0.0005, and pv = 0.0002, we obtain the EE as £,(50969.3,0.0976726,2.13903x 106,13992.2,3.82772) with % = 6.579073. Local asymptotic stability of EE is displayed in Figure 5.2.


Local asymptotic stability of the DFE E0.


Local asymptotic stability of the EE, E, with /?„ = 0.0005 and /lv = 0.0002.

Optimal Control Analysis

The authors [17] used the Pontryagin's maximum principle to determine the necessary conditions for the optimal control of the disease. The problem is to minimize the number of ZIKV-infected human hosts and the cost, which is involved in the mass treatment and insecticide controls. The objective function is taken as

where d],d2, and d3 denote the weighting constants for prevention, treatment, and insecticide efforts, respectively. The costs of the prevention, treatment, and insecticide are non-linear and assume quadratic function forms. Optimal controls fi,Цъ and qj are to be determined such that

The model system with control is

The necessary conditions that an optimal solution must satisfy come from the maximum principle of Pontryagin et al. [92]. This principle converts (5.11)-(5.12) into a problem of minimizing pointwise the following Hamiltonian H, with respect to ць ц2, and ц3.

where As„ ,A1h/ARh ,ASv,Xtv are the adjoint variables or co-state variables. The system of equations is found by taking the appropriate partial derivatives of the Hamiltonian (5.13) with respect to the associated state variables. We summarize the results of the authors [17]: given optimal controls з and solutions SH,IH/RH,SV,IV of the corresponding state

system (5.11)-(5.12) that minimize / (цл, ц2, /Т) over U, there exist adjoint variables As„, A(„, Ar„ , As„, Ajv satisfying the equation

As,, (f/) = A/„ (f/) = ARh [tf) = ASv, (f/) = A,v. (tf) = 0, where

Corollary 4.1 of Fleming and Rishel [45] gives the existence of an optimal control due to the convexity of the integrand of J with respect to Цг,ц2, and ц3, a priori boundedness of the state solutions, and the Lipschitz property of the state system with respect to the state variables. The differential equations governing the adjoint variables are obtained by differentiation of the Hamiltonian function, evaluated at the optimal control.

Then, the adjoint equations can be written as [45]:

Authors [17] have presented the simulation results of the model (5.12) for different cases of control mechanism including the case when all control mechanisms (щ,ц2,Цз) are used to optimize the objective function /. That is, all the prevention, treatment, and use of insecticide controls are optimized. A significant difference in the number of infected humans lH in the controlled cases and the cases without control can be observed from Figure 5.3a. The result depicted in Figure 5.3b suggests that this strategy is very efficient and effective in the control of the number of infected mosquitoes Iv. Numerical results on the optimal controls suggest that the best strategy is to combine all controls, that is, preventive, treatment, and insecticides to control the disease.


Effect of all controls on virus transmission, with (a) /?„ = 0.2 and (b) fiv = 0.09. [From Bonyah, E., Okosun, K. O. 2016. Mathematical modeling of Zika virus. Asian Рас. ]. Trap. Dis. 6(9), 673-679 [17]. Copyright 2016. Reprinted with kind permission from Elsevier.]

Model 2: Zika Virus SEIR Transmission Model

Bonyah et al. [16] modified Model 1, described in Section 5.3, by incorporating the exposed population and called it an SEIR Zika epidemic model. The authors have considered the human-to-human infection as well as the vector (mosquito)-to-human transmission. The total human population NH(t) is taken as NH(t) = SH + EH + IH+ RH, where SH(t), En(t), /„(f), and Rft(t) are the susceptible humans, exposed humans, infected humans, and recovered humans respectively. Similarly, the mosquito population Nv(t) is partitioned as Nv =SV +EV + IV, where Sv(t),Ev(t), and Iv(t) are the susceptible vector, exposed vector, and infected vector respectively. The model system is given by [16]:

All the model parameters are positive and the initial conditions are given by S„(0) > 0,E„(0) > 0, J„(0) > 0,R„(0) > 0,SV(0) > 0,EV(0) > 0,/v(0) > 0.

The parameters are defined as follows:

н, pv: Effective contact rate between susceptible humans and infected mosquitoes and the transmission rate from infected humans to susceptible vectors respectively.

pH, Pv- Natural mortality rates due to each subpopulation of human and vector groups respectively.

AH,AV: Recruitment of susceptible human and susceptible mosquito populations respectively.

p: Effective contact rate between infected humans and susceptible humans that can result in infection.

Y, if. Natural and treatment rates.

Xn,8y: Progression rates from exposed to infected human and mosquito populations respectively.

Let the total dynamics of the human population be given by Hence, N'H(t) + /uHNH H.

Integrating the above inequality, the following inequality is obtained [14]

where NH (0) = N„ (S„ (0) + £„ (0) + IH (0) + RH (0)).

As t —> <*>, we obtain 0 < NH < Анн.

The total dynamics of vector population is given by N'v (f) = Av - /uvNv.

The solution of this equation is

where bly (0)= bly (Sy (0) + Ry (0) + Iy (0)j. As t —) we obtain Ny (f) = А у / /Лу.


which is positively invariant, dissipative and the global attractor is attained in П.

Analysis of equilibrium points and the basic reproduction number^: The DFE point is given by £0 = н//ли,0,0,0, Av/ цу, 0,0). To obtain the basic reproduction number, form the matrices

where L = риАц//Лн> h = PvAy/Цу, k = Ци + кг = Цн + У + and к$ = Цу + 8у. Taking = РнАнХн, f<> = fiySyAy, the basic reproduction number of the model, which is the spectral radius of the matrix р(£У~М, is given by

We find that R2 > Ri, and Ro > 2RX.

Endemic equilibrium: The EE point of the system is given by Ex =(S’h,Eh, Гн, R’h, S’v, E‘v, I'vj, where

where a = pkxk2kipHpvPv > 0, and 'Ro = Ro + 2'R (1 - 7?o).

Note that Ro > 2R. We have a > 0. Now, c > 0 when R0 < 1, and c < 0, when Rtl > 1. The signs of b and c decide the existence of a positive solution. Consider the case c> 0, that is 'Ro < 1: (i) when b> 0, then the signs of the coefficients are +,+,+. In this case, there is no positive root, (ii) When b< 0, then the signs of the coefficients are +,-,+. In this case, there are two positive roots if b2 - 4ac > 0. That is, two positive equilibrium points exist. For 'Ro = 1, c = 0, and if b < 0, then a unique positive solution I‘H = -b/a exists. Now, consider the case c < 0, that is Ro > 1: irrespective of the sign of b, there is a positive root, that is, a unique EE point exists.

Using MATLAB 2013, we have simulated the system (5.14) taking the parameter values as %H = 0.0022= 0.0009, pH = 0.01, S„ = 0.3, juv = 0.003^,, = 1.3,Лн = 0.4,/7 = 0.11, p = 0.029,y = 0.0614799, for the two cases: (i) pH = 0.00005. The value of Ro is obtained as 'Ro = 0.505787. (ii) pH = 0.0005. The value of Ro is obtained as 'R0 = 1.59964. Initial conditions are taken as [100,2,3,1,500,5,70]. For this set of parameter values, we obtain

The time series is presented in Figure 5.4a for 'Ro < 1 and in Figure 5.4b for Ro > 1.

Stability analysis: We can obtain the conditions on the coefficients of the characteristic equation so that the DFE point E0 = (H/pH ,0,0,0, Av/pv,0,0) of the model system (5.14) is locally asymptotically stable (Problem 5.1, Exercise 5).

The endemic equilibrium point E{ of the model system (5.14) is asymptotically stable for Ro> 1 (Problem 5.2, Exercise 5 [5]). The authors have derived the Jacobian matrix without considering RH and obtained a characteristic equation, which is of order six. We have derived the Jacobian matrix considering RH and obtained a characteristic equation of order seven. The results in both the cases are the same.


Time series for the model system (5.14): (a) 7?,, < 1 and (b) > 1.

< Prev   CONTENTS   Source   Next >