A geometric analysis of the SIRS epidemiological model on a homogeneous network

We study a fast–slow version of an SIRS epidemiological model on homogeneous graphs, obtained through the application of the moment closure method. We use GSPT to study the model, taking into account that the infection period is much shorter than the average duration of immunity. We show that the dynamics occurs through a sequence of fast and slow flows, that can be described through 2-dimensional maps that, under some assumptions, can be approximated as 1-dimensional maps. Using this method, together with numerical bifurcation tools, we show that the model can give rise to periodic solutions, differently from the corresponding model based on homogeneous mixing.

a tractable trade-off between simplicity, which allows for more in-depth analysis, and realism, which allows to make more precise predictions.
In particular, compartment models build on the core idea that the population can, at any time, be portioned into compartments characterized by a specific state with respect to the ongoing epidemic. The first of such models divides the population into Susceptible, Infected and Recovered individuals, from which the SIR acronym is used. A Susceptible can become Infected (S → I ) by making contact with an already infected individual, and can then either Recover (I → R) or die, if we assume the disease to be characterized by permanent immunity after a first infection. If we do not make such an assumption, and allow recovered individuals to become susceptible again (R → S), we obtain a so called SIRS model. Many more models, with different compartments, have been proposed and analysed in the past, see e.g. Dafilis et al. (2012) and Hethcote (2000).
Classical compartmental models are based on the homogeneous mixing assumption, i.e. the assumption that any individual in a population may have contacts with any other. Such an assumption, however, is quite unrealistic for many situations in which the observed population is large, and possibly divided in classes, families or generally sub-populations. One possible extension is to subdivide the population into groups, assuming homogeneous mixing within each group, but representing inter-group interactions through a contact matrix (Mossong et al. 2008). Another possible approach is to take into account the network structure of contacts. Often, epidemic dynamics on a network is analysed only through simulations (López-García 2016; Smilkov et al. 2014;Zhang et al. 2013;Castellano and Pastor-Satorras 2010;Volz 2008;Ganesh et al. 2005). The method of pair approximations, introduced in epidemiology by Satō et al. (1994) and Keeling et al. (1997), allows to build a system of differential equations that retains some aspects of the network structure. The ideas and some applications of the methods are presented in detail in the monograph by Kiss et al. (2017). However, not much analytical progress has been made in the study of the resulting systems, possibly because they are generally rather complex.
This paper aims at introducing methods from Geometric Singular Perturbation Theory (GSPT) to analyse these systems, building on the ideas introduced in Jardón- Kojakhmetov et al. (2021). The difference in time-scales between epidemic spread and demographic turnover, which can be observed in many diseases, is the motivation for the use of techniques from GSPT. We refer to Jardón-Kojakhmetov et al. (2021) for a brief introduction of the techniques we use, or to the references therein, and in particular to Jones (1995) and Kuehn (2015), for a more detailed explanation. In particular, we will exploit the entry-exit function (De Maesschalck 2008;De Maesschalck and Schecter 2016) to analyse the behaviour of the system on its critical manifold, which is characterized by a change in stability over a hyperplane.
In this work, we assume homogeneity of the network, in order to obtain analytical results, before validating them numerically. Even with such an assumption, the additional complexity brought by the network structure must be treated properly. In fact, in order to completely describe the evolution of a network in time, one needs to have an equation for each possible state of its nodes, one for each possible state of its edges (along which the epidemic spreads), one for each possible state of triples, i.e. three nodes connected by two edges, and so on. This procedure, however, would generate very large system of ODEs, which would once again be hardly treatable with analytical tools. In order to overcome this difficulty, one can apply the so-called moment closure (Kiss et al. 2017;Kuehn 2016), i.e. approximation formulas which allow us to truncate the dimension of the objects we want to analyse. If we truncate at the node level, we lose the network structure, and we recover a homogeneously mixing system. Instead, we truncate at the edge level, using the pair approximation discussed above, and analyse the system which derives from this choice.
To our knowledge, there are relatively few articles in which GSPT has been applied rigorously to epidemics models (Rocha et al. 2016;Jardón-Kojakhmetov et al. 2021;Heesterbeek and Metz 1993;Zhang et al. 2009;Brauer 2019;Wang et al. 2014). However, for most infectious diseases, the presence of different time scales is natural. Moreover, though a SIR model on networks has been studied with moment closure already (Bidari et al. 2016;Kiss et al. 2017), the SIRS extension has not. Likewise, a thorough bifurcation analysis on compartment models such as the one we analyse in this paper is not present in the literature. The additional feature of the network structure, even in its most simplified version, i.e. homogeneous network, unravels new dynamics for the SIRS system we study. Indeed, we show that there exists a set in the parameter space which allows the system to exhibit a stable limit cycle, a situation that does not arise in the classical compartmental SIRS model. To complement the bifurcation analysis, we extend the geometrical argument from Jardón- Kojakhmetov et al. (2021) to the higher dimensional system we study, providing additional justification for the existence of stable limit cycles. It is worth noticing that the model we study is not globally in fast-slow standard form; as in Jardón-Kojakhmetov et al. (2021); Kuehn and Szmolyan (2015); Kosiuk and Szmolyan (2016), the fast-slow dynamics are only evident in specific regions of the phase space, in which a local change of coordinates brings the system to a standard two time scales form. In particular, we refer to the very recent monograph (Wechselberger 2020), in which the properties of singularly perturbed systems in non-standard form are thoroughly analysed.
This article builds on the analysis on a homogeneous mixing SIRS model we carried out in Jardón-Kojakhmetov et al. (2021), generalizing it to a more complex setting. The main novelty we introduce is the network structure, which increases both the dimensionality and the complexity of the ODE system we study. This additional feature, which increases the realism of the model, allows us to unveil additional dynamics, namely stable limit cycles. Lastly, a further challenge is posed by the fact that the rescaled system close to the critical manifold is characterized by multiple fast variables, as compared to a single fast variable which characterized all the systems studied in Jardón-Kojakhmetov et al. (2021).
The paper is structured as follows: in Sect. 2, we recall the derivation of the model, and introduce the moment closure technique. In Sect. 3, we obtain analytical results on the model, in particular on the fast and slow limit systems and on the application of the entry-exit function. In Sect. 4, we perform a bifurcation analysis and numerical exploration of the model. Finally, in Sect. 5, we conclude with a summary of the results, and with possible research outlooks.

