key: cord-0047051-6n7r5gi4
authors: Danbaba, UA; Garba, SM
title: Stability Analysis and Optimal Control for Yellow Fever Model with Vertical Transmission
date: 2020-07-06
journal: Int J Appl Comput Math
DOI: 10.1007/s40819-020-00860-z
sha: b8efcb07852a1c82335c306755c797b89c88bc7f
doc_id: 47051
cord_uid: 6n7r5gi4

In this study, a deterministic model for the transmission dynamics of yellow fever (YF) in a human–mosquito setting in the presence of control measures is constructed and rigorously analyzed. In addition to horizontal transmissions, vertical transmission within mosquito population is incorporated. Analysis of the mosquito-only component of the model shows that the reduced model has a mosquito-extinction equilibrium, which is globally-asymptotically stable whenever the basic offspring number [Formula: see text] is less than unity. The vaccinated and type reproduction numbers of the full-model are computed. Condition for global-asymptotic stability of the disease-free equilibrium of the model when [Formula: see text] is presented. It is shown that, fractional dosing of YF vaccine does not meet YF vaccination requirements. Optimal control theory is applied to the model to characterize the controls parameters. Using Pontryagin’s maximum principle and modified forward–backward sweep technique, the necessary conditions for existence of solutions to the optimal control problem is determined. Numerical simulations of the models to assess the effect of fractional vaccine dosing on the disease dynamics and global sensitivity analysis are presented.

Yellow fever (YF) is an acute viral haemorrhagic fever that is transmitted by mosquitoes of the Aedes and Haemogogus species. It is endemic in Africa, Central and South America, where approximately one billion people in fourty seven countries are at risk. Symptoms of the disease include fever, headache, jaundice, muscle pain, nausea, vomiting and fatigue [5, 44, 47, 48] . In broad terms, YF can either be jungle or urban. Jungle YF occurs in tropical rain-forest. In Africa, it is usually transmitted by Aedes africanus while in South America by Haemagogus species [6, 47, 48] . Urban YF is transmitted by Aedes aegypti, that it characterized by rapid amplification, capacity for international spread and has devastating effect on public health, delivery in general [47, 48] . In order to ensure adequate supply of vaccine especially in high risk regions, the Eliminate Yellow Fever Epidemics (EYE) strategy, steered by World Health Organization, UNICEF and Gavi, the Vaccine Alliance was inaugurated. The vision of EYE is to have a world without YF epidemics. Its mission is to coordinate international action and help countries at risk of the disease to prevent outbreaks and prepare for the inevitable cases, to minimize suffering, damage and spread through early and reliable detection as well as a rapid and appropriate response. The initiative has three strategic objectives: they include protecting at-risk populations, preventing international spread, and containing outbreaks rapidly [47] .

Yellow fever has attracted less modeling studies when compared with other mosquito borne diseases such as malaria, dengue, West Nile and Zika virus. To study population dynamics of Aedes aegypti (mosquitoes responsible for YF transmission), Dye [17] proposed an appropriate continuous time model that described a field population of adult Aedes aegypti mosquito. Recently, Martorano et al., constructed and analysed a compartmental model for the transmission of YF with vaccination [32] , although vertical transmission is not accounted for in their model, both aquatic and non-aquatic stages of mosquito development were considered. An urban YF epidemic model was also formulated and used to study the 2016 YF outbreak in Luanda, Angola by Zhao et al. [50] . Their model successfully fits the time series of weekly reported YF cases and deaths during the epidemic in Angola [50] . Monica et al. [34] looked at a YF model in a human-vector-primate setting with vaccination in human population. To estimate the incubation periods of YF virus in both human and mosquito populations, four statistical models of incubation periods were fitted with historical data in [26] .

In this work, we extend the model in [32] by incorporating vertical transmission in mosquito population. In addition to vaccination of susceptible humans, the proposed model also incorporates the use of treated bed nets, larvicides and adulticides in mosquito control. The work is organized as follows: Introduction and short review of relevant literature is presented in "Introduction" section. A deterministic model for the transmission dynamics of YF is constructed and analyzed for its basic dynamical features in "Model Formulation" section. Analyses of the full model is provided in "Analysis of the Full Model" section. Threshold quantities and stability analysis of equilibria are also explored in this section. Optimal control form of the model is presented and analyzed in "YF Model for Optimal Control" section. Sensitivity analysis and numerical simulations are presented in "Sensitivity Analysis and Numerical Simulation" section.

