On the early epidemic dynamics for pairwise models.

The relationship between the basic reproduction number R 0 and the exponential growth rate, speciﬁc to pair approximation models, is derived for the SIS, SIR and SEIR deterministic models without demography. These models are extended by including a random rewiring of susceptible individuals from infectious (and exposed) neighbours. The derived relationship between the intrinsic growth rate and R 0 appears as formally consistent with those derived from homogeneous mixing models, enabling us to measure the transmission potential using the early growth rate of cases. On the other hand, the algebraic expression of R 0 for the SEIR pairwise model shows that its value is aﬀected by the average duration of the latent period, in contrast to what happens for the homogeneous mixing SEIR model. Numerical simulations on complex contact networks are performed to check the analytical assumptions and predictions.


Introduction
The success of invasion processes governed by short-range interactions crucially depends on the early development of local densities around the invading individuals, which constitute the so-called invading clusters [3]. These local spatial patterns of individual states define the environmental conditions experienced by the invaders, so determining their final fate. Invading clusters have been observed in lattice models in Ecology [10,14,16] as well as in simple epidemic network models [3,17,18] where invaders correspond to infectious individuals introduced into a totally susceptible population. In both settings, pairwise models offer a useful approach to the study of their dynamics. In such models, local densities are described by either the conditional probabilities that an invader has a neighbour in a particular state [14,16], or the mean fraction of neighbours of a given type around an infectious individual [3,18,40], or by the average number of neighbours of a given type per infectious individual [17]. In [3,18,40] correlations among infection states are defined in terms of local densities and are used to describe the early development of spatial patterns.
The early dynamics of local densities around infectious individuals have been studied from their limit equations [3,17,18] for the susceptible-infected-susceptible (SIS) and susceptible-infectious-recovered (SIR) models. These equations are derived by assuming that the expected number of susceptible individuals, [S], and of susceptible-susceptible links, [SS], are approximated to the total number of individuals (nodes), N , and to the total number of links (non-ordered pairs), L, respectively. Under these assumptions, the expected number of susceptible-infectious pairs (infectious links), [ i.e., the mean number of susceptible, recovered, and infectious individuals around an infectious individual, respectively. Conversely, the so-called mean-field models, like the deterministic SIS, SIR, and SEIR (E for Exposed) models or unstructured competition models in Ecology, cannot capture such an important feature of the early dynamics of a spatially extended system. All these models assume fully homogeneous mixing of individuals or, when dealing with populations of sessile organisms, long-range interactions. So, their predictions about the early dynamics are quite different from those obtained from pairwise models.
The aim of this paper is to analyse the early dynamics of an epidemic by using a pair-approximation model. We will compute the expression of the basic reproduction number, R 0 , from the equilibrium of the limit equation for the local density [SI]/[I] at the early phase of an epidemic outbreak. Moreover, we will obtain the relationship between R 0 so derived and the largest eigenvalue, λ 1 , of the Jacobian matrix of the full pairwise model around the disease-free equilibrium at the beginning of an epidemic. As we will see, when it is positive, λ 1 corresponds to the exponential growth rate of the cumulative number of cases during the early stages of an epidemic outbreak, once local correlations around initially infected individuals (primary cases) have been developed. This interpretation of the early dynamics given by linearisation of deterministic pairwise models is in agreement with previous descriptions given in terms of the development of invading epidemic clusters around primary cases [3,18]. From the point of view of public health management, it could help to obtain more accurate early estimates of R 0 , which are needed to evaluate the transmission potential of a new disease and assess the impact of different intervention strategies [5]. With this respect, it is interesting to note that the same relationships between R 0 and the initial growth rate, r, of an epidemic have been derived from homogeneous mixing models, and have been used to compute R 0 by estimating first r for the cumulative number of cases during the exponential-growth phase of an epidemic outbreak [5,22,29,44,45]. Even though, for populations with heterogeneous mixing, the description of the initial epidemic phase as well as the estimate of r offered by these models can greatly differ from the ones given by pairwise and network models [27].
The SIS, SIR, and SEIR deterministic models assume that the infectious and latent periods are exponentially distributed and populations are homogeneously mixed. For these models, the algebraic expression of R 0 is the same, namely, βN/γ with β being the transmission rate across an infectious contact, N is the population size, and γ the recovery rate [6]. The same expression of R 0 with a general mean duration of the infectious period instead of 1/γ follows for the age-of-infection models, in which infected individuals are classified according to the time since infection (age of infection) [46]. On the other hand, the initial exponential growth rate r of an epidemic varies from one model to another and, consequently, so does the relationship between R 0 and r [33,44,45]. In particular, it is known that the latent period affects the growth rate of the epidemics but, by contrast, it plays no role in R 0 [6], even in age-of-infection models [45,46]. Precisely, for the SIS and SIR deterministic models, R 0 = 1 + r/γ with r = β − γ, whereas, for the SEIR deterministic model, R 0 = (1 + r/θ)(1 + r/γ) where r now also depends on θ, the rate of leaving the exposed class [5,22]. The same quadratic relationship between R 0 and r has been obtained from the so-called Euler-Lotka equation for the SEIR model in [33,44]. Here we prove that, although the same linear and quadratic relationships hold if one formally replaces r by λ 1 , for pairwise models the algebraic expression of R 0 differs from one model to another. The latter implies that different R 0 estimates will be obtained by using methods based on algebraic expressions of R 0 [5].
This study is carried out for an extension of the SIS, SIR, and SEIR pairwise models that includes random rewiring of susceptible individuals. It is assumed that each susceptible individual breaks off a connection to one of its infectious/exposed neighbours at given rates, and reconnects to a randomly chosen susceptible or recovered individual. An SIS pairwise model with such a rewiring was introduced for the first time by Gross et al. in [13], as a way to include behavioural responses to epidemics, a key ingredient to understand the effects of individual behaviour on epidemic dynamics [9]. A detailed analysis of the impact of rewiring in the early dynamics of such an SIS pairwise model under different types of network topologies has been recently done in [17]. Other reconnection mechanisms have been recently proposed in the literature (see, for instance, [20,32,34,35,41,42,43,47]).

The SIS-ω and SIR-ω pairwise models
To simplify the presentation of the results, we start this section considering a unique system of equations that encompasses both the SIS and the SIR pairwise models with rewiring, henceforth referred as SIS-SIR-ω model. Later on, when dealing with the early dynamics, we will treat each of them separately.
Let β be the transmission rate across an infectious contact, and let γ and µ be the rates of recovery with and without immnunity of an infectious individual, respectively. Finally, let ω be the rewiring rate of susceptible individuals as described in the Introduction.
Using expected numbers as state variables and according to the previously introduced notation, let [S], [I], and [R] be the expected number of susceptible, infectious, and recovered individuals, respectively, which satisfies [S] + [I] + [R] = N . Analogously, [ij] stands for the expected number of (non-ordered) pairs of connected individuals in states i and j, respectively, and [ij l] for the expected number of (non-oriented) triples whose individuals are in states i, j, and l (i, j, l ∈ {S, I, R}). Using these definitions, the differential equation for the number of infectious individuals can be written as where [SI]/[I] is the mean number of susceptible neighbours per infectious individual, and taking µ = 0 for the SIR model and γ = 0 for the SIS model. Now, let us compute the mean number of neighbours in state i (i-neighours) around a j-individual that already has at least one l-neighbour, that is, that belongs to a (j, l)pair, averaging over all the (j, l)-pairs. Denote this number by Q(i|jl). If l ̸ = i, each triple of type (i, j, l) contributes with one i-neighbour of a j-individual that already has an l-neighbour. Then, Q(i|jl) = [ijl]/[jl]. When l = i, then each triple of type (i, j, i) contains two (i, j)-pairs and it is counted twice when we average over all the (i, j)-pairs. Therefore, the mean number of i-neighbours of j-individuals in (i, j)-pairs is where it is used that a j-individual that belongs to an (i, j)-pair has one i-neighbour for certain, in addition to the mean contribution to i-neighbours from (i, j, i)-triples per (i, j)-pair. The factor 2 comes from the fact that the triples are not oriented (cf. [3,39] for the oriented case).
Taking these averages into account, the system of differential equations governing the dynamics of the SIS-SIR-ω model reads in the third and fifth equations comes from the expression of Q(i|ji), and takes into account that the total number of non-ordered (S, I)-pairs decreases by two when the infection reaches the central S of an (I, S, I)-triple as well as that the number of (I, I)-pairs increases by two. Remarkably, this factor is missing in the previous formulations of the SIS-ω model [7,12,13,17,32,47].
For any rewiring rate w ≥ 0, the mean degree of the network, k = 2L/N , is constant over time because rewiring keeps constant the total number of links L. However, the nodal degree distribution, p k , does change with time due to the rewiring of susceptible nodes. This means that the mean degree of susceptible (k S ), infectious (k I ), and recovered nodes (k R ) will also change over time. Therefore, the number of equations in system (3) cannot be reduced by using constraints on the total number of links leaving nodes that are in a given state like, for instance, 2[SS] + [SI] + [SR] = k S [S], because k S = k S (t) is unknown.
As usual in pair approximation, we close this system by assuming the statistical independence at the level of pairs. This implies that the expected number of (i, j, l)triples is approximately equal to the number of neighbours of those j-nodes in (i, j)-pairs, (q − 1) [ij], times the conditional probability that a j-node has a neighbour in state l. In uncorrelated heterogeneous networks and in absence of clustering, this amounts to the following closure for the expected number of non-oriented (i, j, l)-triples: where q = k 2 /k is the expected degree of an individual reached by following a randomly chosen link. Similarly to what we observed in (2), the factor 1/2 in the third expression of (4) appears because, when we count (non-oriented) triples that run through a given pair, those triples with endpoints having the same state are counted twice. In particular, This triple closure in terms of q is the natural generalization of the one used in regular networks, where all the nodes have the same degree, to networks with heterogeneous degree distributions. Since the probability that a randomly chosen link points to a node with degree k is proportional to k, the expected degree of a neighbour in uncorrelated networks is q = ∑ k kq k with q k = kp k /k [6]. Note that q = k for regular networks, which leads to the triple closure usually assumed in pair approximation [19,30]. For uncorrelated networks with Poisson degree distribution, q = k + 1 and, hence, (q − 1)/k = 1, which is the value used in the closure (4) in [13]. For uncorrelated networks with exponential degree distribution with k ≥ 0, q = 2k which amounts to (q − 1)/k = 2 − 1/k. A further refinement of the closure (4) should take into account that the degree distributions p S k , p I k , and p R k of susceptible, infectious, and removed nodes, respectively, change over time [2] and noticeable differences between these distributions arise when the variance of p k is high enough. For instance, it is known that, even for networks with a low variance in their initial degree distribution, p I k and p S k become broadened at the endemic equilibrium with the mean degree of susceptible nodes, k S := ∑ k kp S k , being greater than k, whereas the mean degree of infectious nodes, k I := ∑ k kp I k , is less than k [13]. Therefore, the precise form of the triple closure will depend on the state of the central node because it determines the degree distribution that must be used to compute q and k. However, we will assume that the degree distribution has been changed by neither the epidemic progress nor rewiring because we will restrict our analysis to the initial stages of the epidemic dynamics.
On the other hand, if the initial network configuration has no clustering, as it is the case for networks generated with the configuration model, the random disconnection of susceptible nodes from infectious ones and their uniform rewiring to susceptible (and removed) nodes are not expected to generate clustering in the network architecture. So, the correction term accounting for clustering will not be taken into account [15,18,19,39].
Upon introducing the triple closure (4) into system (3) and neglecting the last two equations, the SIS-SIR-ω model becomes where z := (q − 1)/k. Note that these equations depend on neither [SR] nor [IR] which allowed us to neglect the last two equations in (3). Moreover, when γ = 0 (SIS-ω model) the first and last equations in (5) are redundant and this system is reduced to only three equations.
Remark : Note that, when (2) and (4) 3 are introduced into the equations, the factors 2 and 1/2 compensate each other and, so, they do not appear in the final formulation of the pairwise model. This explains why previous works dealing with the SIS-ω model and using non-oriented triples, arrive at the equations given by (5) with γ = 0, although both factors are missing in their derivation [7,12,13,17]. depends on q, that is, it contains the information about the expected degree of the initially infectious individuals, called secondary cases. For this reason, we talk about average infectious individuals or average secondary cases. For highly heterogeneous degree distributions and low rewiring rates, these individuals will have degrees higher than the mean degree k in the population. Moreover, it is known that early spatial correlations between the infection states of neighbouring individuals (invading clusters) are formed within a couple of generations [8,17].

The early epidemic dynamics
This expected number of newly produced infections (secondary cases) generated by a single infectious individual once local correlations are established is, indeed, the definition of the basic reproductive number R 0 given for network models [8,18]. We adopt the average duration of an infectious period, (µ + γ) −1 , as a measure of one generation time [38] In this sense, note that very high rewiring rates will produce a lot of isolated infectious individuals which will lead to ( [I] Just after the introduction of the first infectious individuals, one has → 0, and the corresponding limit system is d dt where q = k 2 /k is computed from the degree distribution at time t = 0.

Early dynamics of the SIR-ω model: µ = 0
Introducing µ = 0 into system (7) which has a unique strictly positive equilibrium In this case, x * is globally asymptotically stable and, substituting it into Eq. (6), we have For regular networks (q = k 0 ) and no rewiring (ω = 0), this expression of the basic reproductive number reduces to the one given in [18].
From the conditions x * > 0 and R 0 < 1, it follows that the relationship ω < β(q −2) < γ + ω guarantees both the epidemic extinction and a positive mean number of susceptible individuals surrounding the initially infectious individuals as the epidemic dies out (i.e., 0 < R 0 < 1). Following [17], we call this situation "transmission-dominated scenario" of the epidemic extinction. Conversely, for high-enough rewiring rates, x * = 0 is the only nonnegative equilibrium. This means that the number of susceptible neighbours goes to 0 faster than the total number of infectious individuals in the population. This situation corresponds to the so-called "rewiring-dominated scenario" of the extinction and implies R 0 = 0.

Early dynamics of the SIS-ω model: γ = 0
For the sake of completeness, in this section we present some results for SIS-ω model which have been obtained in [17]. In this model, both equations in (8) are coupled and two phase portraits are possible depending on whether the nullclines intersect to each other in the first quadrant or not. Let us denote the solution to system (8) by (x, y). If β q ≤ ω, there is no such intersection and the only nonnegative equilibrium is the origin, which is globally asymptotically stable. When β q > ω, (0, 0) is unstable and a positive equilibrium (x * , y * ) appears at the intersection point of the nullclines. From the phase portrait if follows that this equilibrium is globally asymptotically stable [17]. In particular, from the expression of the nullcline y ′ = 0, we obtain y * < 2. Moreover, introducing the expression of y * obtained from the nullcline y ′ = 0 into the equation of the nullcline x ′ = 0, if follows that x * is the positive solution of the equation (1 + βx/µ)(ω + βx − β(q − 2)) = 2β. Note that the strict positivity of both terms in the left-hand side (lhs) implies that x * > q − 2 − ω/β. Therefore, the value of ([SI]/[I]) * for the SIR-ω model is always less than the one for the SIS-ω model, which implies a lower value of R 0 .
For β q > ω, introducing the expression of the unique positive solution x * of previous equation into Eq. (6) with γ = 0, we obtain the basic reproduction number for the SIS-ω model, namely, Note that the condition R 0 ≥ 1 for an initial spread of an epidemic can be rewritten as which, for ω = 0 and q = k 0 , corresponds to the well-known expression of the epidemic threshold in lattices and regular random networks.

R 0 and the local stability analysis of an epidemic outbreak
Once the description of the early dynamics of the epidemics is given in terms of the mean number of susceptible neighbours per infectious individual, we want to relate these results to the behaviour of the absolute numbers around the disease-free equilibrium (DFE) of system (5) The eigenvalues of this matrix are real. One of them is −(µ + γ) < 0. The other two, λ 1 and λ 2 , are the eigenvalues of J 0 DF , the 2 × 2 submatrix obtained by deleting the first column and the first row of J DF . They satisfy λ 2 < 0 and λ 1 > λ 2 . So, the largest eigenvalue of J DF will be greater than or equal to 0 if and only if λ 1 ≥ 0, that is As expected, the conditions for the epidemic spread obtained from (9) and (11) for both the SIR-ω and the SIS-ω models are special cases of (13). Therefore, for these models the linear approximation of system (5) leads to the same epidemic threshold that the limit system (8) for the local densities. However, one thing that is missing is the precise relationship between R 0 obtained from the limit system (8) and λ 1 .
Instead of working with the cumbersome expressions of R 0 and λ 1 , an easy way to find this relationship is to compute the ratio of [SI] to [I] during the initial exponential growth of the epidemic spread given by system (5) and, then, see if this ratio reaches a stationary value in this early stage. From the linear approximation of system (5)  [SI](t) = K 1 e λ 1 t + K 2 e λ 2 t with K i being constants depending on the initial condition and on the eigenvectors associated with λ 1 and λ 2 . Upon introducing this expression for [SI] into Eq. (5) 2 , we obtain a non-homogenous linear differential equation whose general solution is with K 3 being an integration constant. Finally, from these expressions we get [I] as t → ∞. This convergence is exponentially fast for all λ 1 ̸ = −(µ + γ), or proportional to t −1 otherwise. From this limit and the definition of R 0 given by (6), we arrive at the expression which can be algebraically checked for λ 1 > −(µ+γ) with µ = 0 (SIR-ω model) and γ = 0 (SIS-ω model) using (9) and (10), respectively. This relationship implies that the threshold condition R 0 = 0 for a rewiring-dominated scenario is equivalent to λ 1 = −(µ + γ).
This simple result tell us how we must understand the linearisation of deterministic pairwise epidemic models around the DFE. In particular, it says that the initial growth rate of the epidemic spread given by the largest eigenvalue of the Jacobian matrix J DF measures the growth of the infected population once the local densities around infectious individuals are established.

SEIR-(ω 1 , ω 2 ) model
Let us now consider a situation where the dynamics at the early stage of the epidemic and the relationship between R 0 and the largest eigenvalue of the linearization around the DFE are a little bit more complex. The reason is the existence of another class of individuals, the Exposed (E) ones, who are infected but not yet able to transmit the infection. After a latent period, exposed individuals become infectious at a rate θ. The relationship between the average duration of the latent period, 1/θ, and that of the infectious period, 1/γ, will define two different initial scenarios for the evolution of [SI]/[I].
Assuming that some symptoms of the disease can be present in exposed individuals, we allow S-individuals to disconnect from one of their E-neighbours and reconnect to a randomly chosen S or R at a rewiring rate ω 2 ≥ 0 (ω 2 = 0 if exposed individuals are asymptomatic). As with the SIR-ω model, S-individuals also disconnect from one of their infected neighbours and reconnect to an S or R at a rate ω 1 . According to these rewiring processes, once the triple closure (4) is introduced into the equations of the SEIR-ω model,

The early epidemic dynamics
As in the analysis of the SIS-SIR-ω model, we wish to obtain the initial dynamics of the local density [I] Therefore, there are four possible equilibrium values ζ * , namely, ζ * 1 = 0, ζ * 2 = γ/θ − 1, and the solutions ζ * 3 and ζ * 4 to βθ(q − 1) − (θ + ω 2 + θζ * − γ)(β + ω 1 + θζ * ) = 0, and the equilibria of system (16) are Since we are only interested in nonnegative equilibria, two classes of phase portraits are possible. For γ < θ, we have ζ * 2 < 0 and P 1 is the only equilibrium on the boundary of R 3 + . For γ > θ, there are two equilibria on the boundary of R 3 + , P 1 and P 2 . To find the conditions for the existence of positive P 3 and P 4 , it will be convenient to compute the Jacobian matrices of system (16) around P 1 and P 2 . If J i denotes the Jacobian matrix at P i , then Using the expression of J 1 , it is easy to see that ζ * 3 and ζ * 4 correspond to the eigenvalues of J 1 that are not λ 3 = γ − θ divided by θ, that is, where tr(J 0 i ) and |J 0 i | denote, respectively, the trace and the determinant of the 2 × 2 submatrix obtained by the intersection of the two first rows and the two first columns of J i . From these expressions and those of J 0 i , it follows that: • If γ < θ, then tr(J 0 1 ) < 0. Hence ζ * 4 < 0 and ζ * 3 > 0 if |J 0 1 | < 0. This means that P 3 is the only positive interior equilibrium and that it bifurcates from P 1 when the latter becomes unstable (|J 0 1 | = 0). • If γ > θ, P 1 is unstable (λ 3 > 0) and we need ζ * > γ/θ − 1 in order to have x * , y * > 0. However, it is easy to see that ζ * 4 < γ/θ −1 and, moreover, that ζ * 3 > γ/θ −1 ⇐⇒ |J 0 2 | < 0. This implies that P 3 is the only positive interior equilibrium and that it bifurcates from P 2 when the latter becomes unstable (|J 0 2 | = 0). Note that having P 3 > 0 means ([SI]/[I]) * > 0 and, hence, R 0 > 0. Precisely, in this case, from Eq. (6) we have ) .
So, there are two possible ways for a transition between a transmission-dominated scenario and a rewiring-dominated scenario depending on the equilibrium from which P 3 bifurcates: P 1 if γ < θ, and P 2 if γ > θ.

R 0 and the local stability analysis of an epidemic outbreak
Proceeding along the same lines as in the SIR-SIS-ω model, now we want to derive a relationship between R 0 obtained from the limit system (16) and the linear stability analysis of the full system (15).
Denoting the largest eigenvalue of J 0 by λ 1 and the other one by λ 2 , it follows, first, that λ 2 < 0 and, second, that λ 1 > 0 if and only if This condition shows that if ω 2 = 0 (no rewiring of susceptible individuals from their exposed neighbours) or θ → ∞, i.e., when the average duration of the latent period tends to 0, we recover the condition for the instability of the DFE in the SIR-ω model. In general, taking β as a tuning parameter, this condition says that the epidemic threshold given by the value of β such that λ 1 (β) = 0 is higher in the SEIR-(ω 1 , ω 2 ) model than in the SIR-ω model.
From the Jacobian matrix of system (15) around the DFE, it follows that the growth of [SI] during the initial exponential phase is only determined by the eigenvalues of J 0 , that is, [SI](t) = K 1 e λ 1 t + K 2 e λ 2 t . For the sake of simplicity, we shall assume that λ 1 ̸ = {−θ, −γ} and θ ̸ = γ. Introducing [SI](t) into the second equation of (15), this becomes a non-homogeneous linear differential equation for [E](t) whose general solution is Then, upon introducing this expression into the third equation of (15) and solving the resulting non-homogeneous differential equation for [I](t), we have where K 3 and K 4 are integration constants.

Now, using the previous expressions of [SI](t) and [I](t)
, the mean number of susceptible individuals surrounding an infected one at the early stage of the epidemic is given by: as t → ∞. Note that this convergence is exponential because we are assuming λ 1 ̸ = {−θ, −γ}. For the limit case λ 1 = max{−θ, −γ}, one can see that the convergence to 0 is proportional to t −1 if θ ̸ = γ, or proportional to t −2 if θ = γ. Finally, recalling the definition of R 0 given by (6), we obtain As with the SIR-ω model, this equality can be algebraically checked for λ 1 > max{−θ, −γ} using the expression (18) obtained from the equilibrium of the limit system for [SI]/[I].
From the expression of λ 1 (see Appendix) and taking now θ as a tuning parameter, it follows that λ 1 (θ) → β(q − 2) − γ − ω 1 =: λ SIR 1 , the largest eigenvalue of the matrix J 0 DF obtained from (12), as θ → ∞. Moreover, since λ 1 (θ)/θ → 0 as θ → ∞, we also have that R 0 → R SIR 0 := 1 + λ SIR 1 /γ, the value of R 0 for the SIR-ω model, when θ → ∞ (see Appendix for details). Interestingly, (23) remains true even without rewiring (ω 1 = ω 2 = 0), which implies that the difference in the value of R 0 between the SIR and SEIR pairwise models arises because of the time evolution of the local densities around infected individuals. In this special case of no rewiring, however, the epidemic threshold R 0 = 1 is the same for both models.

Simulations
To show the accuracy of the analytical predictions about the initial spread of an SEIR epidemic and how they depend on the network architecture, we performed stochastic simulations on networks of size N = 10000 with Poisson, exponential, and power-law degree distributions, the latter with exponent 3 (p k ∼ k −3 ). Given a mean degree k, Poisson networks were generated by connecting, with probability p = k/N , every pair of nodes in the network. Exponential and power-law networks were generated from random sequences of N degrees drawn from the corresponding probability distributions. Each node was given a degree k i from the sequence and, then, was randomly connected to k i nodes (configuration model). Power-law degree sequences where generated having minimum degree k 0 = 5 and maximum degree (cut-off) k c = k 0 N 1/2 = 500, which avoids the presence of degree-degree correlations [4] and allows to compute the value of q for this kind of degree distribution. For these values of k 0 and k c , the mean degree is 9.86. The mean degree of Poisson and exponential distributions is 10.
On these networks, the stochastic time evolution of the infection spread is simulated by means of the Gillespie algorithm [11]. The initial number of infected nodes is 10 (0.1% of the population size) which are drawn from the whole set of nodes with the same probability p = 1/N . For each network and parameters combination, we run 100 simulations with different initial sets of infected nodes in order to average the outputs.
The main output of the simulations will be the time evolution of the effective reproductive number R t := β/γ · [SI](t)/[I](t), which gives the mean number of infections produced by a single infectious individual during the course of an epidemic. Comparing both panels in Figure 2, one can clearly see the existence of a stationary value of R t , predicted from the limit equation for the dynamics of [SI]/[I], as long as the growth of [E] and [I] is in its initial exponential phase and the depletion of susceptible individuals is negligible. In particular, the right panel depicts how R t fluctuates around R 0 , the value of R t for [SI]/[I] given by (17) 1 , the first component of the equilibrium of the limit system (16). This figure also shows that, as expected, I(t) attains its maximum, the epidemic peak, when R t = 1. From this moment onwards, each infectious individual produces, on average, less than one infection (R t < 1) and the epidemic declines.
When positive rewiring rates are considered, the number of susceptible neighbours per infected individual is necessarily lower than without rewiring. This implies that R t will fluctuate at lower values and, if rewiring is large enough, its expected initial increase will not happen. Figure 3 illustrates how, for the same rewiring rates, the initial behaviour of R t for a Poisson network differs from the one in the other two types of networks. This figure also shows how the amplitude of the fluctuations increases with the variance of the degree distribution, being the prediction for the Poisson network the most accurate (a similar accuracy is observed for other values of the parameters). Finally, differences in the early time evolution of R t are clearly noticeable when simulations are performed on different networks generated according to the same degree distribution, except in networks generated from a Poisson degree distribution due to its low variance.
Stochastic simulations on heterogeneous networks show an initial increase of R t from the value kβ/γ for the randomly-chosen primary cases to values close to R 0 predicted from the expression of [SI]/[I] given by (17) 1 , followed by a decrease of R t (middle and bottom panels in Figure 3). This behaviour reflects the impact of highly connected individuals during the initial exponential-growth phase, and agrees with what was actually seen in the early growth of SARS in Hong Kong in 2003, where two clusters of cases generated by "super-spread" events were reported [31]. It also agrees with the transmission of influenza A (H1N1) in Japan, with a first phase corresponding to the outbreak of May 2009, where infections were mainly confined to school settings, i.e., to closed clusters with very high average contact rates and high values of R t . The second phase occurred from May 29 to July 14. It consisted of the secondary cases generated by school clusters which lead to lower estimates of R t [28].

Discussion
Using the pair approximation with the triple closure introduced in [17], we have analysed and compared the early dynamics of the SIS, SIR, and SEIR epidemic models defined on adaptive networks where connections between infected, exposed and susceptible individuals are broken off by the latter through a rewiring mechanism. Following [13,17], for the SIS and SIR models it is assumed that each susceptible individual that looses a connection with an infectious neighbour establishes a new one with another susceptible or recovered individual chosen at random. For the SEIR model, it is allowed that susceptible individuals also break off links with exposed neighbours. The used triple closure allows to deal with networks with different architectures and to study the corresponding initial behaviour of the epidemic spreading. In particular, Poisson, exponential and power-law degree distributions have been considered for numerical simulations. This choice of network topologies accounts for the variability of contact patterns found in several social networks, which ranges from the homogeneous mixing in Poisson networks to the high contact heterogeneity in networks with power-law degree distributions [2].
For these pairwise models, R 0 has been computed from the equations for the early dynamics of the average number of susceptible neighbours of an infectious individual, [SI]/[I] (see Eqs. (9), (10), (18)). Its expression includes components of the natural history of the disease (mean duration of the latent and infectious periods, transmission and recovery rates) as well as summary statistics of the contact network and the rewiring rates. Even without rewiring, the expressions of R 0 for the SIR and SEIR models are different from each other, in contrast to what happens to their homogeneous mixing versions for which R 0 = βN/γ [6,23,46]. Only when rewiring does not affect exposed individuals, the epidemic threshold given by the equation R 0 = 1 coincides in both models. However, even in this case (ω 2 = 0), the values of R 0 are different to each other under any other choice of the parameters values.
The importance of contact heterogeneity has been long recognized in the epidemiology of sexually transmitted infections [1], for which contact tracing is a useful control measure [25]. Information about contact patterns is also relevant for the control of infectious diseases passed by casual social contact, as it was shown for respiratory-spread infectious diseases like SARS [27] and influenza A (H1N1) 2009 in [28]. For these illnesses, initial estimates of R 0 were based on models that assume fully mixing and, some of them, on transmission data from crowded settings, like hospitals, schools, and apartment buildings, with unusually high contact rates. Consequently, extrapolating these R 0 estimates to the population at large overestimated the infectiousness of the disease because common contact rates are much lower outside these closed settings [28]. Therefore, caution must be exercised in the use of very early parameter estimates due to their inherent high uncertainty. For the 1918 influenza pandemic, for example, a significant reduction in the width of the confidence interval for R 0 was achieved in [5] by increasing the length of the fit during the exponential-growth phase. Other estimates can be obtained by fitting more sophisticated models as, for instance, the sort of pairwise models we have considered or structured models based on the so-called next generation matrix [6,28]. In this sense, although the population contact structure is unknown for commonly gathered epidemiological data, the parametrization of models involving this kind of information can be feasible by collecting the required data by means of survey-based or device-based techniques. An alternative consists of inferring the contact heterogeneity at the population level from commonly gathered data. This approach has been recently considered in [37] for different types of contact networks.
The second contribution of the paper is the relationship between R 0 and the largest (non-zero) eigenvalue, λ 1 , of the Jacobian matrix of the pairwise model around the diseasefree equilibrium with the proviso that λ 1 > min{−θ, −γ} (see Eqs. (14) and (23)). In [8], working with a generalization of pairwise models, R 0 is expressed in terms of the dominant eigenvalue of a matrix that contains all of the information about the types of partnerships present in the network (see [15,21,26] for similar models). In that work, Eames and Keeling comment that "In practice R 0 is calculated from the initial growth rate of an infinitesimal infection in an otherwise susceptible population. However, when the population is structured, the growth rate may depend on which class of individual is infected. We therefore allow the level of infection to equilibrate between the classes (such that high-risk individuals are more likely to be infected), before calculating the number of secondary cases. Correspondingly, for network models, R 0 should be calculated only once the early spatial correlations (which develop within a couple of generations) have formed". What we have shown is that this information about the early epidemic growth contained in R 0 computed in this non-trivial way is already included in the Jacobian matrix of pairwise models around the disease-free equilibrium. In other words, this result tells us how we must understand the linearization of deterministic pairwise models around a disease-free equilibrium.
Interestingly, replacing λ 1 by the early epidemic growth rate r in Eqs. (14) and (23) leads to the same linear and quadratic relationships between R 0 and r than those derived for homogeneous mixing SIS, SIR, and SEIR models [22,23,29,33,45]. However, one must be careful about how to apply these relationships when they are used to estimate R 0 , especially in populations with different levels of mixing, because greatly differing estimates of the initial epidemic growth rate can be obtained from the cumulative number of initial cases [5,27,28]. This typically happens when the primary cases are not representative of the population of interest because, for instance, they belong to close settings like hospitals or schools with very high contact rates [28,36].
When initial infections (primary cases) occur at random in a network model, the variability in the number of secondary cases produced by primary cases is expected to be higher than in the next generations of cases, i.e., once invading clusters have been built up. After the second generation of cases, individuals with the highest degree are very likely to be infected, and it is then when the estimates of λ 1 should be computed from the cumulative number of cases. This heterogeneity in the degree of the infected individuals is reflected in the time evolution of the effective reproduction number, R t , during an outbreak. The stochastic simulations of the evolution of R t during the exponential-growth phase of the epidemic for the SEIR model with rewiring show a very good agreement with the predicted value of R 0 for Poisson networks (see top panels in Figure 3). This good agreement is consistent with the one observed for the SIS model with rewiring in [17] for the same type of networks. For initial degree distributions with higher variance (exponential and scale-free contact networks), the accuracy of the R 0 estimates is not so good. However, for these networks and low rewiring rates, the stochastic simulations depict an initial increase of R t towards values close to (but lower than) the predicted R 0 , which is due to the infection of highly connected individuals, followed by a decrease of R t corresponding to the secondary infections caused by the latter. This behaviour was observed in the SARS and influenza A (H1N1) epidemics, where "super-spreading" events were reported [24,27,28].    , exponential (middle) and power-law (bottom) networks of size N = 10000 and mean degree 10, 10, and 9.86, respectively. In each realization, the epidemic starts by infecting 10 individuals selected uniformly at random. Results for two different networks are shown in all panels. Dashed line corresponds to the value of R 0 predicted by (18). Parameters: β = 0.05, γ = 1/7, θ = 3/5, ω 1 = ω 2 = 0 (left panels), and ω 1 = ω 2 = 0.1 (right panels).