2. The transmission rate is periodic to study impacts of the heterogeneity networks and seasonal diseases.
3. There exists a unique periodic solution which is globally asymptotically stable. 4 . The heterogeneity is important in accelerating disease spreading and increasing the amplitude of solution.
forcing; Basic reproduction number; Amplitude of the infected population
Seasonality, as a typical performance of periodicity, is an important aspect of epidemiology when studying vector-borne and air-borne diseases [1] . Outbreaks of seasonal diseases have been modeled and analyzed widely on the populations with homogeneous mixing in which the contact rate of all people in the population is the same, [2] [3] [4] [5] [6] . London et al. [2] demonstrated that sufficiently large seasonal variations in transmission generate periodic solutions. Aron et al. [3] showed that seasonal models of epidemics cannot account for small-amplitude recurrent epidemics of arbitrary periodicity. Bacaër et al. [7] defined the basic reproduction number of vector-borne diseases in the presence of seasonality. In 2007, Zhang et al. [8] obtained classical results for periodic epidemic models in patchy environment. Later in 2008, Wang et al. [5] provided an approach to calculate the basic reproduction number for such models. Also Liu et al. provided a threshold condition for seasonal models with a pulse or constant control schemes [9] . The impact of seasonal variation patterns on epidemic dynamics also was studied by [6] .
The homogeneous assumption used in the previous works to derive models and analyze the periodic behavior of seasonal epidemics is inadequate in real world as there is a large heterogeneity in contact rates among people in the population. Depending on the infection spreading process, many biological contagions require capturing a contact pattern in the model [10] . The contact patterns can be described by Network Theory [11] , where the nodes of the network represent people and the edges represent interaction/contact among people. The number of edges connected to nodes captures the heterogeneity of contact in population. A wide range of modeling epidemics on networks from Degree Based Mean Field models (DBMF) [12] [13] [14] , which is based on the assumption that all nodes with the same degree are statistically equivalent [13] , to pair-approximation models [15] , and then with more complex Agent based network models [16, 17] have been introduced to study epidemic spreading in heterogeneous population, with a focus on a single outbreak of epidemic. There is a few literature on studying the seasonality and periodic behavior of disease dynamics on networks with a little attention on network topology [18] [19] [20] . These studies used a discrete time Markov chain approach to implement recurrent infections such as seasonal influenza on networks.
Our goal of this study is to fill this gap by providing a simple network modeling approach, DBMF, for periodic phenomena and a complete stability analysis on the system. More specifically,
• We develop a non-autonomous heterogeneous DBMF SIS modeling approach, in which the time-dependent periodic transmission rate captures periodicity of the system.
• Then, we conduct stability analysis of the model to obtain a threshold condition of the existence and uniqueness of the periodic solution for the proposed model.
• Finally, in order to confirm our stability analysis result, we apply our framework and results to a Scale-free network as a case study for simulations. The selection of Scale-free network is because of its heavy-tailed statistical distributions and therefore, its better representation of real world systems [10] .
We consider an infection spreading on a general network with a degree distribution P k , k = 1, 2, · · · , n, and average degree k = n k=1 kP k . Denote N k = S k + I k as the number of individual with degree k, and S k , I k as the respected susceptible and infected individuals with degree k. Let P (k|k) be the probability that individual with degree k is connected to an individual of degreek, and Ik Nk be the probability that individual of degreek is infected, that is, we assume all infected individuals spread uniformly in the population. Then P (k|k) Ik Nk describes the probability that an individual of degree k is connected to an infected individual of degreek.
Assuming all the contacts are independent and summing up these contact probabilities over all possible degrees valuesk, the function Θ k = k P (k|k) Ik Nk defines the approximation of contact probability with a typical infected individual with an arbitrary degree.
The associated Degree-based mean-field SIS (DBMF-SIS) is described as follows [12] [13] [14] :
where β(t) is the transmissibility and β(t)kS k (t)Θ k (t) describes the number of individuals of degree k who become infected at time t, γ is the recovery rate for an infected person, and µ is death rate. Keeping birth rate rate equal to death rate, we assume total population is constant during course of epidemic. Also with the assumption that network is static, total population of people of degree k, N k , is constant. If we assume that network is uncorrelated [21] :
Nk . To take seasonality into account, we define the transmission rate β as a periodic function β(t) = β 0 (1 + δcos(2πt)), where δ is the strength of seasonality [3] . The general periodic transmission function β(t) can be represented as
where t 0 is a non-negative number specifying the initial periodic force, periodic transmission rate at time zero. Because the total population with degree k, N k is constant, we let s k = Sk Nk , and i k = Ik Nk as the fraction of susceptible and infected individuals with degree k. Thus, we can have the following rescaled model (2.2) from the original model (2.1):
where Θ = 1 k n k =1 k P k i k . For convenience, we will perform our further analytical and numerical studies on the proposed model (2.2).
First we prove the existence and positivity of the solution for the model (2.2).
is a solution of model (2.2) corresponding to initial condition x(0) = (s 1 (0), ..., s n (0), i 1 (0), ..., i n (0)) that satisfies:
Then the solution x(t) satisfies similar condition:
That is , the set Ω of all the solutions of model (2.2) defined as below is positively invariant
Proof. We show that if the initial value x(0) ∈ Ω, then for any t > 0, x(t) ∈ Ω. Summing up the two equations of model (2.2), we have
we only need to show s k (t), i k (t) ≥ 0 for any t > 0. First we show for all k i k (t) ≥ 0. Using proof by contradiction, assume there exist k ∈ {1, 2, · · · n} and t * , such that t * = inf {t > 0|i k (t) = 0, di k ds s=t < 0}. Therefore, for any 0 ≤ t < t * we have i k (t) ≥ 0. From the model we have
which shows that s k (t * ) < 0, therefore,
Therefore, for all 0 ≤ t ≤ t * , s k (t) ≥ 0, which is a contradiction indicating i k (t) ≥ 0 for any t > 0 and all ks. Because i k (t) ≥ 0 then we have
Thus, for any t ≥ 0, s k (t) ≥ 0.
Our model (2.2) is autonomous when δ = 0, because in that case we have β(t) = β 0 [1 + δsin(t 0 + ωπt)] = β 0 becomes a constant function of time. Therefore, by using next generation approach [22] , we can derive basic reproduction number in the absence of seasonality, shown by
. First, the model (2.2) always has a disease-free equilibrium E 0 = (1, 1, · · · , 1, 0, 0, · · · , 0).
By setting up the next generation procedure we have
and deriving their Jacobian matrices at E 0 as
where I is identity matrix, basic reproduction number R 0 δ=0 is computed as spectral radius of
However, our model (2.2) is non-autonomous, as β is a periodic function of time, therefore, matrix F is the functions of time, β(t) and F (t). In this case, general basic reproduction number shown by R 0 is spectral radius of a linear map L on the Banach space of ω−periodic functions defined in Equation (2.7) of [5] . We define the vectors F and V as F =
where,
Then we introduce the following linear ω-periodic equation:
and let W (t, s, z), t, s ∈ R t ≥ s, be the evolution operator of model (2.2) on R n . Then ∀t ≥ 0 we
. Let x(t) = (s 1 (t), ..., s n (t), i 1 (t), ..., i n (t)), we define assumptions (A1)-A(7) in [5] as follow:
x(t)) (i th element of vectors) are non-negative and continuous on R × R 2n + and continuously differential with respect to x(t). That is, each element of mentioned vectors is a directed non-negative transfer of individuals.
(A2) There is a real number ω > 0 such that for each
are ω-periodic in t. This condition describes a periodic environment.
.., 2n} , then V − = 0 for i = n + 1, ..., 2n. It means that if a compartment is empty, then there is no transfer of individuals out of the compartment.
(A4) F i = 0 for i < n, that is, the incidence of infection for uninfected compartments is zero.
.., 2n. Therefore, the population remains disease-free if it is disease-free at the beginning.
If we can show the conditions (A1)-(A7) holds then the reproduction numer R 0 would be the root of equation W (ω, 0, z) = 1. It is straightforward to see the conditions (A1)-(A5) of [5] hold.
For the condition (A6) we need to prove ρ(Φ(M (ω))) < 1 where Φ(M (ω)) is the monodromy matrix of the linear A-periodic system dX dt = M (t)X with:
Because matrix M is stable as −µ < 0 then ρ(Φ(M (ω))) = e −µω < 1. Therefore, the condition (A6) holds.
Similarly, the condition (A7) states that ρ(Φ(−V (ω))) < 1, where Φ(M (ω)) is the monodromy matrix of the linear ω-periodic system
Obviously we have ρ(Φ(−V (ω))) < 1, where Φ(M (ω)) is the matrix of the linear ω-periodic
Again, matrix −V is stable as −γ − µ < 0 then ρ(Φ(−V (ω))) = e −(γ+µ)ω < 1. Therefore, the condition (A7) holds.
Although having explicit formula for R 0 is not possible here, but using the theorems in [5] we can analyze the stability of equilibria conditioned to value of R 0 in the remaining part of this Section, and then in next Section, we verify our result by calculating R 0 numerically. First, we prove the global stability of disease-free equilibrium point E 0 for the case that R 0 < 1:
Proof. Based on Theorem 2.2 in [5] , if R 0 < 1, then E 0 is locally asymptotically stable and
Therefore, we only need to show that E 0 is globally attractive. By Lemma 3.1
we have s k (t) ≤ 1 for all k and t > 0. Therefore,
is globally asymptotically stable. By the comparison principle, in model (2.2), i k (t) → 0 and therefore s k (t) = 1 − i k (t) → 1 as t → ∞. It implies that E 0 is globally attractive, and therefore, globally asymptotically stable. Now we aim to prove the existence and uniqueness of the positive periodic solution of model (2.2) when R 0 > 1.
We define Ω 0 as interior set of Ω, that is, Ω 0 := {(s 1 (t), · · · , s n (t), i 1 (t), · · · , i n (t)) ∈ Ω| i k (t) > 0, k = 1, 2, · · · , n} and ∂Ω as its boundary set, that is, ∂Ω = Ω c 0 = Ω\Ω 0 . Denote that u(t, x 0 ) as the solution of model (2.2) with the initial values x 0 = (s 1 (0), · · · , s n (0), i 1 (0), · · · , i n (0)). Let
be the Poincare map associated with model (2.2). Following the fundamental existenceuniqueness theorem in [23] , u(t, x 0 ) is the unique solution of model (2.2) corresponding to initial condition u(0, x 0 ) = x 0 . From Lemma 3.1, Ω is positively invariant and P is point dissipative. Lemma 3.3. If R 0 > 1, and there exist δ > 0 such that for any initial condition x 0 = (s 1 (0), · · · , s n (0), i 1 (0), · · · , i n (0)) ∈ Ω 0 we have
Proof. Because R 0 > 1, then by Theorem 2.2 of [5] ρ(Φ F −V (ω)) > 1. We choose small enough value > 0 and define matrix M
for any δ > 0. Without loss of generality, assume that d(P m (x 0 ), E 0 ) < δ for all m ≥ 0. Because of continuity of periodic solution u for model (2.2), we have
For any t ≥ 0 we can define greatest integer m ≤ t ω as m := [ t ω ] such that t = mω +t for t ∈ [0, ω]. Therefore, we have
Assuming that u(t, P m (x 0 )) = (s 1 (t), ..., s n (t), i 1 (t), ..., i n (t)), this implies that 1− ≤ s k (t) ≤ 1 + , t ≥ 0, k = 1, 2, · · · , n. If x 0 − E 0 ≤ δ. Therefore,
Consider the linear system
By Lemma 2.1 of [8] , there exists a positive ω-periodic functionĝ(t) such that g(t) = e qtĝ (t) for
) > 1 and g(0) > 0, therefore, g(t) → ∞ as t → ∞. By comparison principle and for i k (0) > 0, (k = 1, 2, · · · , n), we then have i k (t) → ∞ as t → ∞, which contradicts the bounded condition in (3.5) . Therefore, the lemma is proved. Now we are ready to show the following theorem for the model (2.2) when R 0 > 1: Theorem 3.4. If R 0 > 1 then the model (2.2) has a unique positive periodic solution that is globally asymptomatically stable.
Proof. First we show that Ω 0 and ∂Ω are positively invariant.
If x 0 = (s 1 (0), · · · , s n (0), i 1 (0), · · · , i n (0)) ∈ Ω 0 , then by solving the model (2.2) for initial value
and
This proves that Ω 0 is positively invariant. We know that ∂Ω is a closed subset of Ω that includes P is uniformly persistent with respect to (Ω 0 , ∂Ω). By Theorem 1.3.6 in [24] , P has a fixed point (s * 1 (0), · · · , s * n (0), i * 1 (0), · · · , i * n (0)) ∈ Ω 0 . On the other hand, from model (2.2) and for any t ∈ [0, ω], we have
The periodicity of s * k (t) shows that for all t > 0, s * k (t) > 0. Following the processes as in Equation (3.8) and for all t > 0, we have i * k (t) > 0 for k = 1, ..., n. Therefore, E * = (s * 1 (t), · · · , s * n (t), i * 1 (t), · · · , i * n (t)) is a positive ω-periodic solution of model (2.2). To show that E * is the only solution -uniqueness of solution-we use the theorem 3.2 in [25] .
Because total population of individual of degree k is constant, we only consider the model for
and suppose Ω I is the projection of Ω over infected class, that is,
Let g = (g 1 , g 2 , · · · , g n ) : Ω I → R n be the right side of Eq.(3.9). We can easily see that, (i) g is
for all x ∈ Ω I . In addition by defining
for all x ∈ Ω I . Therefore, g is cooperative and Dg is irreducible. This proves that conditions (1) and (3) of theorem 3.2 in [25] hold. Now let ε ∈ (0, 1), and x 0 is in Ω I . For j = 1, 2, · · · , n,
This proves that g is strictly sub-linear, and thus, the condition (2) of theorem holds. Since admits a unique positive periodic solution E * which is globally asymptotically stable in Ω 0 .
In this Section we perform numerical simulations for the spreading of an infection via the model (2.2) on BA Scale-free networks. The degree distribution of our networks follows a power law equation P k = ck −3 , in which the factor c is a normalization factor, c = 1 n k=1 k −3 . All of the simulations are obtained with the model baseline parameters in Table ( 
In our first simulation and to validate the theorems 3.2 and 3.4, we illustrate the global stability of disease free equilibrium or the periodic solution for different status of R 0 and for different heterogeneity level of the network. We define total fraction of infected people at time t as i(t) = n k=1 i k (t)P k where i k (t) is the fraction of infected individuals with degree k and P k is the probability of hiving degree k. Fixing the parameters at their baseline values in Table (1) but selecting β 0 to keep R 0 ≤ 1, the steady state infection fraction i(t) is the disease free equilibrium, Figures 1(a) , 1(c), and 1(e). Increasing β 0 to have R 0 ≥ 1 causes that i(t) tends to endemic periodic solution as shown in Figures 1(b), 1(d) , and 1(f). These simulations are confirming the result of global stability of the disease free equilibrium for R 0 ≤ 1, and the periodic solution for R 0 ≥ 1, as stated in Section 4.
In our next inspection, we study the impact of β(t) = β 0 [1 + δsin(t 0 + ωπt)] on R 0 by varying the parameters β 0 and δ over their corresponding range, and for networks with different Table(1) and for n = 100, the growth of R 0 with β 0 is around four times faster than that of case n = 1.
It was clear to expect the linear correlation between β 0 and R 0 as R 0 is a coefficient of β 0 , where this coefficient is some function of other infection parameters and network structure [13] . 
We end up this Section by studying the impact of transmission rate -or specifically β 0 and δ-on the amplitude of infection fraction when it reaches its steady state, i * . The Figure ( and considering the wider range for β 0 , we observe a non-monotonic behavior for amplitude of infection fraction, Figure (4) . For small values of β 0 , i * increases to its maximum point, which this increment is bigger for more heterogeneous networks (network with bigger n). But when β 0 increases more i * start decreasing slowly, independent of network structure.
We developed a degree-based mean field SIS model with a seasonal transmission rate to understand the impact of periodic characteristics on infection thresholds and dynamics during epidemics. We analyzed the basic reproduction number for the model and used it to obtain the asymptotical and global stability conditions of the equilibrium points. Our theoretical results in Section 4 were investigated and verified via numerical simulations in Section 5 on Scale-free networks as highly heterogeneous network.
By analyzing the basic reproduction number R 0 , we obtained the asymptotical stability condition of the equilibria: if R 0 ≤ 1, then the disease-free equilibrium point is globally asymptoti- cally stable, otherwise existence of a unique, periodic and globally stable solution of the system is guaranteed. Conducting numerical simulation to illustrate the theoretical result, we observed that when R 0 ≤ 1 then the solution stabilizes at disease-free equilibrium point, while for R 0 ≥ 1 the non-zero solution for the system is periodic.
To understand how seasonal forcing affects the growth of infection spreading, we plotted basic reproduction number R 0 as a function of non-periodic transmission rate β 0 and strength of seasonality δ. Obviously R 0 increases as non-periodic transmission rate β 0 increases, and this increment is linear which its slope increases as the heterogeneity of network increases. On the other hand, basic reproduction number R 0 has a negative and nonlinear correlation with strength of seasonality δ, that is, R 0 decreases as δ increases. This reveals the fact that infection will spread slower for a shorter period of infection spread, and similar to β 0 this variation is more intense for more heterogeneous network, Figure ( 2).
Beside the basic reproduction number R 0 , the amplitude of periodic solution at steady state is affected by transmission rate factors β 0 and δ. Through simulation, we found out that this amplitude linearly increases by increasing the periodic strength of seasonality parameter, δ. The level of this correlation increases as the network heterogeneity increases, which is similar to results in Barthélemy et al. [26] . On the other hand, the impact of β 0 on the amplitude for small values of this parameter is almost negligible, but for bigger value of β 0 the amplitude of infection start to increase almost linearly. Similar to results for δ, more heterogeneous network amplifies the correlation between β 0 and amplitude of periodic solution, Figure ( 3).
Although the basic reproduction number R 0 is a monotonic function of seasonal force, the periodic solution shows a different behavior. More precisely, the amplitude of periodic solution at steady state is first increasing function of β 0 up to its maximum point, but then start decreasing when β 0 becomes bigger. This non-monotonic trend is independent of the value of strength of seasonality or the heterogeneity level of network, Figure (4) .
Periodic solution generated by sufficiently large seasonal forcing is result of coupling effect of β 0 and δ [2] . Since, control of periodic infection is dependent on amplitude of the endemic, it is necessary to control the strength of seasonality δ as well as transmission rate in the absence of periodic solution β 0 . In addition, the network heterogeneity can impact the prevalence of infection and therefore, needs to be controlled in the course of outbreak.
When the network is static, that is, it is not changing over time, the degree-based mean field network model provides an strong approximation of infection dynamics [10] . But intrinsically, periodic infections such as seasonal diseases evolve at a very slow time scale and even may slower than the time scale of topological changes in network, and therefore, it may not be suitable to have static networks. In this limit, if we assume the network is constantly rewired while preserving its degree and joint degree distribution-annealed network approximation [27, 28] -then we can extend the degree-based mean field network model to be a more suitable approximation of the real world situation.
Though, the proposed degree-based network model provides a natural first step of realistic network models, it cannot account for behavior change such as social distancing and or isolation which destroy degree distribution. When nodes (individuals) become more connected the transmission of infection speeds up and causes an outbreak. Therefore, via behavior change the "activity" of the network should be appropriately reduced to reduce the harm of the epidemic [29] . The proper example is the novel coronavirus pneumonia outbreak in China, in which many towns and communities are recently being isolated. These measures to close the district effectively reduce the scale of the outbreak [30] . These scenarios require richer classes of model structure and hierarchies of model properties such as degree-corrected block model [31] or
Agent-based network model.