Following compartmental modeling approach, the total human population at time t, denoted by, N H (t), is divided into five mutually exclusive compartments of susceptible (S H (t)), vaccinated (V H (t)), exposed (E H (t)), infected (I H (t)) and recovered (R H (t)) humans. So that

Mosquito population is split into aquatic (immature) and non-aquatic (mature) stages. For mathematical tractability, different development stages of the aquatic mosquito population (eggs, larvae and pupae) are lumped into a single compartment A(t). Furthermore, to incorporate vertical transmission, the aquatic mosquito population are further divided into infectious (A I (t)) and non-infectious (A N (t)) mosquitoes. So that the total mosquito population at the aquatic stage at time t, is given by

Similarly, the total mosquito population at non-aquatic stage (adult) at time t, denoted by N V (t), is sub-divided into susceptible (S V (t)) and infectious (I V (t)) mosquitoes. So that N V (t) = S V (t) + I V (t).

The frequency-dependent (standard) incidence function is the most widely used form of incidence in vector borne disease models. Infection from mosquitoes to humans occur after an infectious mosquito bites a susceptible human at a rate b V H , let ρ V H be a transmission probability from an infectious mosquito to susceptible human, then the infection rate of humans is β V H = ρ V H b V H . Therefore the force of infection in humans is given by

Similarly, let β H V = ρ H V b H V be the rate at which susceptible mosquitoes acquire infection from infectious human, where ρ H V is the probability of transmission from an infectious human to a susceptible mosquito and b H V is the biting rate of a susceptible mosquito. Then the force of infection in mosquito population (due to horizontal transmission) is given by

Since mosquitoes bite both susceptible and infected humans, for the total number of bites to be conserved, it is assumed that the total number of bites by the mosquitoes must be equal to the total number of bites received by humans (and this depends on the total sizes of the populations of humans and mosquitoes), see [8, 10, 12, 21, 22, 36] . Thus

Substituting Eq. (3) into Eq. (1), we have The population of aquatic mosquitoes (eggs, larvae and pupae) increases through oviposition by reproductive mosquitoes (infected and non-infected) at a rate φ V . The aquatic mosquitoes mature to adulthood at a rate b V , die naturally at a rate μ A and due to the use of larvicides at a rate c L = r L e L (where r L is the rate of applying larvicides and e L is the efficacy of larvicides).

The population of susceptible adult mosquitoes is generated by maturation of noninfectious mosquitoes from the aquatic stage at the rate b V and decreases by infection and move to infectious class at a rate (1 − r B B )λ V , they die naturally at a rate μ V and due to the use of adulticides at a rate c A = r A e A (where r A is the rate of applying adulticides and e A is the efficacy of adulticides). Finally, The population of infectious adult mosquitoes is generated by maturation of infected mosquitoes from the aquatic stage at the rate b V and by infection of susceptible mosquitoes. This population is reduced by natural death at the rate μ V and due to the use of adulticides at a rate c A .

The time independent YF transmission model with vaccination (and vertical transmission in mosquitoes) is represented by the following system of non-linear ordinary differential equations (a flow diagram of the model is depicted in Fig. 1 and the state variables and parameters of the model are described in Table 1) : 

It is assumed that, all the model parameters are positive and initial conditions are non-negative. In addition, let A N + A I = A so that

which is the standard formulation for aquatic mosquitoes with density dependent and independent death rate, see for instance [15, 16] .

The following biologically feasible region of the model (5)

is positively-invariant and attracting.

Proof It is easy to see that solution to the system (5) exists locally and it is unique (system (5) is C 1 in R 10 + ). Observe that A N + A I ≤ K, thus A N (t) ≤ K and A I (t) ≤ K. Therefore by Gronwall's lemma we have

which are bounded and hence solution exists for all t ≥ 0. In addition,

Consequently, solution of the system (5) with initial condition in remains in for all t > 0 (the ω-limits set of the system are contained in ).

Having obtained the positively-invariant and attracting domain for the system (5), it is sufficient to consider the asymptotic properties of dynamics of the flow generated by the system.

Consider the mosquito component of the model given by (5) in the absence of interaction with humans. By direct computation, we obtained a threshold termed as the basic offspring number (N 0 ) given by

. (9) It is defined as the average number of offspring produced by a female mosquito in her entire lifespan in the absence of interaction with humans. It can be interpreted as follows. The average time spent by mosquito in the aquatic stage is given by

