Periodic solutions for chikungunya virus dynamics in a seasonal environment with a general incidence rate

: The chikungunya virus (CHIKV) infects macrophages and adherent cells and it can be transmitted via a direct contact with the virus or with an already infected cell. Thus, the CHIKV infection can have two routes. Furthermore, it can exhibit seasonal peak periods. Thus, in this paper, we consider a dynamical system model of the CHIKV dynamics under the conditions of a seasonal environment with a general incidence rate and two routes of infection. In the ﬁrst step, we studied the autonomous system by investigating the global stability of the steady states with respect to the basic reproduction number. In the second step, we establish the existence, uniqueness, positivity and boundedness of a periodic orbit for the non-autonomous system. We show that the global dynamics are determined by using the basic reproduction number denoted by R 0 and they are calculated using the spectral radius of an integral operator. We show the global stability of the disease-free periodic solution if R 0 < 1 and we also show the persistence of the disease if R 0 > 1 where the trajectories converge to a limit cycle. Finally, we display some numerical investigations supporting the theoretical ﬁndings.


Introduction
The chikungunya virus (CHIKV) is an arbovirus (i.e., a virus transmitted by arthropods) whose vectors are female mosquitoes of the genus Aedes which are identifiable by the presence of black and white stripes. The two implicated species are Aedes aegypti and Aedes albopictus. Aedes albopictus is present in the south of France and Aedes aegypti in areas including Antilles, Guyana, French Polynesia and New Caledonia. These two mosquitoes are also implicated in the transmission of other arboviruses, including dengue, yellow fever and Zika virus. This disease has been prevalent on the African and Asian continents for more than 50 years. Over the past 30 years, mathematical modeling has been applied to study several epidemic diseases [1][2][3][4]. Recently, the study of vector-borne diseases has gained considerable attention and mathematics have become a useful tool for such studies; also, several temporal deterministic models have been proposed for diseases like dengue, malaria, CHIKV, etc. Several mathematical models have been developed and studied to explain a variety of features influencing the transmission of CHIKV [5][6][7][8][9][10][11][12].
However, several infectious diseases including all diseases caused by the transmission of a pathogenic agent exhibit seasonal peak periods. Studying the population behaviors associated with a seasonal environment becomes a necessity for predicting the risk of transmission of such a disease . In [13], the authors discuss the periodic "SIR" epidemic model. In [14], the authors study an "SEIRS" epidemic model with periodic fluctuations. They calculated the basic reproduction number R 0 by using the time-averaged system and proved a sufficient but unnecessary condition such that the disease could not persist. In [15], the authors consider a class of SIQRS models with periodic behavior of the contact rate; they proved the existence of periodic trajectories. The authors of [16,17] studied an "SEIRS" epidemic model in a seasonal environment and proved some sufficient conditions for both the persistence and the extinction of the disease. In [18], the authors give the definition of the basic reproduction number for seasonal environments. The basic reproduction number R 0 is defined in [19] for several compartmental periodic epidemic models; the authors show that R 0 is a threshold value to prove the stability of the disease-free periodic solution. In [20], the authors propose an extension of the "SVEIR" model by taking into account the seasonal environment.
The aim of this work is to propose a new class of dynamical systems that models CHIKV dynamics under the conditions of a seasonal environment with a general incidence rate, and where the adherent cells are the main target for CHIKV. We will establish the existence, uniqueness, positivity and boundedness of a periodic solution. Therefore, we will study the global dynamics with respect to the basic reproduction number that will be calculated by using the spectral radius of an integral operator. The global stability of the disease free periodic solution will be proved for R 0 < 1; however, the persistence of the disease will be proved for R 0 > 1 by proving that the trajectories will converge to a limit cycle. Finally, some numerical simulations will be given to confirm the theoretical results.

