Global Stability of Zika Virus Dynamics

The few mathematical models available in the literature to describe the dynamics of Zika virus are still in their initial stage of stability and bifurcation analysis, and they were in part developed as a response to the most recent outbreaks, including the one in Brazil in 2015, which has also given more hints to its association with Guillain–Barre Syndrome (GBS) and microcephaly. The interaction between and the effects of vector and human transmission are a central part of these models. This work aims at extending and generalizing current research on mathematical models of Zika virus dynamics by providing rigorous global stability analyses of the models. In particular, for disease-free equilibria, appropriate Lyapunov functions are constructed using a compartmental approach and a matrix-theoretic method, whereas for endemic equilibria, a relatively recent graph-theoretic method is used. Numerical evidence of the existence of a transcritical bifurcation is also discussed.


Introduction
Mathematical models on disease epidemics are ubiquitous in the literature (see e.g. [1,[8][9][10][11][12][14][15][16]21,22]), covering a variety of infectious diseases (e.g. cholera, dengue, influenza and malaria), under various considerations and variations, fully describing the disease from a biological and physiological point of view, and providing rigorous and detailed mathematical analyses, with the goal of helping understand the dynamics of the disease, the conditions under which an outbreak can occur, and how the disease would spread. These models are important in fully understanding how the disease evolves with time, and therefore in expanding emergency management of the disease, in creating new and better health-care resources, and ultimately in preventing the initiation and spread of the disease. Mathematical models on the dynamics of Zika virus are relatively new, and they are mainly focused on the estimation of the main epidemiological parameters of the model-most importantly, the basic reproduction number; they also discuss the generation time (the time between sequential infections), the critical population size (under which the disease will most likely die out), or the impact of short-term mobility between two interconnected communities [5,6,13,18]. The basic reproduction number R 0 , as defined in [3,17], is the average number of secondary infections that an infected individual would cause in a completely susceptible population.
We refer the reader to those articles and the references therein for a detailed biological and epidemiological description of the Zika virus. Here we give a brief summary of the main facts about this virus, as presented in their findings. Since it was first isolated from humans in 1954, there have been Zika virus outbreaks in Micronesia (2007), Polynesia (2013), and Brazil (2015), and the virus has been associated to microcephaly and to GBS. See [4,6] and the references therein for further details. This virus is mainly transmitted to humans by mosquito bites of the Aedes genus, a type of mosquito that is known to also transmit dengue and other diseases. Indeed, some of the results in the articles mentioned above suggest that the dynamics of Zika virus might be similar to that of dengue; in particular, large sporadic outbreaks, and a large proportion of asymptomatic cases are typical in both diseases. The disease can also be sexually transmitted but to a much smaller scale, so that by itself it might not be able to initiate or even sustain an outbreak, but it does increase the risk of infection. The type and the intensity of the transmission varies according to the stage of the infection. During the incubation period, infected humans can infect susceptible mosquitoes and humans. After this period, symptomatic humans are more contagious to mosquitoes, and finally at the convalescence stage they can no longer infect mosquitoes, but they can still spread the infection to humans through sexual activity. It is assumed that Zika virus generates lifelong immunity, which might explain a delay of at least a decade before large Zika virus epidemics recur.