where b V is the rate at which aquatic mosquitoes mature into an adult mosquito, so that the probability that an aquatic mosquito develops into an adult female mosquito is given by

The average life expectancy of an adult female mosquito is given by 1 μ V +c A , so that the average eggs laid by an adult female mosquito throughout her life span is given by

where φ V is the oviposition rate of a female mosquito. Thus, the product of (10) and (11) gives (9), the basic offspring number of the mosquito-only population model. The mosquito component of model (5) has an extinction disease-free equilibrium obtained when N 0 ≤ 1, denoted by E 0 , given by

and non-extinction disease-free equilibrium obtained when N 0 > 1, denoted by E 1 , that is given by

The mosquito extinction equilibrium, E 0 , is globally-asymptotically stable (GAS) when N 0 ≤ 1 and unstable otherwise. The equilibrium E 1 exists and it is locallyasymptotically stable (LAS) when N 0 > 1.

The proof of the theorem is given in "Appendix A". The epidemiological implication of Theorem (2.2) is that, if the basic offspring number can be brought to a value below unity, then the mosquito population goes to extinction and horizontal transmission can be avoided. It is worth mentioning that this result is not attainable.

The disease-free equilibrium of the model given by (5) depends on N 0 . If N 0 ≤ 1 a mosquito extinction DFE (E 2 ) is obtained, while a mosquito persistent equilibrium (E 3 ) is obtained when N 0 > 1. Thus

For the case, when N 0 < 1, the associated reproduction number obtained by linearizing the system about E 2 is given by

The detail computation is provided in "Appendix B".

The mosquito extinction DFE given by, E 2 , is locally-asymptotically stable if the vectorial vertical transmission reproduction number R vv = η V ≤ 1 and unstable otherwise [46] .

The DFE, E 2 , can be shown to be globally-asymptotically stable under the same condition. It is worth mentioning that, the mosquito extinction DFE, E 2 , is less tractable (due to the absence of mosquito in the population). For the DFE, E 3 , obtained in the presence of mosquitoes (N 0 > 1), the vaccinated reproduction number (R 0v ) is given by

where c B = r B B is the rate of reducing contact between humans and mosquitoes through the use of bed nets. The threshold quantity, R 0v is the average number of new secondary cases that one infected individuals can produce in a totally naive population, where a fraction of the population is vaccinated. The computation of the threshold is presented in "Appendix C". The disease can be controlled in a community through appropriate measures such as reducing the carrying capacity of mosquito population, use of larvicides, adulticides or repellents to the extent of lowering R 0v to a value below unity. Effective mosquito control strives to prevent large swarms of adult mosquitoes in an environment through the application of chemical substances called adulticides. They can be applied either aerially or on the ground. Droplets of the chemicals that make physical contact with mosquitoes usually kill them, while large droplets that missed target and settle on surfaces may cause undesirable harm. To achieve small droplets, adulticides are mostly applied in the air as a very fine ultra low-volume (ULV) droplet spray from a truck or aircraft, it is usually organophosphate insecticides and/or synthetic pyrethroids and their combinations [40, 42] . Through aerial or ground applications of larvicides, large population of aquatic mosquitoes can be killed. This method is often more effective and environmentally friendly than the use of adulticides. It should be noted that, larvicides are only administered at identified suspected breeding sites, where they are expected to clear populations of aquatic mosquitoes, conventional larvicides kill aquatic mosquitoes at every stage and therefore, they can be applied whenever necessarily. Some larvicides agents are specific to mosquitoes and when used according to directions will have relatively little impact on the environment and human health. They can prevent the emergence of adult mosquitoes for up to 1 month, which decreases labour costs [47] . In essence, larvicides can be specific for mosquitoes, they have minimal impact on other organisms and often penetrate even dense objects.

Here we analyse the potential impact of a single dose and a fractional dosing of vaccine.

Since not all vaccine have positive impact in a population, it is therefore instructive to first of all assess the impact of vaccine.

In the absence of vaccination (S * H = N * H when V * H = 0), the vaccinated reproduction reduces to

Thus, vaccination of individuals will have positive impact in the community by reducing the value of the associated reproduction number R 0 . Furthermore, the impact of vaccination can be analysed qualitatively by differentiating R 0v with respect to the fraction of vaccinated individuals (V ). It can be shown that

Thus, R 0v is a decreasing function of V . Since the reproduction number measures disease burden, the vaccination will have a positive impact in disease control.