Proposed mathematical model
CHIKV is spread as follows. Mosquitoes contract the virus by biting animals or humans infected with it. They then spread the virus by biting uninfected people. CHIKV infects macrophages and adherent cells and it can be transmitted via a direct contact with the virus or with an already infected cell. Thus, the CHIKV infection can have two routes, i.e., CHIKV-to-cell and cell-to-cell infections (see Figure 1). We adopt a general non-linear incidence rate for both routes of infection. Furthermore, CHIKV dynamics can exhibit seasonal peak periods, which is why all parameters of the models are T -periodic functions where T > 0 is the period. The considered epidemic model for the CHIKV dynamics in a seasonal environment with a general form of the incidence rate is given by the following four-dimensional ordinary differential equation system which generalizes the models given in [21,22].
with the positive initial condition (S 0 , I 0 , P 0 , A 0 ) ∈ R 4 + . S (t), I(t), P(t), and A(t) describe susceptible cells, infected cells, CHIKV and antibodies, respectively. The uninfected cells are generated by a periodic rate m(t)S in (t), die at a periodic rate m(t)s(t) and become infected by the virus and infected cells at a periodic rate β 1 (t)P(t) f (S (t)) + β 2 (t)I(t) f (S (t)), where β 1 (t) and β 2 (t) are the periodic incidence rates. The periodic variables m(t), m p (t), and m a (t) represent, respectively, the periodic mortality rates for the infected cells, CHIKV and antibodies. δ(t) is the rate of periodic production of CHIKV from infected cells. Antibodies attack the CHIKV at a periodic rate of r(t)A(t)P(t). Once an antigen is encountered, the antibodies expand at a periodic rate of m a (t)A in (t) and proliferate at a periodic rate of k(t)r(t)A(t)P(t). All of the parameters of the model are positive periodic functions. We give in Table 1 more epidemiological significance of the model parameters.
Diagram describing an epidemic compartmental model that takes into consideration a seasonal environment for the CHIKV dynamics (inspired from [23, Figure  2]). Compartments S , I, P and A are described by circles; transition rates between compartments are described by arrows and labels.
Throughout the rest of the paper, we will use the following assumptions: (1) f is an increasing, non-negative C 1 (R + ) concave function such that f (0) = 0.
Assumption 1 expresses that the CHIKV-to-cell and cell-to-cell incidence rates increase with increasing number of susceptible cells. Assumption 1 affirms also that no CHIKV-to-cell or cellto-cell infection can take place in the absence of susceptible cells. All of the model parameters are T -periodic functions influenced by the seasonal environment. We assume also that the instantaneous CHIKV mortality rate is greater than the instantaneous antibody loss rate.
Proof. Let x, x 1 ∈ R + , and the function g Rate of instantaneous production of the virus from infected cells r(t) Rate of instantaneous attack of the virus by the antibodies k(t) Instantaneous proliferation rate for antibodies β 1 (t), β 2 (t) Instantaneous incidence coefficients 3. Mathematical analysis for the autonomous system As a first step, we consider the autonomous form of the model (2.1) where its right-hand side is independent of t i.e., all parameters are constants. In this section, we will use the following additional assumption on the incidence rate.
The given condition of Assumption 2 is a mathematical artifice and has no biological meaning; here, it is only used to prove the existence and uniqueness of the positive steady state.

Basic properties
In this subsection, we give some basic properties for the model (3.1) such as the existence, positivity and boundedness of the trajectories of (3.1), as well as the existence of an invariance set of all solutions of system (3.1). Proof. Note that if S = 0 thenṠ = mS in > 0; if Hence Then Thus, Ω is invariant for the model (2.1) since all variables are non-negative.

