Influence of seasonality on Zika virus transmission

1.   Introduction

Zika virus infection is a recurring mosquito-borne flavivirus that it is transmitted through mosquito bites ^[1,2,3]. Zika virus was first detected in Uganda in 1947 ^[4]. According to the World Health Organization, about 86 countries were affected by Zika virus since the outbreak began ^[5]. In 2015 and during two years more than 4, 000 pregnant women were infected with Zika virus in Brazil which affect their new babies born ^[6,7].

The mathematical modeling in epidemiology began in the late 19th century and played an important role in studying, predicting, and proposing optimal control strategies for infectious diseases. A large number of mathematical models were proposed for a variety of infectious diseases ^[8,9,10,11]. In particular, several mathematical models predicting the transmission of Zika virus were proposed ^[12,13,14]. Many diseases prove seasonal comportment and thus taking account of seasonally in diseases modeling is important. For example, periodic fluctuations has the main impact in the evolution of disease transmissions which affect the contact rates that will change seasonally. Furthermore, periodic changes can affect birth rates of populations and thus, vaccination programs change seasonally. Variants of mathematical models are extensively used to model seasonally recurrent diseases. The mathematical models that describe these diseases are seasonally forced. Therefore, the seasonality of infectious diseases is very repetitive ^[15], and several mathematical models in epidemiology considering the impact of seasonality were analyzed ^{[14,16,17,18,19,20,21,22,23,24,25]}. When considering the seasonality in a mathematical model, the basic reproduction number can be approximated either trough the time-averaged model as in ^[26,27,28] or other ways as in ^{[29,30,31,32,33]}. In ^[34], the authors studied a periodic reaction-diffusion mathematical model for Zika virus transmission with seasonal and spatial heterogeneous structure, in ^[35], studied a partial differential equation model with periodic delay and in ^[36], the authors studied the impact of weather seasonality on the spread of Zika fever. Our aim is to consider the impact of the seasonality on the dynamics of ZIKV with a generalized incidence rate. The basic reproduction number, $\mathcal{R}_0$, was defined using the next generation matrix method in the case of the fixed environment and by using an integral linear operator for the case of seasonal environment. We perform the global analysis of the proposed system. It is deduced that the disease-free solution is globally asymptotically stable if $\mathcal{R}_0 < 1$. However, for the case where $\mathcal{R}_0 > 1$, we proved the persistence of disease. The theoretical findings were confirmed by intensive numerical example.

The structure of this manuscript is as follows. In the next Section, we describe a generalised compartmental model for ZIKV dynamics in a seasonal environment. In Section 3, we consider in the first step the case of a fixed environment, we calculate $\mathcal{R}_0$, and we study the local and global stability of the equilibria of the system. It is deduced that the disease-free steady state is stable if $\mathcal{R}_0 < 1$; however, the endemic steady state is stable if $\mathcal{R}_0 > 1$. In section 6, we focus on the influence of the seasonality. We prove that the virus-free periodic solution is stable if $\mathcal{R}_0 < 1$; however, the disease will persist if $\mathcal{R}_0 > 1$. We give in Section 7 several numerical tests confirming the theoretical results. We finish by giving some concluding remarks in section 8.

2.   Generalised Zika epidemic model

The ZIKV transmission follows the following steps. Mosquitoes get the virus when biting infected humans. Later, infected mosquitoes spread the ZIKV when biting uninfected humans. It should be noted that infected mosquitoes remain infected until they die. However, an infected human can recover and become immune against the disease. Thus, the model that we proposed here uses an SI compartmental model to predict the virus transmission in the mosquitoes population and a SIR-compartmental model to predict the virus spread within the human population ^[37]. Thus, the proposed model is a compartmental one generalizing the ones given in ^{[38,39,40,41]} and described by the following five dimensional dynamics of ordinary differential equations.

with the positive initial condition $(X_s^h(0), X_i^h(0), X_r^h(0), X_s^v(0), X_i^v(0)) \in \mathbb{R}^{5}_+$. The susceptible human are denoted by $X_s^h$, the infected human population are denoted by $X_i^h$ and the recovered human are denoted by $X_r^h$. Similarly, the susceptible mosquito are denoted by $X_s^v$ and the infected mosquito are denoted by $X_i^v$. More details on the meaning of the parameters are resumed in Table 1. The susceptible human catches up with the infection at a rate $\beta_h f_h(X_i^v)X_s^h$, with $\beta_h$ describing the contact rate of uninfected human and infected mosquito, and $f_h$ is the infected mosquito to uninfected human incidence rate. In the mosquito population, the susceptible mosquito catches up with the infection at a rate $\beta_v f_v(X_i^h)X_s^v$, where $\beta_v$ is the contact rate of uninfected mosquito and infected human, and $f_v$ is the infected human to uninfected mosquito incidence rate. The bilinear incidence rates in epidemiological models are intensively used ^[42]. When considering real data of disease dynamics, incidence rates are more appropriate with nonlinear forms ^[32].