Based on the available clinical data [47] , the minimum standard dose administered should preferentially contains 3000 international units (IU)/dose, but no less than 1000 IU/dose. Let

be the fraction of the vaccinated individuals at steady-state (when standard dose of YF vaccine is issued). Then solving for R 0v = 1 we obtained

Thus, for the vaccination to be effective in bringing R 0v < 1, the fraction of vaccinated individuals (V c ) at steady-state must be greater than the vaccinated threshold ratio (V > V c ) defined in Eq. (15) . Notice that if η V = 1, that is the case when all eggs laid by infected female mosquitoes are infected, then it is highly unlikely for any vaccine coverage to bring R 0v to a value less than unity (in this settings) since the critical vaccination rate reduces to V c = 1 .

The proof follows from Theorem 2 of [46] and the fact that R 0v < 1 if and only if η V < 1 and V > V c .

Consider the model (5) with N 0 > 1. We claim the following result (the proof of the Theorem is presented in "Appendix D").

The yellow fever model (5) undergoes backward bifurcation at R 0v = 1 whenever the bifurcation coefficient a, given by (D.3) is positive.

The epidemiological implication of the phenomenon of backward bifurcation is that the classical requirement of R 0V < 1 is, although necessary, no longer sufficient for the effective control of the disease in the population [3] . Hence, the presence of backward bifurcation makes the feasibility of the effective control of YF in a population difficult. In the next section, the possible cause(s) of this phenomenon is(are) explored.

Here, a global asymptotic stability of the non-extinction equilibrium, (E 3 ), is presented To confirm the absence of backward bifurcation in the model (5) for a special case when disease induced mortality in human is negligible (δ H = 0). Consider the model (5) with N 0 > 1 and δ H = 0. (7) given by

We claim the following result Theorem 3.5 The non-extinction equilibrium (E 3 ) of the model (5) is globally-asymptotically stable in the positively invariant set * if δ H = 0 and R 0v ≤ 1.

The proof of the Theorem is presented in "Appendix E". This result shows that, in the absence of disease induced death, the DFE of the model (5) is GAS. Hence, the classical requirement of R 0v ≤ 1 is necessary and sufficient condition for disease elimination from the community provided δ H = 0 (thus, YF will be effectively controlled or eliminated from the population if R 0v ≤ 1). This result is consistent with that in [21] , which suggests that disease induced death in humans is the main cause for the emergence of the backward bifurcation phenomenon in this setting.

The best way to stretch vaccine supplies and protect as many people as possible to stop the spread of yellow fever in emergency situations is by using fractional dosing. Based on the available evidence, the Strategic Advisory Group of Experts (SAGE) on Immunization affirms that a fractional dose can be used as part of an exceptional response when there is a large outbreak and a shortage of vaccine [47, 48] . In the case of dose fractionation, a smaller amount of antigen would be used per dose in order to increase the number of persons who can be vaccinated with a given quantity of vaccine. Studies show that the yellow fever vaccine given as one fifth of the regular dose, still provides full immunity against the disease for at least 12 months and likely longer [47, 48] . This strategy was previously proposed Figure 2 shows the simulation of the vaccinated reproduction number as a function of vaccine efficacy with single dose and fractionated threefold. Although when = 0, all the simulations have the same value of R 0v = R 0 (about 1.143), the vaccinated reproduction number becomes less than unity when the vaccine efficacy, > 0.4 (for a single dose), while for threefold fractionated vaccine, a higher vaccine efficacy, is required to possibly bring R 0v to a value below unity. This result is consistent with those in [47] , which stated that a fractional YF vaccination does not meet YF vaccination requirements under the International Health Regulations (IHR) (Fig. 3) .

For homogeneous population, controlling the basic reproduction number will be sufficient for disease control. In the case of heterogeneous populations (usually vector borne diseases) with more than one host types, control is often targeted at one host. The type-reproduction number (T) is a threshold quantity that correctly determines the critical control effort for a heterogeneous populations [25] . A method that is used to estimate the required effort(s) for controlling disease by targeting a specific sub-population of hosts, under the premise that infection may pass through other sub-populations before causing secondary infections is described in [25, 41] . If K is the next generation matrix with large domain and hosts 1, 2 and 3 represent the populations of E H , A I and I V . The type i reproduction number is given by where I is an identity matrix, P is a projection matrix and e is a unit vector with all elements equal to zero except the ith. Let