Basic reproduction number and steady states
When studying a disease, an important health question is whether the disease is spreading in the population. The response can be obtained by calculating the average number of people an infectious person infects while they are contagious. This average is known as the basic reproduction number whose calculation is complex. It can be determined by using the values of the model parameters; thus, we can determine if the disease spreads. The basic reproduction number can be calculated by using several methods. For example, in the case of a single infected compartment, the basic reproduction number is simply the product of the infection rate and its mean duration [24]. It can be calculated by using graph, or network, theory [25,26]. If there are several compartments representing infectious individuals as in our case here, the next-generation matrix method introduced by Diekmann et al. [24] and developed later in [27,28], divides the population into two compartments whose first compartment represents the infected individuals. The goal is therefore to see the rate of change in the population established in each of these compartments. The matrix approach calculating the basic reproduction number explains the relationship between compartmental models and population matrix models. For the dynamics given by (3.1), we shall calculate the basic reproduction number by using the next generation matrix method. Therefore, we will calculate the steady states by proving their existence and uniqueness.
The determinant of V is given by det(V) = m a m(rA in + m p ) > 0; thus, and the next-generation matrix is given by Therefore, the basic reproduction number (i.e., spectral radius of FV −1 ) is given by: .
Proof. Let E = (S , I, P, A) be an equilibrium point satisfying the following , We define the function .
Then, we obtain One deduces, therefore, that The derivative of the function g is given by (3.8) By Assumption 1, we have that g (P) ≤ 0 ∀ P ∈ (0, m a kr ). Then, the function g(P) admits a unique root P * ∈ (0, m a kr ). Thus, we obtain Thus, the infected equilibrium E * = (S * , I * , P * , A * ) exists if R 0 > 1.
From the equilibrium point conditions of E * , we have that mS in = mS * + β 1 P * f (S * ) + β 2 I * f (S * ) → mS * + mI * = mS in ; then, S * + I * = S in . Furthermore, we have that m p kP * = δkI * + m a A in − m a A * ; then. m a kP * +m a A * ≤ m p kP * +m a A * = δkI * +m a A in < δkS in +m a A in , which means that kP * +A * < A in + δkS in m a .

Local stability
In this subsection, we aim to study the local stability of the steady states of system (3.1) by using the linearization method with the Jacobian matrix.
Proof. The Jacobian matrix at point E 0 is given by: J 0 admits four eigenvalues; λ 1 = −m < 0 and λ 2 = −m a < 0. λ 3 and λ 4 are eigenvalues of the sub-matrix The trace of S j 0 is given by and the determinant of S j 0 is given by Then, E 0 is locally asymptotically stable if R 0 < 1, and it is unstable if R 0 > 1.
Proof. The Jacobian matrix at a point E * = (S * , I * , P * , A * ) is given by: The characteristic polynomial is then given by: The characteristic polynomial P(X) = 0 if, and only if , the left-hand side satisfies the following condition: and the right-hand side satisfies the following condition: This is impossible; then, Re(X) < 0 and E * is locally asymptotically stable.

Global stability
In this subsection, we aim to study the global stability of the steady states of system (3.1) by using the Lyapunov theory. Let us define the function G(z) = z − 1 − ln z that will be used throughout this subsection.

Periodic system
Let us reconsider the periodic dynamics given by (2.1) where our aim is to prove that the system admits a bounded positive T -periodic solution.
For a continuous, positive T -periodic function g(t), we set g u = max t∈[0,T ) g(t) and g l = min t∈[0,T ) g(t).

Preliminary
Let (R m , R m + ) be the ordered m-dimensional Euclidean space associated with the norm · . For . Consider a T -periodic m×m matrix function denoted by C(t) which is continuous, irreducible and cooperative. Let us denote by φ C (t) the fundamental matrix, which is the solution of the following systeṁ (4.1) Let us denote the spectral radius of the matrix φ C (T ) by r(φ C (T )). Therefore, all entries of φ C (t) are positive for each t > 0. Let us apply the theorem of Perron-Frobenius to deduce that r(φ C (T )) is the principal eigenvalue of φ C (T ) (simple and admits an eigenvector y * 0). For the rest of the paper, the following lemma will be useful. Lemma 4. [30]. There exists a positive T -periodic function y(t) such that x(t) = y(t)e kt is a solution of system (4.1) where k = 1 T ln(r(φ C (T ))).
Let us start by proving the existence (and uniqueness) of the disease free periodic trajectory of model (2.1). Let us consider the following subsysteṁ with the initial condition (S 0 , A 0 ) ∈ R 2 + . Equation (4.2) admits a unique T -periodic solution (S * (t), A * (t)) with S * (t) > 0 and A * (t) > 0 which is globally attractive in R 2 + ; thus, system (2.1) has a unique disease-free periodic solution (S * (t), 0, 0, A * (t)).
Let us introduce the following result.