The Original Model
Our main goal in this work is to extend previous research on the Zika virus by providing a rigorous global stability analysis of a mathematical model describing the dynamics of this virus. We focus on the deterministic SEIR model presented by Gao et al. [6], where the human population is divided into susceptible S h (t), exposed E h (t), symptomatically infected I h1 (t), convalescent I h2 (t), asymptomatically infected A h (t), and recovered R h (t) populations. The mosquito population is divided into susceptible S v (t), exposed E v (t), and infectious I v (t) populations. The total populations of humans and mosquitoes, that is, The model proposed in [6] is given bẏ (2.1) In this model, susceptible humans, S h , can get infected from sexual contact with E h , I h1 , or I h2 , at a rate β, or from being bitten by infected mosquitoes, I v , with biting rate a and transmission probability b. The parameters 0 < κ < 1 and 0 < τ < 1 denote the transmissibility of E h and I h2 to S h respectively. Once in contact with the virus, a susceptible human has a θ probability of entering the exposed population, otherwise they will enter the asymptomatic population. Asymptomatic humans remain in that stage through the duration 1/γ h of the asymptomatic infection, then enter the recovered population. Exposed humans move to being symptomatically infected, then convalescent, before being considered recovered.
Susceptible mosquitoes can be infected by coming into contact with either E h or I h1 , with transmission probability c. Once in contact with the disease, mosquitoes become exposed, and then they become infectious, (with incubation period 1/ν v ), and they will remain so for the duration of their lifespan 1/μ v . There is some evidence that asymptomatically infected individuals may be capable of transmitting the disease, but it is not well documented, and therefore not considered in this model.
The feasible region for model (2.1): is compact and invariant with respect to system (2.1). As for an equilibrium point, with all parameters being positive, setting all equations in (2.1) equal to zero implies that the state variables must be of the form Thus, no endemic equilibrium exists, and in principle a disease-free equilibrium point (DFE) is not necessarily unique. If we let S * h := S h (0) > 0, then R * h := R h (0) = N h − S * h is uniquely determined, and so is the DFE of the form Here we provide global stability analysis of the DFE E * 0 for completeness, and as an extension of the work in Gao et al. [6] Global Stability Under the assumption above, and using a matrix-theoretic method from [21], here we establish a global stability result of the equilibrium point E * 0 of model (2.1), in terms of the basic reproduction number R 0 . Theorem 1 If R 0 < 1, then the disease-free equilibrium point E * 0 of (2.1) is globally asymptotically stable in 1 .
Proof We start by writing the system in compartmental form [21] by splitting the variables in system (2.1) into a disease compartment x ∈ R 5 , and a disease-free compartment y ∈ R 4 , where We remark that following [17], here we are considering exposed mosquito population E v as part of the "diseased" population. Then, we define with F i representing the rate of new infections, and V i representing transition terms (see [21] for details on these functions). This implies that the system (2.1) can be split as Then, the next generation matrix [23] is given by and the basic reproduction number is given by the spectral radius ρ of F V −1 : We then define (see [21]): The inequality above is true because clearly S h (0) ≥ S h (t), for all t ≥ 0; that is S * h ≥ S h . Now observe that if R 0 < 1, then Q ≤ 0 in 1 , and thus Q is a Lyapunov function for system (2.1). In addition, Q = 0 implies that ω T x = 0 and therefore x = 0 (no disease), but one can readily show that the largest (and only) invariant set in the disease-free system . This equilibrium point is globally asymptotically stable (GAS) in the disease-free system. Applying LaSalle's invariance principle, one concludes that R 0 < 1 implies that the DFE E * 0 is globally asymptotically stable in 1 .

Remark 1 The basic reproduction number in (2.3) can be expressed as
where R hh and R hv represent the basic reproduction numbers of human and vectorial transmission respectively. These quantities R hh and R hv are equivalent to those given in [6], after rescaling by S h * N h , and are directly related to the next generation matrix and the disease-free equilibrium point E * 0 .