Notice that k i j is the expected number of cases of type i produced by one infected individual of type j, so that from (16) the type-reproduction number for exposed humans is

.

Observe that R 0v < 1 implies

Thus T 1 < 1 implies R 0v < 1. Similarly, it can be shown that the infected aquatic mosquito type-reproduction number (T 2 ) satisfies T 2 < 1 and the infectious adult mosquito typereproduction number (T 3 ) also satisfies T 3 < 1 whenever R 0v < 1. 

Usually, incidences of YF and other vector borne diseases are seasonality dependent with their peaks during warm and rainy seasons, therefore it is reasonable to integrate time dependent controls in the model, the goal of which is to show the possibility of implementing time dependent controls while minimizing implementation cost. Let the time dependent effort in preventing human-mosquito contacts through the use of treated bed nets be u 1 (t), so that the contact rate between mosquitoes and humans reduces by a factor (1 − u 1 (t) ) where 0 ≤ u 1 (t) ≤ 1. The effort in vaccinating humans is u 2 (t) : 0 ≤ u 2 (t) ≤ 1. Similarly, the effort in the application of larvicides is u 3 (t) : 0 ≤ u 3 (t) ≤ 1, while that of spraying adulticides is u 4 (t) : 0 ≤ u 4 (t) ≤ 1. For instance, there is no any effort in controlling mosquitoes when u 3 (t) = u 4 (t) = 0, while aquatic and mature mosquitoes die at maximum possible rates c L and c A respectively, when u 3 (t) = u 4 (t) = 1. Maximum control is attained by the use of bed nets and vaccination when u 1 (t) = 1 and u 2 (t) = 1, respectively and no effort invested when u 1 (t) = u 2 (t) = 0. The autonomous system given by (5) is extended to include the aforementioned time dependent controls. Let (t) = ω H (1 − u 2 (t)), then the non-autonomous version of the model (5) is given by

Mosquitoes

Following the non-autonomous system given by (18) , an optimal control problem is formulated with the following objective (cost) function.

The interval [0, T ] represents the time through which various control measures are implemented. The cost incurred due to human infection of YF (which is proportional to the number of infected individuals) over the period of intervention is given by

where B 1 and B 2 are positive weight constants associated with exposed and infected humans, respectively. Similarly, the cost due to the presence of mosquitoes in the community, which is proportional to the number of aquatic and adult mosquitoes is given by

where B 3 and B 4 are positive weight constants. Because of the short supply of YF vaccine, in order to optimize the available vaccines, there is need to minimize the total number of vaccines used over the period of intervention, thus the integral

which measures the total number of vaccinated individuals during the period of intervention included in the objective functional, with B 5 being a positive weight constant. The positive terms D 1 , D 2 , D 3 and D 4 are weight constants for efforts in the use of bed nets, vaccination, larvicides and adulticides, respectively, and regularize the optimal control. D 1 u 2 1 , D 2 u 2 2 , D 3 u 2 3 , and D 4 u 2 4 describe the cost associated with the aforementioned prevention and control measures. The degree of the cost functions follow from the non-linearity of controls and the convexity of quadratic functions [43] . Similar assumption has been used in optimal control problems in epidemiology, see for instance [7, 29, 35, 39, 43] and some of the references therein.

The aim is to minimize the number of infected humans and total mosquito population while optimizing limited vaccines and keeping the cost of vaccination, use of treated nets and application of pesticides low. Therefore we seek to optimize u * 1 , u * 2 , u * 3 and u * 4 such that

where

is the control set. The impact of each control does depends on adherence and effort, if for example u 1 = 1, production and distribution of bed nets is at maximum, but its impact also depends on c B , likewise the remaining control functions.

The existence of optimal control solution can be established using Theorem 4.1 and Corollary 4.1 of [19] . , u 2 , u 3 , u 4 ) over G.

Proof Clearly the set of controls and state variables are non-empty and the control set G is closed and convex. The integrand of the objective functional is convex on G. Furthermore, the model is linear in the control variables and bounded by a linear system in the state variables, thus, the existence of an optimal control is guaranteed [7, 19] .

The necessarily conditions that optimal controls and their corresponding states must satisfy are derived using Pontryagin's Maximum Principle [35, 38] , where the problem of finding time-dependent control variables u * 1 (t), u * 2 (t), u * 3 (t) and u * 4 (t) that minimize J is equivalent to the problem of minimizing the Hamiltonian function defined as (t, x, u) where g (t, x, u) is the integrand of the objective functional (19) and λ(t) is the adjoint vector such that