(t)Ṗ(t) +Ȧ(t) =k(t)(δ(t)I(t) − m p (t)P(t)) + m a (t)A in (t) − m a (t)A(t) =k(t)δ(t)I(t) − k(t)m p (t)P(t) + m a (t)A in (t) − m a (t)A(t)
which implies that Ω u is a forward invariant compact absorbing set for (2.1). Let N 1 (t) = S (t) + I(t) and N 2 (t) = k(t)P(t) + A(t) be the sub-population sizes at time t. Next, let y 1 (t) = N 1 (t) − S * (t), t ≥ 0. Then, it follows thatẏ 1 (t) = −m(t)y 1 (t), which implies that lim Next, in Subsection 4.2, we define R 0 , the basic reproduction number and we will prove that the disease free periodic trajectory (0, 0, S * (t), A * (t)) is globally asymptotically stable (and therefore, that the disease dies out) once R 0 < 1. Then, in Subsection 4.3, we will prove that I(t) and P(t) exhibit uniform persistence (i.e., the disease persists) once R 0 > 1. Therefore, we deduce that R 0 is the threshold parameter between the uniform persistence and the extinction of the disease.

Disease free periodic solution
We start by giving the definition of the basic reproduction number for model (2.1) by using the theory given in [19] where Our aim is to check the conditions (A1)-(A7) in [19, Section 1]. Note that system (2.1) can have the following formẊ (4.5) The first five conditions (A1)-(A5) are fulfilled.
The system (4.5) admits a disease free periodic trajectory X * (t) = where f i (t, X(t)) and X i are the ith components of f (t, X(t)) and X, respectively. By an easy calculation, we get that −m a (t) and then that r(φ M (T )) < 1. Therefore X * (t) is linearly asymptotically stable in the subspace Γ s = (0, 0, S , A) ∈ R 4 + . Thus, the condition (A6) in [19, Section 1] is satisfied.
Now, let us define F(t) and V(t) to be two by two matrices given by and where F i (t, X) and V i (t, X) are the i-th components of F (t, X) and V(t, X), respectively. By an easy calculation, we obtain the following from system (4.5): .
Consider Z(t 1 , t 2 ) to be the two by two matrix solution of the system d dt Z(t 1 , t 2 ) = −V(t 1 )Z(t 1 , t 2 ) for any t 1 ≥ t 2 , with Z(t 1 , t 1 ) = I, i.e., the two by two identity matrix. Thus, condition (A7) is satisfied.
Let us define C T to be the ordered Banach space of T -periodic functions defined on R → R 2 , associated with the maximum norm . ∞ and the positive cone C + T = {ψ ∈ C T : ψ(s) ≥ 0, for any s ∈ R}. Define the linear operator K : Let us now define the basic reproduction number, R 0 , of model (2.1) by using R 0 = r(K). Therefore, we conclude the local asymptotic stability of the disease free periodic solution E 0 (t) = (S * (t), 0, 0, A * (t)) for (2.1) to be as follows.
Therefore, E 0 (t) is unstable if R 0 > 1 and it is asymptotically stable if R 0 < 1.
Proof. Using Theorem 5, we have that E 0 (t) is locally stable once R 0 < 1 and it is unstable once R 0 > 1. Therefore, we need to prove the global attractivity of E 0 (t) when R 0 < 1.
For the following subsection, we consider only the case in which R 0 > 1.