A Modified Model
The mathematical model studied in "The original model" clearly reflects the fact that for a basic reproduction number satisfying R 0 < 1, the disease is supposed to die out (all solutions curves converge to the DFE point) so that an outbreak does not happen. However, the simplicity of the model does not allow for a mathematical statement on conditions under which such an outbreak could occur. Studying stability conditions for a disease-free equilibrium point is an important part in the understanding of disease dynamics, but studying the possible existence of an endemic equilibrium (EE) and its stability properties is at least as equally important, especially when developing strategies for prevention and control of a disease.
In this section, we modify the original model (2.1) to allow the possibility of existence of an endemic equilibrium, by including some additional terms in the interaction of the state variables. This new model reflects in part some of the biological context in the model studied in [17], but keeps most of the dynamics of (2.1). Here we study the model Comparing these equations with the original model (2.1), we observe the introduction of birth and death rates for humans and mosquitoes. In particular, μ h is the human birth and death rate, which can be understood as the baseline mortality of humans. For mosquitoes, v is the natural birth rate, r v is the intrinsic growth rate, with r v = v − μ v , and K v is their carrying capacity. As explained in [17], the expression used to model the recruitment rate of susceptible mosquitoes models the fact that they are recruited from pupal stage, and that it should account for egg-laying rate, as well as survival of eggs, larvae and pupae. The term is density dependent because it is assumed that mosquito population depends on availability of egg-laying sites and competition between larvae. Clearly, this is more realistic than just using the term μ v in the corresponding equation of system (2.1). All other parameters are defined as in the original model (2.1). Letting , the feasible region is given by which is compact and invariant with respect to the system (3.1). This system has a DFE point as well as an EE point. The DFE point is

Global Stability of DFE
Following a similar approach as in "Global stability", here we establish global stability of the DFE point E 0 (Fig. 1).
Theorem 2 If R 0 < 1, then the disease-free equilibrium point E 0 of (3.1) is globally asymptotically stable in 2 .
Proof We split the variables into the same compartments as in the proof of Theorem 1. The corresponding functions F and V are now given by: The next generation matrix F V −1 has exactly the same pattern as the one in (2.2), this time with With these values, the basic reproduction number for (3.1) is given by As in the proof of Theorem 1, a left eigenvector ω > 0 of V −1 F associated to R 0 can be computed, and the function Similarly as before, if R 0 < 1, Q is a Lyapunov function for system (3.1). In addition, Q = 0 implies that x = 0, the largest (and only) invariant set in the disease-free system is the disease-free equilibrium point, which is GAS in the disease-free system, and by applying LaSalle's invariance principle, one concludes that R 0 < 1 implies that E 0 point is globally asymptotically stable in 2 .

Remark 2 As before, the basic reproduction number in (3.2) can be expressed as
where R hh and R hv represent the basic reproduction numbers of human and vectorial transmission respectively. One can actually write Note: Observe that the basic reproduction number in (2.4) is a particular case of (3.4) with μ h = 0 and K v = N v .

Global Stability of EE
In "Global stability of DFE", we proved that if R 0 < 1, the disease-free equilibrium is globally asymptotically stable. When this condition is reversed, then such equilibrium loses stability and one can prove the existence of an endemic equilibrium (EE) as well as its global stability properties. We start with a statement on the existence of such equilibrium point, which we denote as Theorem 3 Consider system (3.1), and let F and V be defined as above. If R 0 > 1, then the system has an endemic equilibrium.
Proof First observe that the function f in (3.3) satisfies f (x, y 0 ) = 0. This implies that the Lyapunov function Q satisfies Q = (R 0 − 1)ω T x > 0, indicating that the diseasefree equilibrium E 0 is unstable, and therefore that solutions move away from E 0 , which is equivalent to uniform persistence of (3.1) (see Theorem 4.3 in [7], where E 0 plays the role of the invariant set N ). Then, according to Theorem D.3 in [20], the invariance of 2 and the uniform persistence guarantee the existence of the endemic equilibrium.
To prove a global stability result for the endemic equilibrium (EE) of system (3.1), we use a graph-theoretic approach introduced by Shuai and Van der Driessche [21]. Here we briefly summarize the main ideas; for details, the reader is referred to the aforementioned article.
Given a weighted digraph with p vertices we define as usual the p × p weight matrix A with a i j > 0 if it exists, and a i j = 0 otherwise, and we will denote such weighted digraph as (A). The Laplacian of (A) is defined as Let c i be the cofactor of l ii . We will use the following combinatorial identities: If a i j > 0 and the out-degree of vertex j satisfies d + ( j) = 1, for some i, j, then We state a theorem that will be used in the construction of an appropriate Lyapunov function. For details and a proof of this theorem, see [21].