The optimality equation is given by

and transversality conditions as

where λ S H , ...λ I V are adjoint functions.

Given an optimal control (u * 1 , u * 2 , u * 3 , u * 4 ) and the corresponding state solutions of the non-autonomous system given by (18) , there exist adjoint functions satisfying

with final time condition as λ i (T ) = 0, i = 1, . . . , 9. In addition, the optimal control u * j , j = 1, 2, 3, 4 are given by

also u * 3 is given by

and u * 4 is given by

Proof The non-autonomous system given by (18) together with the objective functional given by (19) and (20) are converted into a problem of minimizing the Hamiltonian, H , defined by (21) . Therefore applying Pontryagin's Maximum Principle [35, 38] , the proof follows.

In this section, global sensitivity analysis using partial rank correlation coefficient (PRCC) for the basic offspring number and vaccinated reproduction number are conducted. Numerical simulations for the optimal control model given by (18) is also presented.

Local sensitivity analysis is used to provide direct information on the effect of small parameter perturbation, it evaluates the relative change in a function due to change in a single parameter, where other parameters are kept at constant values. It does not indicate the effect of simultaneous large perturbations in all model parameters. Thus, the need for a more robust form of sensitivity analysis for a multidimensional parameter space. The PRCC is a robust sensitivity measure for a non-linear but monotonic relationships between inputs and output, with little correlation between the inputs [31] . The PRCC of the basic offspring number (N 0 ) and that of the vaccinated reproduction number (R 0v ) are computed with parameter ranges as presented in Table 1 Fig. 4 Partial rank correlation coefficient plots of the various parameters of the model (5) using N 0 as the output function of the parameters of the YF model, using N 0 as the response function are depicted in Fig. 4 . The Fig. show that the top parameters that most influences the values of the threshold quantity, N 0 , are mosquito death rates due to adulticides and larvicides (c A and c L , respectively) and the vertical transmission rate (η V ). It can be be seen from the figure that c A and c L are negatively correlated and φ V is positively correlated to N 0 . On the other hand, R 0v is most positively correlated to η V , which is followed by β H V then K, and it is most negatively correlated to c A , c B and then b H as presented in Fig. 5 .

Using the forward-backward sweep method as described in [30] , solution of the optimal control problem can be obtained numerically. An initial guess for the optimal control is used to solve the state system in forward time, after which the guessed optimal control and the obtained solution to the state system are used as input to the adjoint system, which is solved numerically in backward scheme using the transversality condition. The controls are then updated using convex combination of the previous controls and the value from the characterizations.

The following numerical values for the model parameters are used as in Table 1 show that the populations of exposed and infected humans reaches the DFE at a faster rates with control than without control. Figure 9a shows the population of recovered humans where the population with control is larger than those without control, this may be attributed to the fact that recovery confers permanent immunity. Simulation of the model for the population of infectious mosquitoes (both aquatic and non-aquatic) is presented in Figure 9b Similarly, control profiles of the model are depicted in Figs. 10 and 11. The simulation of the model for population of exposed humans shows similar dynamics as that of Case 1 as presented in Fig. 12a . Simulation of model for infected humans, where population with control shoot up before falling to the DFE is presented in Fig. 12b . The population of recovered humans shows wider margin between cases with and without control in this instance compared to Case 1 as depicted in Fig. 13a , also, population of infected mosquitoes sporadically bumps up before reaching DFE as shown in Fig. 13b . Notice that, less effort and less costs are expended (in this case) in comparison to Case 1.

A deterministic model for the transmission dynamics of yellow fever in a population is constructed and rigorously analysed. The model with standard incidence formulation incorporates the use of treated bed nets and vaccination as forms of prevention in humans, while the use of larvicides and adulticides are used in controlling mosquito population. Some of the key results obtained include:

• The mosquito-only model has a threshold quantity called the basic offspring number (N 0 ) which described the extinction or persistence of mosquito population. • The model has two disease-free equilibriums, the mosquito-extinction equilibrium (E 2 ) which is globally-asymptotically stable (GAS) when the basic offspring number (N 0 ) is less than unity and the non-mosquito-extinction equilibrium, (E 3 ), which is locallyasymptotically stable when R 0v ≤ 1. • The YF model undergoes the phenomenon of backward bifurcation (a dynamic phenomenon characterized by the co-existence of two stable attractors when the associated reproduction number of the model is less than unity). A condition for the emergence of this phenomenon have been identified. • The DFE, E 3 , of the model is shown to be GAS in the positively-invariant region ,