Case in which all parameters are constants
In the first step, we performed numerical simulations for the system (3.1) when all parameters are constant. Thus, the model is given by with the positive initial condition (S 0 , E 0 , I 0 , R 0 ) ∈ R 4 + . We give the results of some numerical simulations confirming the stability of the steady states of system (5.2).
In Figure 2, the approximated solution of system (5.2) approaches asymptotically to E * once R 0 > 1. In Figure 3, the approximated solution of the given model (5.2) approaches the equilibrium E 0 , which confirms that E 0 is globally asymptotically stable once R 0 ≤ 1.

Case in which all parameters are constant with a periodic seasonally forced function
In the second step, we performed numerical simulations for the system (2.1) where only the Tperiodic seasonally forced functions β 1 and β 2 are dependent on time. The other parameters were set to be constant. Thus the model is given by with the positive initial condition (S 0 , I 0 , P 0 , A 0 ) ∈ R 4 + . We give the results of some numerical simulations confirming the stability of the steady states of system (5.3). The basic reproduction number R 0 was approximated by using the time-averaged system.
In Figure 4, the approximated solution of system (5.3) asymptotically approaches the periodic solution with the persistence of the disease. In Figure 5, we show a magnified view of the limit cycle when R 0 > 1. In Figure 6, the approximated solution of system (5.3) approaches the disease-free trajectory once R 0 < 1.

Case in which all parameters are periodic functions
In the third step, we performed numerical simulations for the system (2.1) where all parameters were set as T -periodic functions. Thus the model is given by I(t) = ηβ 1 (t)P(t)S (t) κ + S (t) + ηβ 2 (t)I(t)S (t) κ + S (t) − m(t)I(t), P(t) = δ(t)I(t) − m p (t)P(t) − r(t)A(t)P(t), A(t) = m a (t)A in (t) + k(t)r(t)A(t)P(t) − m a (t)A(t), (5.4) with the positive initial condition (S 0 , I 0 , P 0 , A 0 ) ∈ R 4 + . We give the results of some numerical simulations confirming the stability of the steady states of system (5.4). The basic reproduction number R 0 was approximated by using the time-averaged system.
In Figure 7, the approximated solution of system (5.4) approaches asymptotically to a periodic solution with the persistence of the disease once R 0 > 1. In Figure 8, we provide a magnified view of the limit cycle when R 0 > 1. In Figure 9, the approximated solution of system (5.4) approaches the disease-free periodic trajectory E 0 (t) = (S * (t), 0, 0, A * (t)) once R 0 ≤ 1.

Conclusions
When studying the CHIKV dynamics, it is important to consider that the contamination of uninfected cells can be realized via contact with CHIKV (CHIKV-to-cell transmission) and by contact with infected cells (cell-to-cell transmission). Moreover, disease spread can have seasonal peak periods and it is important to consider this when modelling the dynamics. In this paper, we have proposed an extension of the CHIKV epidemic model already considered in [5,21,22] by taking into account the seasonal environment. In the first step, we studied the case of an autonomous system where all parameters are supposed to be constants. We calculated the basic reproduction number and the steady states of the system. We gave the existence and stability conditions for these steady states. In the second step, we considered the non-autonomous system, gave some theoretical results and defined the basic reproduction number, R 0 through the use of an integral operator. We show that if R 0 ≤ 1, all trajectories converge to the disease-free periodic solution; however, the disease persists once R 0 is greater than 1. Finally, we gave some numerical examples that support the theoretical findings, including those for the autonomous system, the partially non-autonomous system and the fully non-autonomous system. It has been deduced that if the system is autonomous, the trajectories converge to one of the equilibriums of the system (2.1) according to Theorems 3 and 4. However, if at least one of the model parameters is periodic, the trajectories converge to a limit cycle according to Theorems 6 and 7.

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