We suppose that the parameters of the considered system are non-negative continuous bounded and $T$-periodic functions. We assume also that a susceptible human catches up with the infection only in the presence of an infected mosquito and similarly, a susceptible mosquito becomes infected only in the presence of an infected human and that transmission rates increase with the infected human and infected mosquitoes. Therefore, the model (2.1) satisfied the assumption given hereafter.

Assumption 1. (1) $f_h$ and $f_v$ are non-negative $C^1(\mathbb{R}_+)$, increasing concave functions satisfying $f_h(0) = f_v(0) = 0$.

(2) $\Lambda_h(t)$, $\Lambda_v(t)$, $\beta_h(t)$, $\beta_v(t)$, $m_h(t)$, $m_v(t)$, $r_h(t)$ and $u(t)$ are continuous, bounded and $T$-periodic non-negative functions.

Lemma 1. $Xf_h'(X) \leq f_h(X) \leq Xf_h'(0)$ and $Xf_v'(X) \leq f_v(X) \leq Xf_v'(0)$, $\forall X \in \mathbb{R}_+$.

Proof. For $X, X_1 \in \mathbb{R}_+$, let $g_1(X) = f_h(X)-Xf_h'(X)$. By using Assumption 1, we have $f_h'(X) \geq 0$ and $f_h''(X) \leq 0$. Then, $g_1'(X) = -Xf_h''(X) > 0$ and $g_1(X) \geq g_1(0) = 0$ which leads to $f_h(X)\geq Xf_h'(X)$. By the same way, let $g_2(X) = f_h(X)-Xf_h'(0)$ then $g_2'(X) = f_h'(X)-f_h'(0) < 0$ and $g_2(X) \leq g_2(0) = 0$ then $f_h(X)\leq Xf_h'(0)$. The proof is the same for the function $f_v$. □

3.   Case of autonomous system

We start by studying the case of constant parameters and thus we obtain the following system considered already in ^[41].

such that $(X_s^h(0), X_i^h(0), X_r^h(0), X_s^v(0), X_i^v(0)) \in \mathbb{R}^{5}_+$.

We begin by giving some basic properties of the system (3.1) as follows.

3.1.   Basic properties

Lemma 2. The dynamics (3.1) admits an invariant attractor set given by

Proof. Since $\dot{X_s^h}\mid_{X_s^h = 0} = m_h \Lambda_h > 0$, $\dot{X_i^h}\mid_{X_i^h = 0} = \beta_h f_h(X_i^v)X_s^h\geq 0$, $\dot{X_r^h}\mid_{X_r^h = 0} = (u+r_h) X_i^h \geq 0$, $\dot{X_s^v}\mid_{X_s^v = 0} = m_v \Lambda_v > 0$, and $\dot{X_i^v}\mid_{X_i^v = 0} = \beta_v f_v(X_i^h) X_s^v\geq 0$. Therefore, $\mathbb{R}_+^5$ is invariant by the model (3.1). Let us denote by $T_h = X_s^h+X_i^h+X_r^h$ and $T_v = X_s^v+X_i^v$ to be the sizes of the total human and mosquitoes populations, respectively. From Eq (3.1) we have $\dot T_h = m_h \Lambda_h - m_h T_h$. Hence $T_h = \Lambda_h$ if $T_h(0) = \Lambda_h$. Similarly, $\dot T_v = m_v \Lambda_v - m_v T_v$. Hence $T_v = \Lambda_v$ if $T_v(0) = \Lambda_v$. □

Let us now discuss the existence and uniqueness of equilibrium points of system (3.1).

3.2.   Steady states: existence and uniqueness

We start by calculating the basic reproduction number of our system (3.1) denoted by $\mathcal{R}_{0}$ ^[10,11]. We consider the matrices

and

Then,

and $\mathcal{R}_{0}$ is given by

Lemma 3. ● If $\mathcal{R}_0\leq 1,$ then the system (3.1) admits an equilibrium point denoted by $\mathcal{E}_{0} = \left(\Lambda_h, 0, 0, \Lambda_v, 0\right)$.

● If $\mathcal{R}_0 > 1,$ then the system (3.1) admits two steady states; $\mathcal{E}_0$ and an endemic equilibrium point denoted by $\bar{\mathcal{E}}$.

Proof. To prove the existence and uniqueness of the equilibria according to the values of the basic reproduction number, let $E(X_s^h, X_i^h, X_r^h, X_s^v, X_i^v)$ be any steady state satisfying

which is equivalent to

and

Now, using the second equation of (3.2), one has