• For standard dosing of YF vaccine, a high vaccine efficacy (at-least 40%) would be required to reduce R 0v to value below unity. • Fractional dosing of YF vaccine does not meet YF vaccination requirements of yellow fever elimination. This result is consistent with that of International Health Regulations (IHR). However, the vaccine will always have a positive impact in a community. • Using Pontryagin's maximum principle and modified forward-backward sweep technique, the necessary conditions for existence of solutions to the optimal control problem is determined. • Numerical simulations of the non-autonomous model (with optimal control) using different weight constants show that slight increase in control effort show wide difference in the impact of control. • Using sensitivity analysis, it is shown that the vaccinated reproduction is most sensitive and positively correlated to the rate of vertical transmission, thus making it the most important parameter to target in controlling the disease.

Therefore f (b) ≤ 0 ≤ f (0) provided N 0 ≤ 1, thus by Theorem (5.1), the mosquito component of the system given by (5) defines a positive dynamical system on [0, b], moreover, the equilibrium (E 0 ) is GAS on [0, b]. Because q is arbitrary, b can be chosen such that its bigger than any x ∈ R 2 + . Hence the result holds on R 2 + . The second part of the proof follows by linearization.

For E 2 , the matrix of new infection terms and that of transition terms are respectively given by

The next generation matrix with large domain (K L = F V −1 ) is

Thus using the approach of [14] with an auxiliary matrix E, the NGM (K ) is

Therefore the mosquito extinction basic reproduction number, which is the dominant eigen-

(E.4)

Following the above simplification, the system given by (5) can therefore be re-written in a pseudo-triangular form as Then, the associated DFE is GAS in [15, 27] .

Recall that the model given by (5) is defined on a positively invariant domain given by in (7) . Also straightforward computation shows that the eigenvalues of A 11 (x) are real and negative. Therefore conditions 1 and 2 of (5.4) are satisfied, for condition 3 of (5.4), the following definition is used.

where A 1 and A 3 are square matrices of order at least 1 or if A can be transformed into the form (D.7) by simultaneous permutations of rows and columns [18] . It is irreducible otherwise. Alternatively, A square matrix is irreducible if and only if its associated digraph is strongly connected. Fig. 3 is the associated digraph of the matrix A 22 (x), and it is clear that it is strongly connected. Thus condition 3 is satisfied. Furthermore, since

is an upper bound of A 22 (x). For condition 5 of Theorem (5.4), the following result of [27] is applied. In the case ofĀ 22 (x) defined above, we have

.

(D.10)

and it is stable if

(D.12)

It should be noted that, condition (D.12) is a generalization of condition (D.11), which is also equivalent to

Thus, satisfying condition (D.12) is sufficient for the GAS of the DFE.

Dynamic epidemiological models for dengue transmission: a systematic review of structural approaches

Mathematical modeling of sterile insect technology for control of anopheles mosquito

Backward bifurcation analysis of epidemiological model with partial immunity

Qualitative assessment of the role of temperature variations on malaria transmission dynamics

Yellow fever: epidemiology and prevention

Epidemiology and ecology of yellow fever virus

Backward bifurcation and optimal control in transmission dynamics of West Nile virus

A mathematical model for assessing control strategies against West Nile virus

Optimal bed net use for a dengue disease model with mosquito seasonal pattern

Some models for epidemics of vector-transmitted diseases

Dynamical models of tuberculosis and their applications

Modeling the transmission dynamics of Zika with sterile insect technique

Vertical transmission of the yellow fever virus by Aedes aegypti (Diptera, Culicidae): dynamics of infection in F1 adult progeny of orally infected females

The construction of next-generation matrices for compartmental epidemic models

On a temporal model for the Chikungunya disease: modeling, theory and numerics

Vector control for the Chikungunya disease

Models for the population dynamics of the yellow fever mosquito, Aedes aegypti

Special matrices and their applications in numerical mathematics

Deterministic and stochastic optimal control

First evidence of natural vertical transmission of yellow fever virus in Aedes aegypti, its epidemic vector

Backward bifurcations in dengue transmission dynamics

Mathematical analysis of West Nile virus model with discrete delays

Mathematical analysis of an age-structured vaccination model for measles

Efficacy and duration of immunity after yellow fever vaccination: systematic review on the need for a booster every 10 years

The type-reproduction number T in models for infectious disease control

