# 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) = S* _{H}(t) + I_{H}(t) + Rn(t),* where

*S*and R

_{H}(t), I_{H}{t),_{H}(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

*N*

_{v}(t) = S_{v}(t) + I_{v}(t),where *S _{v}(t)* and

*I*respectively denote the susceptible and infectious mosquitoes. The transmission is governed by the following set of ODEs [17]:

_{v}(t)

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

The parameters in the model are defined as follows:

*/3 _{H},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.

A_{H} : Recruitment rate into susceptible population.

*Т) _{и}* : Recovery rate from treatment.

*A _{v}* : Recruitment rate into susceptible mosquito population.

*r)*: Vector death rate from the use of an insecticide.

_{v}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 E_{0} = *(A _{H}//U*

_{H},

*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 *7Z*_{0}: 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 *R _{0}* (the spectral radius of

*FV*') is obtained as We can write t, = +

_{=}K

_{=}^,

*Цн О **p _{v}* 2

*'R _{0}* 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 *E _{0}* is obtained as

Three of the eigenvalues of *J _{Eo}* are

*-ц*and the remaining two eigenvalues are the solutions of the polynomial

_{Н/}- ц_{Н/}- ц_{у},*P(x) = x*0, where

^{2}+ r^x + =

A sufficient condition that *>* 0 and *r _{0}* > 0, is (2 + R

_{3}) < 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, *S° _{H} = А_{н}/ц_{н}* and S“ =

*A*in equation (5.8) and rearranging, we get

_{v}jn_{v}since *'R _{0} =* '/?.| +

*'Ro*.

*L'(t)<0,*if

*'R*< 1.

_{0}*L'(t)*= 0, if S

_{H}= Sh,S

_{v}

*= S°,I*0. Hence, the largest compact invariant set

_{H}= I_{v}=*(S*e £2:

_{H},I_{Hl}R_{H},S_{V},I_{V})*L'(t) = 0*is the singleton set {E

_{0}}. Thus, E

_{0}is globally asymptotically stable under the above conditions, when all solutions of system (5.5) in

*R*are bounded [65].

^{5}Stability of the endemic equilibrium: Denote *t _{4} = S'_{H}l5_{H},t_{5} = Svpv, t_{6} = /3h(*

The characteristic equation of the Jacobian matrix J_{£l} is given by
where ?7i2 = f_{8} + *tq +* ho = *t _{7} +* f

_{8},

*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: *p _{w}* >0,

*r/u*>0,

*rp*>0. From these conditions, we obtain the following sufficient conditions:

_{2}> 0,ПчЛп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, = a _{2} = a_{3} =* a

_{4}= 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, *p _{H}* = 0.0001, (7 = 0.05,

*ц*1/(365x60), y = l,000, A

_{н}=_{v}= 1,000,

*p*= 0.00002,

_{v}*ц*0.07142. The initial conditions are taken as (100,000, 2, 2, 2,000, 10). The DFE is obtained as E

_{У}=_{0}(2.19 x 10

^{6},0,0,1.4x 10

^{4},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

*p*= 0.0002, we obtain the EE as £,(50969.3,0.0976726,2.13903x 10

_{v}^{6},13992.2,3.82772) with % = 6.579073. Local asymptotic stability of EE is displayed in Figure 5.2.

FIGURE 5.1

Local asymptotic stability of the DFE *E _{0}.*

FIGURE 5.2

Local asymptotic stability of the EE, E, with /?„ = 0.0005 and */l _{v}* = 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],d _{2},* and

*d*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

_{3}*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 A_{s}„ ,A_{1h/}A_{Rh} ,A_{Sv}*,X _{tv}* 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

*S*of the corresponding state

_{H},I_{H/}R_{H},S_{V},I_{V}system (5.11)-(5.12) that minimize / (*ц*_{л}, *ц*_{2}, /Т) *over U,* there exist adjoint variables A_{s}„, A_{(}„, Ar„ , A* _{s}„,* Aj

_{v}satisfying the equation

As,, (f/) = A/„ (f/) = A_{Rh} *[tf*) = A_{Sv}, (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

*ц*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.

_{3},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

*l*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

_{H}*I*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.

_{v}.FIGURE 5.3

Effect of all controls on virus transmission, with (a) /?„ = 0.2 and (b) *fi _{v}* = 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 *N _{H}(t)* is taken as

*N*where

_{H}(t) = S_{H}+ E_{H}+ I_{H}+ R_{H},*S*/„(f), and

_{H}(t), En(t),*Rft(t)*are the susceptible humans, exposed humans, infected humans, and recovered humans respectively. Similarly, the mosquito population

*N*is partitioned as

_{v}(t)*N*where

_{v}=S_{V}+E_{V}+ I_{V},*S*and

_{v}(t),E_{v}(t),*I*are the susceptible vector, exposed vector, and infected vector respectively. The model system is given by [16]:

_{v}(t)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,S_{V}(0) > 0,E_{V}(0) > 0,/_{v}(0) > 0.

The parameters are defined as follows:

*Iв _{н}, p_{v}:* Effective contact rate between susceptible humans and infected mosquitoes and the transmission rate from infected humans to susceptible vectors respectively.

*p _{H}, Pv-* Natural mortality rates due to each subpopulation of human and vector groups respectively.

*A _{H},A_{V}:* 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) + /u_{H}N_{H} H.*

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

where *N _{H}* (0) = N„ (S„ (0) + £„ (0) +

*I*(0) +

_{H}*R*(0)).

_{H}As *t —>* <*>, we obtain 0 < *N _{H} < А_{н}/ц_{н}.*

The total dynamics of vector population is given by *N' _{v}* (f) =

*A*

_{v}- /u_{v}N_{v}.The solution of this equation is

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

Also

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,

*A*0,0). To obtain the basic reproduction number, form the matrices

_{v}/ ц_{у},

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 *R _{2} > Ri,* and

*Ro > 2R*

_{X}.Endemic equilibrium: The EE point of the system is given by *E _{x} =(S’*

_{h},Eh, Г

*where*

_{н}, R’_{h}, S’_{v}, E‘_{v}, I'_{v}j,

where *a = pk _{x}k_{2}kip_{H}p_{v}Pv >* 0,
and

*'Ro = Ro + 2'R*(1 - 7?o).

Note that *Ro > 2R.* We have *a >* 0. Now, *c* > 0 when *R _{0} <* 1, and

*c <*0, when

*R*> 1. The signs of

_{tl}*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

*b*-

^{2}*4ac >*0. That is, two positive equilibrium points exist. For

*'Ro =*1,

*c =*0, and if

*b*< 0, then a unique positive solution

*I‘*exists. Now, consider the case

_{H}= -b/a*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,

*p*0.01,

_{H}=*S„ =*0.3,

*ju*0.003^,, = 1.3,Л

_{v}=_{н}= 0.4,/7 = 0.11,

*p =*0.029,y = 0.0614799, for the two cases: (i)

*p*0.00005. The value of

_{H}=*Ro*is obtained as

*'Ro =*0.505787. (ii)

*p*0.0005. The value of

_{H}=*Ro*is obtained as

*'R*1.59964. Initial conditions are taken as [100,2,3,1,500,5,70]. For this set of parameter values, we obtain

_{0}=

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 E_{0} = *( _{H}/p_{H} ,0,0,0, A_{v}/p_{v},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

*R*and obtained a characteristic equation, which is of order six. We have derived the Jacobian matrix considering

_{H}*R*and obtained a characteristic equation of order seven. The results in both the cases are the same.

_{H}FIGURE 5.4

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