If $X_i^h = 0$, then we obtain an equilibrium point given by the ZIKV-free equilibrium point $\mathcal{E}_0 = \left(\Lambda_h, 0, 0, \Lambda_v, 0\right)$. If $X_i^h\not = 0$, let us define the function $g$ as follows:

The limit of the function $g$ at the origin is

Note that the value of $g$ at $\Lambda_h$ is

Furthermore, the derivative of $g$ on $(0, \Lambda_h)$ is expressed as follows

Therefore, we deduce that the function $g$ is decreasing. Then $g$ has a unique root $\bar{X}_i^h\in(0, \Lambda_h)$. Therefore,

and the endemic steady state denoted by $\bar{\mathcal{E}} = (\bar{X}_s^h, \bar{X}_i^h, \bar{X}_r^h, \bar{X}_s^v, \bar{X}_i^v)$ exists if only if $\mathcal{R}_0 > 1$. □

4.   Local stability

We aim in this section to study the local stability of both equilibrium points $\mathcal{E}_0$ and $\bar{\mathcal{E}}$ with respect to the values $\mathcal{R}_0$.

Theorem 1. For $\mathcal{R}_0 < 1$, $\mathcal{E}_0$ is locally asymptotically stable (LAS).

Proof. The Jacobian matrix for $\mathcal{E}_0$ is

admitting the following three eigenvalues $\lambda_1 = \lambda_2 = - m_h < 0$ and $\lambda_3 = - m_v < 0$. By considering the sub-matrix

where the trace satisfies Trace$(S_{j_0}) = -(r_h+u+ m_h+ m_v) < 0$ and det$(S_{j_0}) = m_v(r_h+u+ m_h)-\beta_h \beta_v f_h'(0)f_v'(0) \Lambda_h \Lambda_v = m_v(r_h+u+m_h) (1-\mathcal{R}_0^2)$. Therefore $J_0$ admits four eigenvalues with negative real parts if $\mathcal{R}_0 < 1$ and then $\mathcal{E}_0$ is LAS for $\mathcal{R}_0 < 1$. □

Theorem 2. For $\mathcal{R}_0 > 1$, $\bar{\mathcal{E}} = (\bar{X}_s^h, \bar{X}_i^h, \bar{X}_r^h, \bar{X}_s^v, \bar{X}_i^v)$ is LAS.

Proof. By calculating the Jacobian matrix at $\bar{\mathcal{E}} = (\bar{X}_s^h, \bar{X}_i^h, \bar{X}_r^h, \bar{X}_s^v, \bar{X}_i^v)$, we obtain :

admitting the following characteristic polynomial:

where

Using the fact that

and

we obtains

By applying the Routh-Hurwitz criterion ^[43,44], we deduce easily that the eigenvalues have negative real parts (see ^[45,46] for an other application). Thus, $\bar{\mathcal{E}}$ is LAS. □

5.   Global stability

Theorem 3. If $\mathcal{R}_0 \leq1,$ then $\mathcal{E}_0$ is globally asymptotically stable (GAS).

Proof. Consider the Lyapunov function $U_{0}(X_s^h, X_i^h, X_r^h, X_s^v, X_i^v)$:

Clearly, $U_{0}(X_s^h, X_i^h, X_r^h, X_s^v, X_i^v) > 0$ for all $X_s^h, X_i^h, X_r^h, X_s^v, X_i^v > 0$ and $U_{0}\left(\Lambda_h, 0, 0, \Lambda_v, 0\right) = 0$. The time derivative of $U_{0}$ is :

If $\mathcal{R}_{0}\leq1,$ then $\dfrac{dU_{0}}{dt}\leq0$ for all $X_s^h, X_i^h, X_r^h, X_s^v, X_i^v > 0.$ Let

It can be easily shown that $W_{0} = \left\{\mathcal{E}_0\right\}$. Applying LaSalle's invariance principle ^[9], we deduce that $\mathcal{E}_0$ is GAS when $\mathcal{R}_0\leq 1.$ □

Define the set

Theorem 4. If $\mathcal{R}_0 > 1,$ then $\bar{\mathcal{E}} = (\bar{X}_s^h, \bar{X}_i^h, \bar{X}_r^h, \bar{X}_s^v, \bar{X}_i^v)$ is GAS in $\Gamma_2$.

Proof. Let us define the function $G(X) = X-1-\ln (X)$ which is a positive function defined on $\mathbb{R}_+^*$ with derivative $G'(X) = 1-\dfrac{1}{X}$. Consider the Lyapunov function denoted by $\bar U(X_s^h, X_i^h, X_r^h, X_s^v, X_i^v)$ and defined as:

Clearly, $\bar U(X_s^h, X_i^h, X_r^h, X_s^v, X_i^v) > 0$ for all $X_s^h, X_i^h, X_r^h, X_s^v, X_i^v > 0$ and $\bar U(\bar{X}_s^h, \bar{X}_i^h, \bar{X}_r^h, \bar{X}_s^v, \bar{X}_i^v) = 0.$ The time derivative of $\bar U$ is :

By using the fact that

we get

Therefore, $\dfrac{d\bar U}{dt}\leq 0$ for all $X_s^h, X_i^h, X_r^h, X_s^v, X_i^v \in \Gamma_2$ and $\dfrac{d\bar U}{dt} = 0$ if and only if $(X_s^h, X_i^h, X_r^h, X_s^v, X_i^v) = (\bar{X}_s^h, \bar{X}_i^h, \bar{X}_r^h, \bar{X}_s^v, \bar{X}_i^v) = 0$. Using the LaSalle's invariance principle ^[9], we obtain the global stability of $\bar{\mathcal{E}}$ in $\Gamma_2$. □

6.   Case of periodic environment

In this section, we return to the main system (2.1) studying the seasonality influence that we write it in the following way:

with positive initial condition $(X_i^h(0), X_i^v(0), X_s^h(0), X_r^h(0), X_s^v(0)) \in \mathbb{R}^{5}_+$. Let $\rho(t)$ to be a continuous, positive $T$-periodic function. Let us denote by $\rho^u = \max\limits_{t\in[0, T)} \rho(t)$ and $\rho^l = \min\limits_{t\in[0, T)} \rho(t)$.

6.1.   Preliminary

Let us consider the two-dimensional system

such that $(X_s^h(0), X_s^v(0)) \in \mathbb{R}_+^2$. System (6.2) admits exactly one $T$-periodic solution $(\bar X_s^h(t), \bar X_s^v(t))$ globally attractive in $\mathbb{R}_+^2$ with $\bar X_s^h(t) > 0$ and $\bar X_s^v(t) > 0$. Then, the main system (6.1) admits exactly one disease-free periodic solution $\mathcal{E}_0(t) = (0, 0, \bar X_s^h(t), 0, \bar X_s^v(t))$.

Proposition 1. The positive compact set

is an invariant and attractor of all solutions of model (6.1) such that

Proof. It is easy to see that

and $\dot X_s^v(t)+\dot X_i^v(t) = m_h(t)[\Lambda_v(t) -(X_s^v(t)+X_i^v(t))] \leq 0, \mbox{ if } X_s^v(t)+X_i^v(t) \geq \Lambda_v^u$. Let us define $B_1(t) = X_s^h(t)+X_i^h(t)+X_r^h(t)$ and $B_2(t) = X_s^v(t)+X_i^v(t)$. Consider $x_1(t) = B_1(t)-\bar X_s^h(t), t \geq 0$, then $\dot x_1(t) = -m_h(t) x_1(t)$, and therefore $\lim_{t \to \infty} x_1(t) = \lim_{t \to \infty} (B_1(t) - \bar X_s^h(t)) = 0$. Similarly, consider $x_2(t) = B_2(t)-\bar X_s^v(t), t \geq 0$, therefore $\dot x_2(t) = -m_v(t) x_2(t)$, and then $\lim_{t \to \infty} x_2(t) = \lim_{t \to \infty} (B_2(t) - \bar X_s^v(t)) = 0$. □

6.2.   Disease-free trajectory

In this section, we shall define the expression of the basic reproduction number; $\mathcal{R}_0$, according to the definition given by the theory in ^[32]. For $X = (X_i^h, X_i^v, X_s^h, X_r^h, X_s^v)$, let

and

and $\mathcal{V}(t, X) = \mathcal{V}^- (t, X)-\mathcal{V}^+(t, X)$. Therefore, the dynamics (6.1) can be written in the following way:

Then, it easy to see that conditions (A1)–(A5) of [32,Section 1] are valid.

The dynamics (6.4) admits a disease-free periodic trajectory $\bar X(t) = \mathcal{E}_0(t) = (0, 0, \bar X_s^h(t), 0, \bar X_s^v(t))$. Let us define

with $f_i(t, X(t))$ and $X_i$ are the $i$-th components of $f(t, X(t))$ and $X$, respectively. An easy calculus gives us

Therefore, $r(\beta_M(T)) < 1$ and the solution $\bar X(t)$ is linearly asymptotically stable in $\Omega_s = \left\{(0, 0, X_s^h, 0, X_s^v) \in R^5_+\right\}$. Therefore, the condition (A6) of [32,Section 1] holds.

Let us define ${\bf A^+}(t)$ and ${\bf A^-}(t)$ to be two matrices defined by

An easy calculus gives us

Consider $Z(s_1, s_2)$, the solution of the dynamics $\dfrac{d}{dt} Z(s_1, s_2) = -{\bf A^-}(s_1)Z(s_1, s_2)$ for any $s_1 \geq s_2$, with $Z(s_1, s_1) = I_2$. Thus, condition (A7) of [32,Section 1] is valid.