Formulation of the SIRS model on a network
In this section we describe and propose an SIRS model for epidemics on graphs, building on the model proposed in Kiss et al. (2017, Sec. 4.2.2). We are interested in the graph generalization of the model studied in Jardón-Kojakhmetov et al. (2021), in order to drop the homogeneous-mixing hypothesis, under which we assumed that each individual in the population could have contacts with any other. We then assume loss of immunity to be slower, compared to the other rates (this is the case e.g. for pertussis (Dafilis et al. 2012;Lavine et al. 2011), and it could potentially be true for the recent SARS-CoV-2 (Kissler et al. 2020; Randolph and Barreiro 2020)); this assumption brings the model to a non-standard singularly perturbed system of ODEs, which we study with techniques from GSPT.

The model
The construction of the model is essentially what is presented in detail in Kiss et al. (2017, Ch. 4), extended to the SIRS case. For ease of reading, we briefly repeat the whole method.
We consider a network of N nodes, with N large, representing the individuals of a population, and we assume this network to be homogeneous, meaning that each node has fixed degree n ∈ N ≥2 , representing the number of direct neighbours each individual has. We assume the network to be undirected and completely connected, meaning that, given any two nodes in the network, there is a finite sequence of edges (or an undirected path) which starts in the first and ends in the second. Each node can be in one of either three states, namely S (susceptible), I (infected) or R (recovered). We will indicate the number of each state at time t with [·](t); we stress the distinction between the notation X , indicating a state, and [X ], indicating the number of individuals in the state X . We indicate the number of edges connecting a node in state X to one in state Y at time t with [XY ](t) for all t ≥ 0. We distinguish between an edge XY , counted starting from a node in state X , and the same edge counted starting from the other node in state Y , for a reason of conserved quantities, namely (7a), (7b) and (7c) to be defined below. For example, we count the number of edges S I by "visiting" each node in state S, and counting all its neighbours in state I , then summing over all the nodes in state S; this implies that, at all times, by definition, [S I ] = [I S]. The edges connecting a node with another in the same state, such as SS, hence, will always be counted twice.
Infection can only spread if a node in state S is connected to a node in state I through an edge S I ; we denote the infection rate with β ≥ 0. Nodes in state I recover, independently from their neighbours, at a rate γ > 0; and nodes in state R lose their immunity, again independently from their neighbours, at a slow rate , with 0 < β, γ . Based upon these modelling assumptions, it is then straightforward to prove using the master equation of the epidemic model, that one obtains the following system of ODEs: In order to fully describe the dynamics of the system, we need an ODE for [S I ] as well. To understand how the number of edges [S I ] evolve in time, we need to consider the role of triples, as exemplified in Fig. 1. A triple is a path of length 2 through a central node in state Y , connected to two nodes in state X and Z , respectively; we indicate such a triple with XY Z. The positions of X and Z are interchangeable, and the most important node is the central one, as we will explain shortly. The only change of the system which depends on the presence of a specific edge is the contagion which brings S I → I I . Direct neighbours of a node in the state S which get infected, i.e. the node X in a triple X S I , see their edge X S change to X I due to their belonging to the triple. The two other possible changes in the system, namely the recovery (a node in state I becoming R, which happens at a rate γ ) and the loss of immunity (a node in state R becoming S, which happens at a rate ) only happen at a node level, so the only nodes which see this change are the direct neighbours of the node changing state, and we do not need to consider their belonging to a triple.
For clarity, we fix a lexicographic order S ≺ I ≺ R for nodes and edges, and write the explicit equations for the edges which follow this order only. If we take into account all the triples with a central node in state S and at least one node I , which could infect the central one (as described in Fig. 2), we obtain the following system of ODEs, which describes the evolution in time of nodes and edges:

Fig. 2
Complete description of the edges dynamics considering edges and triples. Straight lines: infections; wobbly lines: recovery; dashed lines: loss of immunity. The base diagram is the same which appears in Kiss et al. (2017), to visually describe their SIR model; the new, slow dynamics in our model are the dashed blue arrows, symbolizing loss of immunity Notice the 2 which multiplies the right hand sides of edges connecting nodes in the same state: as we mentioned above, they are always counted twice, whether they are created or lost. To fully describe the system, we would then need to have ODEs for triples, quadruples, etc. Instead, we proceed as in Kiss et al. (2017), and apply moment closures.

Moment closures
Moment closure methods are approximation methods used in many contexts, in order to reduce large (or infinite) dimensional systems of equations to a smaller finite dimension (Kuehn 2016). Proceeding as in Kiss et al. (2017, Sec. 4.2), one can approximate the edges as functions of the nodes, or triples as functions of nodes and edges. If we choose the first option, assuming independence between the state of nodes, we can approximate all edges as follows: This implies that we lose the network structure and, up to rescaling the infection parameter byβ = nβ, we recover the SIRS system already studied in Jardón-Kojakhmetov et al. (2021).
Instead, in this work we choose to apply the second order approximation, and hence we approximate each triple with the formula given in equation (4.6) of Kiss et al. (2017), namely This approximation is based on the conditional independence between the states of neighbors of a node, using a counting argument, which for clarity we recall from Kiss et al. (2017

Analysis of the model
In this section we present the pair approximation SIRS model, and give our main analytical results. First, we are going to reduce the dimension of the system, exploiting three conserved quantities. Secondly, we introduce a formulation for the basic reproduction number for the system, and we describe the behaviour of the fast limit system. Next, we are going to derive the equilibria of the system in the biologically relevant region, and we show that the slow manifold of our perturbed system is exponentially close to the critical manifold. Lastly, we rescale the system in an O( )-neighbourhood of the critical manifold, with a scaling similar to the one proposed in Jardón-Kojakhmetov et al. (2021), and we apply the entry-exit formulation. Throughout the analysis, we notice that the parabola [SS] = n[S] 2 , i.e. approximation (4) applied to the edges in state [SS], on the critical manifold is of particular importance for the dynamics.

Fast-slow model
In this section, we derive the system we will study for the remainder of the article, applying moment closure to (3) and reducing its dimension.
Applying approximation (5) to every triple in system (3), we obtain the following singularly perturbed autonomous system in non-standard form: , in which, as from our assumptions, the processes of infection and recovery are fast, and the process of loss of immunity is slow. By construction, the sum of all the edges starting from a node in the state . This can be checked by carefully computing the difference of the derivatives of the right hand side(s) and the left hand side(s) of (7). By doing so, we reduce the dimension of the system, obtaining The basic reproduction number R 0 can be obtained (Kiss et al. (2017, p. 140)) for the limit as → 0 of system (8) as We notice that, for (9) to be well-defined and dependent on the parameters of the system, we need n > 2. The equality n = 2 describes the very special case of a ring network, i.e., a connected network in which all nodes have exactly two neighbours. In the remainder of the paper we assume R 0 > 1 and n > 2.

Remark 1
We notice that the threshold R 0 ≶ 1 in (9) is equivalent to since they all correspond to β(n − 2) ≶ γ . A formula corresponding to R 1 is given in Kiss et al. (2017), shortly after the definition of R 0 .
We notice that R 1 has a much more intuitive biological interpretation than R 0 . Consider a network with all the nodes in susceptible state S, except one in state I . Consider one of the n edges in state I S: this could either transition to RS, at a rate γ , and the epidemics would die out immediately, or spread the infection to the node in state S, at a rate β, and become an edge I I . If the latter happens, with probability β/(β + γ ), (n − 1) new edges move to state S I ; hence, R 1 can be interpreted in the classical meaning of "the number of edges infections caused by one infected edge in an otherwise susceptible population". Recall that the disease spreads only through edges S I (or I S, equivalently), so their number should be the quantity we measure in order to quantify the contagiousness of the disease; an edge I I can not be used to spread the disease. Now we compute the basic reproduction number R 1 for system (8) and > 0 sufficiently small.

Proposition 1
The basic reproduction number R 1 for system (8) is given by Proof We use the method first introduced in Diekmann et al. (1990), and then generalized in Van den Driessche and Watmough (2002) (see also Diekmann et al. (2010)).
We linearize system (6) with the matrix A given by There are clearly many ways of doing that, but the preferred splitting is such that M and V can be interpreted as the transmission (i.e. relative to new infections) and transition matrix (i.e. relative to any other change of state), respectively. Then, we compute where ρ indicates the spectral radius of a matrix. The choice for the two matrices is It can easily be checked, then, that V −1 has non-negative entries, and that, since M V −1 has two rows of zeros, This finishes the proof.

Remark 2
The perturbed R 1 given in (11) has a similar biological interpretation for the perturbed system to the one given for the corresponding R 1 (10) of the limit system as → 0. We need to compute R 1 , which corresponds to the average number of S I edges produced by an S I edge in a totally susceptible population. As in the previous case, an edge S I will become an edge I I with probability β/(β + γ ), producing in this case n − 1 edges S I . However, the original edge I I , after having become I R can become again an I S edge with probability /( + γ ). After having returned to the state S I , the edge will produce other R 1 S I edges, since the pairwise model does not consider higher order correlation and does not "remember" that the neighbours of S had already been infected once. Hence from which one obtains (12). Through this argument, we see that threshold for the SIRS model is different from the one for the SIR model, while in the homogeneous mixing case the two coincide.
The set is forward invariant under the flow of (8), for ≥ 0, so that solutions of (8) are global in time.
Proof Apparently the right-hand side of (8) has a singularity at [S] = 0; however, in the set , the terms [S I ]/[S] and [SS]/[S] are both bounded by n, so that the righthand side is indeed Lipschitz. Hence, system (8) has a local solution. Furthermore, it can be easily checked that the system is forward invariant by showing that the flow is pointing inwards on the boundary of . Hence, solutions of system (8) are global in time.

Fast limit
In this section, we study the fast subsystem (or layer equations) corresponding to the limit of system (8) as → 0 on the fast time scale. Hence, taking the limit → 0 in system (8), we obtain For ease of notation, we introduce In the fast dynamics, the susceptible population can only decrease, and eventually the infected population will not have any more susceptibles to "recruit" and will decrease as well. In particular, we prove the following: . From (14a) and (14d) we see that while from (14a) and (14c) Note that from (15) v is clearly decreasing for n > 2, and we see that Recall Lemma 2, which implies v ≥ 0; if v = 0, then [SS] = 0, and from Eq. (14c) we observe that [SS] will not change, so 0 is its corresponding limit value. Assume then We notice that we can rewrite (17) using (16) and obtain This means that which implies, recalling that which implies Similarly, using (15), we can show that which implies, using (20) and (21), and recalling that In particular, combining (21) and (22), we can write We notice, from (7a), that this implies that [S R] converges to a non-negative limit as well. Combining (14a) and (14b) (14) analogous to (23) holds for all t. Indeed, noticing that

Remark 3 A relation between [SS](t) and [S](t) for system
is a constant of motion for system (14), we observe that for any t ≥ 0 the relation holds.
The equilibria of the limit system are all of the form In particular, λ 5 changes sign on the hyperplane defined by We notice that β(n − 1) > 0, since we assume n > 2. Considering (26), we define the loss of hyperbolicity line on the critical manifold C 0 as We now give a closed formula for the value of [S] ∞ .

Proposition 3 Consider a generic initial condition
Proof We proceed as in (Bidari et al. 2016, Sec. 3). From our assumptions, (14d) and (24), we obtain Multiplying both sides by the integrating factor [S] Integrating (29) from t = 0 to t = +∞, and recalling that, by Proposition 2, Since, by assumption, the left-hand side of (30) is O( ), we ignore it, and we consider the right-hand side only. Hence, we find [S] ∞ by solving , from which we immediately obtain that [S] ∞ is given as a zero of the function H (x) defined in (28). Next, we prove that such a zero is unique.

Remark 6
The loss of hyperbolicity line can be rewritten as R 1 [SS]/(n[S]) = 1, meaning there is a minimum ratio of susceptible edges to susceptible nodes for the epidemic to "explode". In this form, the line can more straightforwardly be related to the threshold R 0 S = 1 in the classic S I R model. Moreover, since by (7a) we derive the constraint [SS] ≤ n[S], if R 1 < 1 the epidemic can never start, as the loss of hyperbolicity in that case lies in an unreachable region of the critical manifold.

Equilibria of the perturbed system
The following proposition discusses the equilibria of system (8).
Proposition 4 For > 0 sufficiently small and R 0 > 1, system (8) has 2 equilibria in the relevant region of R 5 . Disease free equilibrium: Endemic equilibrium: to their first order on the components are given by: Proof The disease free equilibrium is trivial. The endemic equilibrium is computed by expanding the variables in power series of , e.g.

Remark 7
Since we assume R 0 = β(n−2) γ > 1, recall Remark 1 and (11) We notice that the disease free equilibrium belongs to C 0 defined in (25), and by computing the corresponding  (27); hence, it approaches it as → 0. See also Lemma 5 below.

Slow manifold
In this section we provide a multiple time scale description of the disease-free, or near disease-free states.

Proposition 5
The slow manifold of system (8) is exponentially close in to the critical manifold C 0 given by (25).

Proof
The invariant manifold C 0 is an invariant manifold also for system (8) Recall that [S] ∞ and [SS] ∞ are the initial conditions for the slow flow. Solving (32) explicitly yields Lemma 3 The parabola , given by (34), is uniformly attracting for the slow reduced subsystem (32).
Proof Recall (33). The distance between a solution curve of (33) and can be parametrized by the slow time τ . With this in mind, let d(τ ) denote such distance, we then have Proof Notice that, considering Lemma 3, the assumption of starting close to the parabola is not restrictive. The derivation of G(x) is analogous to the derivation of H (x) of Proposition 3, using
In Fig. 4 we compare formula (35) and direct integration of the layer system (14) starting with a small fraction of infected nodes. (23). Since we showed that the parabola (34) is attracting for the slow flow, we can assume that, after the first slow piece of any orbit, [SS] 0 = n[S] 2 0 + O(δ 1 ), where 0 < δ 1 1. We can then rewrite (23) as

Remark 8 Recall
where the ≈ symbol indicates an O(δ 1 ) error. For n large enough, the last factor is close to 1, and the entry point for the slow flow is approximately on the parabola. This rescaling brings the model, after rearranging the variables, to a singularly perturbed system of ODEs, namely [S] , which can be rewritten in a standard form, and rescaled to the slow time scale, denoting now the time derivative with respect to the slow time parameter with an overdot, giving Taking now the lim →0 of (37), we obtain the system of algebraic-differential equationṡ The last three equations of (38) are satisfied for [v] = [Sv] = [vv] = 0. This is exactly the critical manifold of (8), on which the dynamics is described by (32). Using (32), we can show how λ 5 changes in time, in the slow flow, by deriving its formulation (26) with respect to time, obtaininġ

This implies that λ 5 is increasing if [SS] < α([S]), where the function α is defined by
In Fig. 5 we visualize the behaviour of two orbits in the slow dynamics. We note that, even if an orbit enters the slow flow in a point below the purple line but above the green curve, i.e. in the region whereλ 5 < 0, it eventually has to cross the green line before crossing the purple curve, since they represent respectivelyλ 5 = 0 and λ 5 = 0. Hence, any orbit will eventually evolve in the regionλ 5 > 0. We prove the following:

Proposition 6 The subset {([S], [SS]
) ∈ (0, 1) × (0, n)|λ 5 > 0} is forward invariant for system (32). If we take the scalar product of ν with the vector field F given by (32), we obtain meaning that on the curve [SS] = α([S]), this scalar product is negative, hence orbits approaching the curve from below will not cross it. (34) and (39), we notice that the curve α is always above the parabola ; hence, by invariance of and Proposition 6, an orbit starting above the parabola will eventually be "squeezed" between α and .

Remark 9 By comparing
The following Lemma is also insightful: The previous Lemma qualitatively tells us that for sufficiently small and large enough n one expects the endemic equilibrium to be near the intersection point

([S] * , [SS] * ). In this case, if the endemic equilibrium is stable, then ([S] * , [SS] * )
is a good approximation of such an equilibrium, while if there are limit cycles then these are roughly "centered" at ([S] * , [SS] * ). Notice that, in the latter case, limit cycles do not have to be "centered" at ([S] * , [SS] * ) if n is not large enough. See more details in the forthcoming sections.

Entry-exit function
Dividing the last three equations of system (36) by on both sides, we obtain System (40) can be rewritten as where we denote x := [S] [SS] , z := [Sv], and w := [v] [vv] . The critical manifold C 0 = {z = 0, w = 0 0 } is invariant for system (41) both when > 0 and = 0. Recall (26); it is clear that g(x, 0) = λ 5 ≶ 0 when x ∈ C A 0 or x ∈ C R 0 , respectively. To control the relation between the starting point of the slow dynamics and the transition point back to the fast dynamics, we are going to employ the entry-exit function (De Maesschalck and Schecter 2016;Liu 2000;Schecter 2008). We cannot apply (Liu 2000, Theorem 2.5), as was done in Jardón- Kojakhmetov et al. (2021), since that requires that the fast dynamics is one-dimensional, while (41) has three fast variables. Instead, we can apply Theorem 4.7 of Liu (2000); see also Schecter (2008, Theorem 2.4), noting that the map 0 was defined at p. 417. We remark, however, that one of the hypotheses necessary for the application of such formula is the separation of the negative eigenvalues from the eigenvalue which causes the loss of stability. To be more precise, let us recall that the non-trivial eigenvalues of the layer equation on C 0 are λ 3 = −2γ , λ 4 = −γ and λ 5 , see (26). Therefore, the exchange lemma of Theorem 4.7 of Liu (2000) can be applied only to trajectories contained in the portion of C 0 in which λ 5 > −γ . Recalling (26), we have that λ 5 > −γ if and only if Thus, let us now assume that a solution of (32) satisfies (42) we can implicitly compute the exit time T E of an orbit on the slow manifold through the integral which is simply a rewriting of the entry-exit integral (Liu 2000). Figure 12a, b show two trajectories satisfying (42) for which formula (44)  where the last (approximate) equality is obtained as in Remark 8. It is possible that (42) does not hold for some τ > 0 even if holds at τ = 0, but following similar steps as we have just described, one can show that this also happens only if [S] ∞ is sufficiently close to 0. Hence (42) holds for all τ ≥ 0 whenever [S] ∞ is large enough. We recall that [S] ∞ can be computed as zero of the function H (x) given in Proposition 3. Accordingly, let us recall from Proposition 3 that Thus, we see that for x < [S] 0 , H changes proportionally to β. Then, because H increases as β also increases, we can use the implicit function theorem to argue that [S] ∞ is a decreasing function of β, for any fixed [S] 0 ∈ (0, 1). The previous fact can also be seen from (21). Hence, from the arguments described above, and recalling that β > γ n−2 so that R 1 > 0, we conclude that there exists β * such that (42) is satisfied for all τ if γ (n−2) < β < β * . On the other hand, (42) would not hold if β is large enough. In that case, the formula (44) does not provide good approximations for the exit point, see more details in Appendix B. From now on we shall assume that (42) holds.

Lemma 6 The exit time T E is finite for any initial point ([S]
Proof Recall (45). For small positive values of τ , λ 5 (τ ) < 0, since the slow dynamics begins in the attracting region C A 0 . Hence, for small values of τ ≥ 0 the integral τ 0 λ 5 (σ )dσ < 0.
From (46), we observe that hence there exists at least one finite T E which satisfies (44). From our previous analysis, we know that λ 5 (τ ) = 0 only once during the slow flow, and it remains positive afterwards; hence, such T E is unique.

Application of the entry-exit formula to the parabola
As we have remarked so far, the parabola (34) is of particular interest for the dynamics, even more so for large values of n. Hence, we are interested in understanding the entryexit relation on this specific invariant set. We now consider the evolution, under the slow flow, of the point ([S] ∞ , [SS] ∞ ) = (0, 0); with these initial conditions, (33) becomes Being able to write [SS] as a function of [S] allows us to compute the exit point for the origin, which in general is not possible, since λ 5 depends on both slow variables.
Combining (47) and (45) we obtain which can be equivalently rewritten, introducing for ease of notation C := ((n − 2)β − γ )/(β(n − 1)), as  Fig. 6 for a sketch of the function h, and a visualization of the argument of this proof). We observe that h(0) = 1 and h(1) = 0. Deriving h(x), we see that The study of the asymptotic behaviour of system (8) is then reduced to two 2-dimensional maps, from C 0 to itself; specifically, we define 1 ( . We now explain the rationale behind introducing the one-dimensional maps 1 and 2 obtained by restricting the domain and approximating the range of 1 and 2 to the parabola (see Fig. 7). Indeed, Remark 8 and Lemma 3 show that the parabola is close to being invariant for the map 2 • 1 and it is well known that the occurrence of near one-dimensional return maps is an important theme in multiple time scale systems (Bold et al. 2003;Guckenheimer et al. 2006;Kuehn 2011;Medvedev 2005).
Next, consider a point with [S] coordinate [S] 0 , O( ) away from the parabola (34), in the repelling part of the critical manifold. Its image [S] ∞ under the fast flow, which defines the map 1 sketched in Fig. 7, is given by (35). We notice that this value depends on both β and γ , as well as on n. For n large enough, the entry point in the slow flow will be close to the parabola, as argued in Remark 8; hence, we will be able to compute its exit point [S] 1 using (49), which again depends explicitly on all the parameters of the system in a highly non-trivial way. This is different from the SIRWS model studied in Jardón-Kojakhmetov et al. (2021), in which there was a clear separation between fast parameters, which dictated the fast dynamics, and had no influence on the slow one, and slow parameters, which characterised the viceversa. The map 2 in Fig. 7 sketches the relation between the entry point [S] ∞ and its corresponding exit point [S] 1 , i.e. (50).
Depending on the relative position of [S] 0 and [S] 1 , we might be able to deduce the asymptotic behaviour of the system. However, the high dimensionality of the layer  (35)) hinders the analysis of the system with non-numerical tools.
We proceed now to a bifurcation analysis of system (8), and finally, with a technique similar to the one detailed Jardón- Kojakhmetov et al. (2021, Sec. 3.4.1), to numerically investigate the existence of periodic orbits by concatenation of fast and slow pieces. We stress the versatility of the numerical argument we present, which is similar to the one we used in Jardón-Kojakhmetov et al. (2021), applied now to a higher dimensional system.

Bifurcation analysis and numerical simulations
In this section, we carry out a bifurcation analysis for the behaviour of system (8), which will then be verified by numerical simulations and by a geometrical argument. Bifurcation analysis is done on system (8), which for small values of is stiff (as we showed in Proposition 5, the slow manifold is exponentially close to the critical manifold), while the numerical simulation concerns a combination of systems (14) and (32), which are both non-stiff.
It is important to notice that, even though the layer system (14) (46) is fundamental, in this setting, to carry out a meaningful numerical exploration of the model.
Without loss of generality, we set γ , which is the inverse of the average infection interval, to 1; this simply amounts to an O(1) rescaling of time, and we rescale the other parameters accordingly, keeping however the same symbols, for ease of notation. System (8) then has only three parameters, namely , n and β.
Using MatCont (Dhooge et al. (2008)), we are able to completely characterize system (8) through numerical bifurcation analysis. We only consider the first octant of R 3 , for the biological interpretation of the parameters. Numerical analysis shows the existence of a Hopf surface , whose "skeleton" is depicted in Fig. 8. For values of the parameters between the plane = 0 and , the system exhibits a stable limit cycle, while for values above , the system exhibits convergence to the endemic equilibrium (31). Our bifurcation analysis suggests the existence of a value * ≈ 0.18 such that, for > * , the system only exhibits convergence to the endemic equilibrium, regardless of the values of β and n. To make Fig. 8 more readable, we provide intersections of n β Σ Fig. 8 A skeleton of the bifurcation surface . Green (respectively, red and blue) curves correspond to constant values of (respectively, β and n). We notice that, for values of n ≥ 6, system (8)   the surface with some planes n = k (Fig. 9a), β = k (Fig. 9b), and finally = k (Fig. 10).
As in Jardón-Kojakhmetov et al. (2021), we see an expansion of the parameter region which exhibits stable limit cycles as decreases, see Fig. 10. This means that, as decreases, i.e. as the ratio between the average lengths of the infectious phase and the immunity interval decreases, we are more likely to observe occurrence of stable limit cycles in the disease dynamics. We do not observe, however, a divergence in the n direction, as the limit as → 0 of the surface contained in the green curves of Fig. 10 is still bounded.
In order to verify the accuracy of the surface , we investigate the system via a numerical implementation of the same geometrical argument used in Jardón-Kojakhmetov et al. (2021, Sec. 3.4.1). There, we numerically showed the existence of a candidate orbit by concatenating heteroclinic orbits of the layer equation, from the critical manifold to itself, and orbits of the slow flow, truncating each at the corresponding exit time. The system studied in Jardón-Kojakhmetov et al. (2021) was 3-dimensional, but the slow flow evolved on a 2-dimensional plane in R 3 ; as we showed Once we have obtained the candidate value [Ŝ], we define a small interval J 1 in the [SS] coordinate around [ŜS]. Next, the interval J 2 is obtained as the evolution of J 1 under the layer equation, that is J 2 = 1 (J 1 ). Analogously, we then obtain the interval J 3 as the evolution of J 2 under the slow flow for a time precisely given by T E , that is J 3 = 2 (J 2 ). In Jardón-Kojakhmetov et al. (2021) we have shown that if J 3 intersects J 1 transversally, then the perturbed system, for > 0 small enough, exhibits a limit cycle. Furthermore, if the map = 2 • 1 is a contraction, then the limit cycle is locally stable. Here we provide an additional numerical illustration of our argument, see Fig. 11. The black solid and dashed lines depict, respectively, a choice of J 1 and the corresponding J 3 . We see that these two intervals intersect transversally. This hints to the possibility of a singular cycle with exit point close to the intersection point. To be more specific, when taking an interval J 1 corresponding to [S] [ŜS]) such that ( p) = p. That is, p is a fixed point of the map . In this way, the concatenation of an orbit of the layer system that is backward asymptotic to p and an orbit of the slow subsystem with exit point at p is what we call a singular cycle. Regarding stability, since we have shown that the slow flow is a contraction, and that the flow of the layer system is not expanding, we can conclude that the singular cycle is attracting. Figure 12a, b depict the numerical realization of the two limit systems for a choice of (n, β) for which we do expect limit cycles: indeed, J 1 and J 3 intersect transversely, and bifurcation analysis confirms that, for this choice of the parameters, there is a stable limit cycle for > 0 sufficiently small. Figure 12c, is the projection on the ([S], [I ])-plane of an orbit of system (8), starting from a random initial condition. As we expected, for sufficiently small, the perturbed system exhibits a stable limit cycle, as argued from the limiting situation depicted in Fig. 12 a, b.
Moreover, we remark on the necessity of being "small enough", with this condition being represented by the surface , for the above argument to hold. In particular, in Fig. 13, we show how increasing from 0.01 to 0.02 destroys the limit cycle at corresponds, instead, to convergence to the endemic equilibrium. Hence, for this specific realization, = 0.01 is small enough, whereas = 0.02 is already too large. The threshold for is given by the intersection of the surface and the line {n = 3, β = 1.2} in the (n, β, )-space represented in Fig. 8.