7)
and assume that (i) There exist functions D i : E → R, G i j : E → R, and constants a i j ≥ 0 such that where S(C) denotes the set of all arcs in C. Then, there exist constants c i ≥ 0, i = 1, . . . , p (as defined above), such that the function .
We now give the main result on global stability of the EE point E * of (3.1). In the theorem below, we use the fact that N h and N v are constants, and represent the total human and vector populations respectively, for any t ≥ 0. Thus, we have Theorem 5 If R 0 > 1, the endemic equilibrium point E * of (3.1) is globally asymptotically stable in 2 .
Proof We start by defining the functions

Differentiation yields
=: a 10,9 G 10,9 + a 10,8 G 10,8 The associated digraph is shown in Fig. 2. One can readily verify that along each cycle C of such graph one gets (i, j)∈s(C) G i j (z) = 0, where s(C) denotes the arc set of C. Then, Theorem 4 guarantees the existence of constants c i in the definition of a Lyapunov function D, and it also guarantees that such Lyapunov function satisfies D ≤ 0. Such constants are found using the combinatorial identities in (3.5) and (3.6), which give Thus, the Lyapunov function, D, and its derivative (for any i = 1, 2, 3, 4 and any j = 8, 9) are given by From here we see that the only invariant set for (3.1) whereḊ = 0 is the equilibrium point E * . Thus, by LaSalle's Invariance Principle, we conclude that the EE is globally asymptotically stable in 2 , and thus unique.

Bifurcations
In the previous sections we proved that when the basic reproduction number satisfies R 0 < 1, the DFE E 0 is a stable equilibrium of system (3.1), and when R 0 > 1, the endemic equilibrium E * is stable instead. That is, at R 0 = 1 there is a exchange of stability properties between these two equilibria. Even more, for R 0 < 1 there is an unstable "endemic" equilibrium. This is a clear indication of the presence of a transcritical bifurcation when R 0 = 1. Using Matcont [2], we computed a transcritical bifurcation when a 0 = 0.3287922 (the biting rate), which corresponds to R 0 = 1. When a < a 0 , the DFE is stable while the EE is unstable.

Conclusions and Final Remarks
A mathematical model proposed by Gao et al. [6] to study the dynamics of Zika virus has been discussed and then extended in this work. We have explicitly computed the diseasefree equilibrium point of the original model and studied its global stability properties. The basic reproduction number was computed as the spectral radius of the corresponding next generation matrix and is comparable to the one found in [6]. As this original model does not posses an endemic equilibrium, we have modified it so that such equilibrium does exist, but keeping most of the dynamics of the original equations. For this new model, besides proving the existence of an endemic equilibrium point, we have also established theorems that provide conditions for the global stability of both the disease-free and the endemic equilibrium points. A matrix-theoretic and a graph-theoretic method, as presented by Shuai and Van den Driessche [21], were used to construct the appropriate Lyapunov functions. The basic reproduction number obtained for the modified model is in a sense a generalization of the one given in [6], and its variation from R 0 = 1 indicates an exchange of stability between the two equilibria. More generally, the original model (2.1) is a particular case of the modified model (3.1) by letting μ h = 0 and K v = N v . Besides the numerical evidence provided by phase portraits that a transcritical bifurcation exists, we have computed such bifurcation using Matcont, with the mosquito biting rate a as the bifurcation parameter.
Some variations and generalizations of these models would be worth investigating, including the study of multi-patch models as proposed in [18], the human connectivity of communities, as in [1,5,9,19], an individual-level heterogeneity in the probability of infection, and an estimation of the so-called "critical community size", below which the disease will die out [13]. Specific heterogeneities in humans such as gender and socioeconomics have also been proposed as future study by Gao et al. [6]. The introduction of these features into the current models should bring richer dynamics and possibly provide conditions for the existence of other bifurcations, as they are typically found in other models of disease epidemics.