In order to define the basic reproduction number, $\mathcal{R}_0$, of (6.1), we define a linear integral operator as follows

where $C_T$ is the ordered Banach space of $T$-periodic functions defined on $\mathbb{R}$ to $\mathbb{R}^2$. Therefore,

Theorem 5. By using [32,Theorem 2.2], the following statements are verified:

● $\mathcal{R}_0 < 1 \; \Leftrightarrow \; r(\beta_{F-V} (T)) < 1$.

● $\mathcal{R}_0 = 1 \; \Leftrightarrow \; r(\beta_{F-V} (T)) = 1$.

● $\mathcal{R}_0 > 1 \; \Leftrightarrow \; r(\beta_{F-V} (T)) > 1$.

We deduce that $\mathcal{E}_0(t)$ is asymptotically stable if $\mathcal{R}_0 < 1$ and it is unstable if $\mathcal{R}_0 > 1$. Now, we show that if $\mathcal{R}_0 < 1$ then $\mathcal{E}_0(t) = (0, 0, \bar X_s^h (t), 0, \bar X_s^v (t))$ is globally asymptotically stable and thus the disease is extinct.

Theorem 6. $\mathcal{E}_0(t)$ is globally asymptotically stable for $\mathcal{R}_0 < 1$, however, it is unstable for $\mathcal{R}_0 > 1$.

Proof. Since Theorem 5 affirms that $\mathcal{E}_0(t)$ is locally stable for $\mathcal{R}_0 < 1$ and that it is unstable for $\mathcal{R}_0 > 1$, we need to prove the global attractivity for $\mathcal{R}_0 < 1$. We obtained the limits (6.3) in Proposition 1; therefore for $\kappa_1 > 0$, there exists a time $T_1 > 0$ satisfying $X_s^h(t)+X_i^h(t)+X_r^h(t) \leq \bar X_s^h(t)+ \kappa_1$ and $X_s^v(t)+X_i^v(t) \leq \bar X_s^v(t)+ \kappa_1$ for $t > T_1$. Therefore, $X_s^h(t) \leq \bar X_s^h(t)+ \kappa_1$ and $X_s^v(t)\leq \bar X_s^v(t)+ \kappa_1$; and

for $t > T_1$. Let us define the matrix $M_2(t)$ as follows

Since $r(\varphi_{F-V} (T)) < 1$, we can chose $\kappa_1 > 0$ small enough such that $r(\varphi_{F-V+ \kappa_1 M_2}(T)) < 1$. Consider now the following two-dimensional system,