Incubation periods of yellow fever virus

Global asymptotic stability for the disease free equilibrium for epidemiological models

A model of the dynamic of transmission of Malaria, integrating SEIRS, SEIS, SIRS and SIS organization in the host-population

Vaccination and treatment as control interventions in an infectious disease model with their cost optimization

Optimal control applied to biological models

A methodology for performing global uncertainty and sensitivity analysis in systems biology

Equilibrium analysis of a yellow fever dynamical model with vaccination

Comparative safety and immunogenicity of two yellow fever 17D vaccines (ARILVAX and YF-VAX) in a phase III multicenter, double-blind clinical trial

Modelling and stability analysis of SVEIRS yellow fever two host model

An introduction to optimal control with an application in disease modeling

A mathematical model for endemic malaria with variable human and mosquito populations

Pan American Health Organization/World Health Organization: Epidemiological update: yellow fever

The mathematical theory of optimal processes

Optimal vaccination and bednet maintenance for the control of malaria in a region with naturally acquired immunity

Insecticide substitutes for DDT to control mosquitoes may be causes of several diseases

A new method for estimating the effort required to control an infectious disease

Pesticide use for West Nile virus

Optimal control in epidemiology

Global yellow fever vaccination coverage from 1970 to 2016: an adjusted retrospective analysis

Biological and phylogenetic characteristics of yellow fever virus lineages from West Africa

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

A global strategy to eliminate yellow fever epidemics

Yellow fever fact sheet

Fractional dosing of yellow fever vaccine to extend supply: a modeling study

Modelling the large-scale yellow fever outbreak in Luanda, Angola, and the impact of vaccination

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations

Proof Let f : * −→ R 4 be continuous, where * ⊆ R n . Consider a system given bẏ x = f (x), we use the following theorem in [2] in the proof of Theorem 2.2 Theorem 5.1 [2] Let a, b ∈ * be such that a < b, [a, b] ⊆ * and f (b) ≤ 0 ≤ f (a). Thenẋ = f (x) defines a (positive) dynamical system on [a, b] . Moreover, if [a, b] contains a unique equilibrium q then q is globally asymptotically stable on [a, b] .By rewriting the mosquito component of (5) in the form ofẋ = f (x) and considering the interval

For the case of E 3 , applying similar method to that of Appendix B, the NGM with large domain K L is given by(B.1) Thus using the approach of [14] with an auxiliary matrix E, the NGM (K ) isThus, the vaccinated reproduction number which is the dominant eigenvalue of K is given by

Proof We apply a method which is based on the Centre Manifold Theory [11, 46] to prove the existence of backward bifurcation for the model (5) . Let,The transformed model (5) is represented by,The forces of infections in human and mosquito populations are respectively given by,, and, β H V x 4 x 8Thus, the Jacobian matrix (J * ) at the DFE withunder which we obtained the left eigenvector (v) and the right eigenvector (w) corresponding to the zero eigenvalue given byClearly v I ≥ 0, w 1 < 0, w 2 < 0 while w 6 and w 8 can be positive or negative, such choice is justified by Remark 1 of [11] which states;

The requirement that w is non-negative in the theorem is not necessary. When some components in w are negative, we still can apply this theorem, but one has to compare w with the actual equilibrium because the general parametrization of the Centre Manifold before the coordinate change is,provided that x 0 is a non-negative equilibrium of interest (usually x 0 is the disease-free equilibrium). Hence, x 0 − 2bφ V a > 0 requires that w j > 0 whenever x 0 ( j) = 0. If x 0 ( j) > 0, then w( j) need not be positive [11] .It can be verified that vw = 1, thus all the necessarily conditions for the use of the Center Manifold theory are satisfied, and.

The yellow fever model (5) undergoes backward bifurcation at R 0v = 1 whenever the bifurcation coefficient a, given by (D.3) is positiveSince the bifurcation coefficient b is positive, the direction of the bifurcation depends on the sign of a, which can be positive or negative, and a > 0 means the model (5) may undergoes backward bifurcation at R 0V = 1 [11] .

Proof Conditions for global asymptotic stability of the DFE (E 3 ) can be found using a method described in [27] . Similar approach was employed in [15, 16, 28] . Using the property of the DFE, system (5) can be rewritten in a pseudo-triangular form as follows,(E.1)The equation of vaccinated humans is rewritten asSimilarly, the equation of non-infectious aquatic mosquitoes can be rewritten as