Discussion and outlook
We have analysed the behaviour of an SIRS model for epidemics on networks, exploiting Geometric Singular Perturbation Theory (GSPT), similarly to what performed in plane of a numerical simulation of system (8) from a random initial condition, exhibiting convergence to a stable limit cycle predicted by our geometrical argument (color figure online) Jardón-Kojakhmetov et al. (2021), in a way appropriate to a system in a nonstandard singularly perturbed form. From the point of view of applications, the main result found, through bifurcation analysis and geometric numerical arguments, is that, for a significant open subset of the parameter space, the network model exhibits stable limit cycles. This result is in sharp contrast to the global convergence to equilibrium of solutions to the standard (random mixing) SIRS model (Hethcote 1976;O'Regan et al. 2010;Jardón-Kojakhmetov et al. 2021).
Stable periodic cycles are found the value of n, the number of neighbours every individual has, is between 3 and 5 included. It is not surprising that a small value of n is needed, since for large n the model approximates a random mixing model, which, as stated above, exhibits global convergence to equilibrium. On the other hand, n ≤ 5 may appear a strongly restrictive condition for actual contact networks; however, it must be considered that in the model all contacts are assumed to have the same (a) (b) (c) (d) Fig. 13 Numerical illustration of the geometrical argument, with a particular focus on the effect of increasing on the dynamics. a, b Correspond to the singular flows as in Fig. 12, here for the choice of parameters β = 1.2 and n = 3. c, d The projections on the ([S], [I ]) plane of two numerical simulations of system (8) starting from a random initial condition, for = 0.01 and = 0.02, respectively. Notice that since the intervals J 1 and J 3 intersect transversely, we expect to find a stable limit cycle for > 0 sufficiently small. This is indeed true for = 0.01, but no longer true for = 0.02. In other words, for this choice of parameters, = 0.02 is not sufficiently small so that our geometric argument holds strength, while in reality individuals may have a limited number of strict contacts (say, three to five very close friends) plus a number of infrequent contacts. An interesting question would be whether the same result would hold with such a network model with weighted connections. Thus, it is clear that there is a strong motivation to use techniques from GSPT to investigation more complex network models, since the current analysis unveiled asymptotic behaviours which are impossible in the corresponding system under the random mixing assumption.
The results shown here were obtained for a deterministic system representing the pair approximation (Kiss et al. 2017) of a stochastic network model. A natural question is whether the periodic solutions, identified in the pair-approximation model, correspond to a recognizable pattern in the network simulations as well. Preliminary numerical explorations seem to suggest that indeed network simulations fluctuate around a closed orbit when the deterministic model predicts an attracting periodic solution. We plan to pursue further these investigations, that are however hampered, for particularly small, by the necessity of simulating a very large network to avoid the infection dying out in the long time orbits spend near [I ] = 0.
This last observation raises a general issue concerning the applicability of the model. Our analysis (Lemma 6) shows that the slow flow takes a finite time (the time passing between two consecutive outbreaks of the infection) in the slow scale, hence of the order O(1/ ) in the fast time scale. During that period the value of I (representing the fraction of infected nodes) is exponentially (in ) close to 0; hence, in a more realistic stochastic model, complete extinction of the infection would be likely, unless population size N is very large or there exists the possibility of external reintroduction of the infection. In the context of stochastic models, it would probably be interesting considering joint limits in which → 0 while N → ∞, as in the discussion of critical community size in Diekmann et al. (2013).
We stress the versatility of our geometric procedure, which gives us a numerical intuition of the asymptotic behaviour of a stiff system, i.e. system (8) with 0 < 1, without having to actually integrate it, but through simple integration of the corresponding two non-stiff limit systems, which we derived through the use of GSPT. This is particularly important for the high(er) dimensionality of the system, which hinders analytical results on the perturbed system. In particular, the same strategy is likely to generalize to more complicated network-based ODE models derived from moment closure.
We must acknowledge that one of the main tools used in the analysis, the entry-exit map, requires an extra assumption (regarding the order of the negative eigenvalues of the critical manifold) to ensure that the fast dynamics near the critical manifold can be approximately one-dimensional. The assumption holds if the contact rate β is not too large, but fails for large values of β. When this happens (see Appendix B), the formula (44) to predict the exit point is no longer correct, but we see numerically that periodic solutions still exist in that case too. This stresses the need to develop methods to be able to deal with the case when the fast dynamics is not one-dimensional.
As mentioned above, we expect that for n → +∞ the system approaches the simple epidemic model with random mizing; however, this has not been rigorously established. One intermediate step between having two independent perturbation parameters (i.e., → 0 and n → +∞) could be to couple n and , for example taking n = O(1/ α ), for some α > 0. However, this goes beyond the scope of this project, and we leave this as a prompt for future research. by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

(a) (b)
Fig. 14 a Evolution under the layer system (red) of a small interval J 1 . Its image defines the entry interval J 2 on the critical manifold. The evolution of each point of J 2 under the slow flow, stopped at the exit time provided by formula (46), is shown in blue. The evolution of J 2 defines the exit interval J 3 . The black line represents the border of (42): notice that the blue orbits partially lie below that line. b Zoom on the relative position of J 1 and J 3 on the critical manifold. Since these intervals intersect transversely, our geometrical argument predicts a limit cycle with exit point O( ) close to [S] ≈ 0.77 we provide numerical evidence that the formulas derived in Sect. 3.6 fail to reliably predict the exit point whenever the above assumption does not hold. We start by considering a pair of parameters (β = 16, n = 3) for which our bifurcation analysis predicts that orbits converge to a stable limit cycle. In Fig. 14, we depict the same procedure as in in Fig. 12, and we indicate with a black line the set [SS] = n n−1 [S] . Notice that J 2 , and part of the orbits corresponding to the slow flow, lie in the region defined by [SS] < n n−1 [S]. That is, for these orbits the important above assumption does not hold. Therefore, the exit time presented in Sect. 3.6 may not be correct, as we indeed show below.
We now compare the theoretical result we would obtain by applying the entry-exit function to numerical simulations of the perturbed system (8), shown in Fig. 15. The simulation confirms the presence of a stable limit cycle, but the exit point (the largest [S] value in the cycle) is not O( ) close to the one predicted by the entry-exit formula ([S] ≈ 0.77). These numerical simulations suggest that the straightforward application of the entry-exit formula to orbits crossing the line [SS] = n n−1 [S] provides wrong numerical estimates of the exit point. This is likely due to the fact that the ordering of the eigenvalues is not preserved along the orbits. Thus, we conjecture that the interplay between the two attracting but constant fast directions and λ 5 has a non-trivial role in the entry-exit map.
A first step towards resolving the issue described in this Appendix would be to understand how the ordering of the multiple non-trivial eigenvalues of the critical manifold influences the dynamics and especially what the corresponding entry-exit function would be.
(a) (b) (c) Fig. 15 Simulations of system (8) for n = 3, β = 16, random initial conditions, and decreasing values of . In this case the exit point for the limit cycle is [S] ≈ 0.58, which is not close to the one predicted by our entry-exit formulation, compare with Fig. 14