From the theory in ^[47], there exists a positive function $x_1(t)$ that it is $T$-periodic satisfying $w(t)\leq x_1(t) e^{a_1t}$ where $w(t) = \left(X_i^h(t), X_i^v(t)\right)^T$ and $a_1 = \frac{1}{T} \ln\left(r(\varphi_{F-V+\kappa_1 M_2}(T)\right) < 0$. Thus, $\lim\limits_{t\to \infty} X_i^h(t) = 0$ and $\lim\limits_{t\to \infty} X_i^v(t) = 0$. Furthermore, we have that $\lim\limits_{t\to \infty} X_s^h(t)-\bar X_s^h(t) = \lim\limits_{t\to \infty} Z_1(t)-X_s^h(t)-\bar X_i^h(t)-\bar X_r^h(t) = 0$ and $\lim\limits_{t\to \infty} X_s^v(t)-\bar X_s^v(t) = \lim\limits_{t\to \infty} Z_2(t)-X_i^v(t)-\bar X_s^v(t) = 0$. Therefore, $\mathcal{E}_0(t)$ satisfies the globally attractivity for $\mathcal{R}_0 < 1$. □

Now, we show that if $\mathcal{R}_0 > 1$ then $X_i^h(t)$ and $X_i^v(t)$ are uniformly persistent and then the disease persists in the population.

6.3.   The endemic solution

Let us define $X_0 = (X_i^h(0), X_i^v(0), X_s^h(0), X_r^h(0), X_s^v(0))$ and $X_1 = (0, 0, \bar X_s^h(0), 0, \bar X_s^v(0))$ and consider the Poincaré map $Q : \mathbb{R}^5_+ \to \mathbb{R}^5_+$ associated with the model (6.1) such that $X_0 \mapsto u(T, X^0)$ is the solution of system (6.1) with the initial value $u(0, X^0) = X^0\in \mathbb{R}_+^5$. Consider the sets

and

Note that $Q$ is point dissipative. Furthermore, $\Omega$ and $\Omega_0$ are invariant. Through the theory in ^[8,47], we obtain

with $M_\partial\supseteq \left\{(0, 0, X_s^h, 0, X_s^v), \; X_s^h \geq 0, X_s^v \geq 0 \right\}$. It remains to prove that $M_\partial \setminus \left\{(0, 0, X_s^h, 0, X_s^v), \; X_s^h \geq 0, X_s^v \geq 0\right\} = \emptyset$. Consider $(X_0) \in M_\partial \setminus \left\{(0, 0, X_s^h, 0, X_s^v), \; X_s^h \geq 0, X_s^v \geq 0\right\}$.

If $X_i^v(0) = 0$ and $0 < X_i^h(0)$, then $\dot X_i^v(t)_{|t = 0} = \beta_v(t) f_v(X_i^h(0)) X_s^v(0) > 0$. If $X_i^v(0) > 0$ and $X_i^h(0) = 0$, therefore $X_i^v(t), X_s^h(t) > 0$ for any $t > 0$. Then, $\forall\; t > 0$, we have

which implies that $X(t) \not\in \partial \Omega_0$ for $0 < t \ll 1$ and that $\Omega_0$ is positively invariant and thus the satisfaction of (6.9). Therefore, $Q$ admits a fixed point $X_1$ in $M_\partial$. We obtain the following result.

Theorem 7. If $\mathcal{R}_0 > 1$, then the system (6.1) has at least a periodic trajectory satisfying $\exists \; \eta > 0$ such that $\forall \; X_0 \in Int(\mathbb{R}_+)^2 \times \mathbb{R}_+^3$ and $\liminf\limits_{t\to\infty} X_i^h(t) \geq \eta > 0, \liminf\limits_{t\to\infty} X_i^v(t) \geq \eta > 0.$

Proof. The goal is to prove the trajectories of (6.1) are uniformly persistent with respect to $(\Omega_0, \partial \Omega_0)$ using the properties of the Poincaré map, $Q$ as in [8,Theorem 3.1.1]. Since $r(\varphi_{F-V}(T)) > 1$, then $\exists\; \varepsilon > 0$ satisfying $r(\varphi_{F-V-\varepsilon M_2} (T)) > 1$. Let consider the following two-dimensional system

The Poincaré map, $Q$ related to the system (6.10) has a unique fixed point $(\bar X_{s\gamma}^h, \bar X_{s\gamma}^v)$ that it is globally attractive. By using the implicit function theorem, $\gamma \mapsto (\bar X_{s\gamma}^h, \bar X_{s\gamma}^v)$ is continuous. Thus, $\gamma > 0$ can be chosen small enough such that $\bar X_{s\gamma}^h(t) > \bar X_{s}^h(t)-\varepsilon, \;$ and $\bar X_{s\gamma}^v(t) > \bar X_{s}^v(t)-\varepsilon, \;$ $\forall \; t > 0$. Since the solution is continuous with respect to $X_0$, then there exists $\gamma^* > 0$ satisfying $\|X_0-u(t, X_1)\|\leq \gamma^*$; therefore

We aim to prove that

by contradiction. Assume that $\limsup\limits_{p \to \infty} d(Q^p(X_0), X_1) < \gamma^*$ for some $X_0\in \Omega_0$. We can assume that $d(Q^p(X_0), X_1) < \gamma^*, \forall\; p > 0$. Then we obtain

Assume that $t\geq 0$ can be written as $t = pT +t_1$ with $t_1 \in [0, T)$ and $p = \lfloor\dfrac{t}{T}\rfloor$. Therefore

Set $(X_i^h(t), X_i^v(t), X_s^h(t), X_r^h(t), X_s^v(t)) = u(t, X_0)$. Therefore $0 \leq X_i^h(t), X_i^v(t) \leq \gamma, t \geq 0$ and

The Poincaré map, $Q$ associated with the system (6.10) has a fixed point $(\bar X_{s\gamma}^h, \bar X_{s\gamma}^v)$ which is globally attractive where $\bar X_{s\gamma}^h(t) > \bar X_s^h-\varepsilon, \;$ and $\bar X_{s\gamma}^v(t) > \bar X_s^v(t)-\varepsilon$; then, $\exists \; T_2 > 0$ satisfying $X_s^h(t) > \bar X_s^h(t)-\varepsilon$ and $X_s^v(t) > \bar X_s^v(t)-\varepsilon$ for $t > T_2$. Then, for $t > T_2$, we have

Since $r(\varphi_{F-V-\varepsilon M_2} (T)) > 1$, then there exists a $T$-periodic positive function $x(t)$ ^[47] satisfying $J(t) \geq e^{a t}x(t)$ with $a = \dfrac{1}{T}\ln\; r\left(\varphi_{F-V-\varepsilon M_2} (T)\right) > 0$, thus $\lim\limits_{t\to \infty} X_i^h(t) = \infty$ which is impossible since the solution is bounded. Therefore, (6.11) is satisfied and $Q$ is weakly uniformly persistent with respect to $(\Omega_0, \partial \Omega_0)$. Regarding Proposition 1, the Poincaré map, $Q$ admits a global attractor. Therefore $X_1$ is an isolated invariant set of $\Omega$ and $W^s(X_1) \cap \Omega_0 = \emptyset$. Thus any solution in $M_\partial$ should converge to $X_1$ which is an acyclic in $M_\partial$. Applying [8,Theorem 1.3.1 and Remark 1.3.1], we deduce that $Q$ is uniformly persistent with respect to $(\Omega_0, \partial \Omega_0)$. Moreover, using [8,Theorem 1.3.6], $Q$ has a fixed point $\tilde X_0 = (\tilde X_i^h, \tilde X_i^v, \tilde X_s^h, \tilde X_r^h, \tilde X_s^v) \in \Omega_0$ with $\tilde X_0 \in Int(R_+)^2 \times R_+^3$.

Assume that $\tilde X_s^h = 0$ and by inject this in system (6.1), $\tilde X_s^h(t)$ verifies

with $\tilde X_s^h = \tilde X_s^h(nT) = 0, n = 1, 2, 3, \cdots$. From Proposition 1, $\forall\; \kappa_3 > 0$, $\exists\; T_3 > 0$ such that $\tilde X_i^v(t) \leq \Lambda_v^u + \kappa_3$ for $t > T_3$. Therefore, we get

$\exists \; \bar n > 0$ satisfying $nT > T_3, \forall\; n > \bar n$. Then we obtain

for $n > \bar n$ which is impossible since $\tilde X_s^h(nT) = 0$. Therefore, $\tilde X_s^h(0) > 0$ and then $\tilde X_0$ is a $T$-periodic positive solution of system (6.1). □

7.   Numerical investigation

Our objective of this section is to present some numerical simulations regarding the proposed mathematical model (2.1) that consider the influence of periodicity on the dynamics of the Zika virus. This model is a five dimensional compartmental model considering the dynamics of a population consisting of susceptible humans, infected humans, recovered humans, susceptible and infected mosquitoes. Several numerical illustrations will be used to exemplify the suitability and utility of the proposed Zika virus structure in the seasonal environment. All numerical simulations were done using the MATLAB R2024a software.

We used Monod-type functions for modelling both incidence rates:

where $\zeta_h$ and $\zeta_v$ are nonnegative constants. Note that the functions $f_h$ and $f_v$ are continuous, increasing and concave. Many diseases prove seasonal comportment and thus taking account of seasonally in diseases modeling is important. Variants of mathematical models are extensively used to model seasonally recurrent diseases. Seasonality may come from various sources. A famous example of a seasonally forced function can take the following form $k(t) = k_0(1 +k_1\cos(2 \pi(t+\psi)))$, where $k_0 \geq 0$ is the baseline transmission parameter, $0 < k_1 \leq 1$ is the amplitude of the seasonal variation in transmission and $0 \leq \psi \leq 1$ is the phase angle. Therefore, for all numerical simulations, the periodic functions that reflect the influence of seasonality on the dynamics of the Zika virus dynamics are given by

The parameter values employed to generate the figures in this section are as follows: The constants $\zeta_h$, $\zeta_v$, $\Lambda_h^0$, $\Lambda_v^0$, $\beta_h^0$, $\beta_v^0$, $m_h^0$, $m_v^0$, $r_h^0$ and $u^0$ and the phase shift $\psi$ are given in Table 2. Unfortunately, we have no biological data to use for our simulations. Parameters values considered here have no biological meaning and are chosen arbitrarily.

The seasonal cycles frequencies $\Lambda_h^1$, $\Lambda_v^1$, $\beta_h^1$, $\beta_v^1$, $m_h^1$, $m_v^1$, $r_h^1$ and $u^1$ are displayed in Table 3.

Numerical investigations of the considered model are discussed for three cases, namely, autonomous system (the parameters are assumed to be constants), periodic contact between human and mosquito (where only contact rates are assumed to be periodic functions with same period) and full seasonal environment (where all parameters are assumed to be periodic functions with same period).

7.1.   Autonomous system

The numerical examples given in this subsection deal with the case of autonomous system with fixed parameters.

with $(X_s^h(0), X_i^h(0), X_r^h(0), X_s^v(0), X_i^v(0)) \in \mathbb{R}^{5}_+$. We calculated $\mathcal{R}_0$ using the next generation matrix method ^[10,11].

In Figure 1 we present the numerical simulations of the model (7.2) for two values of the basic reproduction number. As it can be seen, the solutions of the system (7.2) converge to endemic equilibrium point, $\bar{\mathcal{E}} = (\bar{X}_s^h, \bar{X}_i^h, \bar{X}_r^h, \bar{X}_s^v, \bar{X}_i^v)$, reflecting the persistence of disease when $\mathcal{R}_0 > 1$ (left), however, it converges asymptotically to the disease-free equilibrium point $\mathcal{E}_0 = (\Lambda_h^0, 0, 0, \Lambda_v^0, 0)$ for the case where $\mathcal{R}_0\leq 1$ (right). In Figures 2 and 3, we consider several initial conditions where all corresponding solutions converge to the same equilibrium point for both cases of the $\mathcal{R}_0$ values. Therefore, Figures 2 and 3 confirm the global stability of $\mathcal{E}_0$ and $\bar{\mathcal{E}}$ for the cases $\mathcal{R}_0\leq 1$ and $\mathcal{R}_0 > 1$, respectively.

7.2.   Case of seasonal contact

The numerical examples given in this subsection deal with the case where only contact between humans and mosquitoes is assumed to be seasonal and then the contact rates, $\beta_h$ and $\beta_v$ are periodic functions.

with $(X_s^h(0), X_i^h(0), X_r^h(0), X_s^v(0), X_i^v(0)) \in \mathbb{R}^{5}_+$. We calculated $\mathcal{R}_0$ using the time-averaged dynamics as in ^[26,27]. Several other approximations of $\mathcal{R}_0$ are used in ^[30,31]. In Figure 4, the solutions of the system (7.3) converge to a periodic orbit reflecting the persistence of disease when $\mathcal{R}_0 > 1$ (left), however, it converges asymptotically to the periodic solution $\mathcal{E}_0(t) = (\bar X_s^h (t), 0, 0, \bar X_s^v (t), 0)$ for the case where $\mathcal{R}_0\leq 1$ (right). In Figures 5 and 6, we consider several initial conditions where all corresponding solutions converge to the same periodic solution for both cases of the $\mathcal{R}_0$ values. Therefore, Figures 5 and 6 confirm the global stability of $\mathcal{E}_0(t) = (\bar X_s^h(t), 0, 0, \bar X_s^v(t), 0)$ and the persistence of the disease for the cases $\mathcal{R}_0\leq 1$ and $\mathcal{R}_0 > 1$, respectively.

7.3.   Full periodic system

The numerical examples given in this subsection deal with the full seasonal environment where all the parameters of the system are periodic functions.

with $(X_s^h(0), X_i^h(0), X_r^h(0), X_s^v(0), X_i^v(0)) \in \mathbb{R}^{5}_+$. We calculated $\mathcal{R}_0$ using the time-averaged dynamics as in ^[26,27]. In Figure 7, the solutions of the system (7.4) converge to a periodic orbit reflecting the persistence of disease when $\mathcal{R}_0 > 1$ (left), however, it converges asymptotically to the periodic solution $\mathcal{E}_0(t) = (\bar X_s^h (t), 0, 0, \bar X_s^v (t), 0)$ for the case where $\mathcal{R}_0\leq 1$ (right). In Figures 8 and 9, we consider several initial conditions where all corresponding solutions converge to the same periodic solution for both cases of the $\mathcal{R}_0$ values. Therefore, Figures 8 and 9 confirm the global stability of $\mathcal{E}_0(t) = (\bar X_s^h(t), 0, 0, \bar X_s^v(t), 0)$ and the persistence of the disease for the cases $\mathcal{R}_0\leq 1$ and $\mathcal{R}_0 > 1$, respectively.

8.   Conclusions

In this research, we devised a reliable Zika virus model considering the impact of seasonality observed in real life. The qualitative analysis of this model is presented in both cases, fixed and seasonal environments. We calculated the basic reproduction number using two different methods, the next generation matrix method in the case of the fixed environment and through a linear integral operator in the case of the seasonal environment. Therefore, we investigated the local and global stability for both cases. It is deduced that if $\mathcal{R}_0 \leq 1$, trajectories of the system approach a disease-free periodic solution and then the disease goes extinct; however, if $\mathcal{R}_0 > 1$, the disease persists and the trajectories of the system converge to a limit cycle. In our case, the solution of the considered model in a seasonal case shows a periodic behavior with an average close to the solution of the model in a fixed environment. This means that the main difference between the autonomous system and the periodic environment case is qualitative.

Author contributions

M. El Hajji: Conceptualization, Methodology, Writing-original draft, Writing-review and diting, Supervision; M. F. S. Aloufi: Conceptualization, Methodology, Writing-original draft, Writing-review and editing; M. H. Alharbi: Conceptualization, Methodology, Writing-original draft, Writing-review and editing, Supervision. All authors have read and approved the final version of the manuscript for publication.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgements

The authors are grateful to the unknown referees for many constructive suggestions, which helped to improve the presentation of the paper.

Conflict of interest

All the authors declare no conflicts of interest.