{'fraction_non_alphanumeric': 0.051930647457911046, 'fraction_numerical': 0.014797442555211213, 'mean_word_length': 4.069780608634112, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 26, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A four compartment epidemic model with retarded transition rates

We study an epidemic model for a constant population by taking into account four compartments of the individuals characterizing their states of health. Each individual is in one of the compartments susceptible (S); incubated - infected yet not infectious (C), infected and infectious (I), and recovered - immune (R). An infection is 'visible' only when an individual is in state I. Upon infection, an individual performs the transition pathway S to C to I to R to S remaining in each compartments C, I, and R for certain random waiting times, respectively. The waiting times for each compartment are independent and drawn from specific probability density functions (PDFs) introducing memory into the model. We derive memory evolution equations involving convolutions (time derivatives of general fractional type). We obtain formulae for the endemic equilibrium and a condition of its existence for cases when the waiting time PDFs have existing means. We analyze the stability of healthy and endemic equilibria and derive conditions for which the endemic state becomes oscillatory (Hopf) unstable. We implement a simple multiple random walker's approach (microscopic model of Brownian motion of Z independent walkers) with random SCIRS waiting times into computer simulations. Infections occur with a certain probability by collisions of walkers in compartments I and S. We compare the endemic states predicted in the macroscopic model with the numerical results of the simulations and find accordance of high accuracy. We conclude that a simple random walker's approach offers an appropriate microscopic description for the macroscopic model.


Introduction
The origin of modern epidemic modeling started with the seminal work of Kermack and Mc Kendrick almost a century ago [1]. They introduced the so called 'SIR model' where SIR stands for three compartments characterizing the states of health of an individual, namely S = susceptible, I = infected, R = recovered. The classical SIR model and various generalizations are able to capture some of the features of epidemic spreading of infectious diseases as observed in influenza, measles, mumps, and rubella. In the meantime a huge burst of compartment models and generalizations of the SIR model have been introduced [2,3].
While the interest in mathematical modeling of epidemic spreading was growing continuously, it is unsurprising that the emergence of the Covid pandemics has launched an enhanced interest and urgent need in advanced epidemic modeling [4]. Many models consider a set of nonlinear ordinary differential equations (ODEs) for the time evolution of the compartment populations where a new direction is about to emerge by combining these models with approaches inspired by network science [6,7,5,8,9,10,11] and fractional dynamics [12,13,14,15,16].
It has turned out that in many cases epidemics including Covid exhibit quasi-periodic patterns and spontaneous new outbreaks even after longer times of inactivity. Persistent oscillatory characteristics in epidemic dynamics was already pointed out a long time ago in the work of Soper [17]. Indeed, a major drawback of classical SIR type models without memory effects lies in their incapacity to capture (persistent) oscillatory behavior.
The present paper aims to tackle this issue and is a generalization of recent works [10,11]. In [11] is introduced a SIRS compartment model which takes into account a random duration of immunity R allowing a delayed transition R → S. What is found that this model is able to capture persistent oscillatory behavior in the number of infected individuals. This feature comes along as oscillatory Hopf instability of the endemic state. In the first part of our paper we extend this model to four compartments as relevant states of health in order to capture a larger variety of epidemics. We consider the four compartments S -susceptible; C -infected but not yet infectious (in the incubation phase); I-infected and infectious (ill); R-recovered and immune. We assume random sojourn times (waiting times) t C , t I , t R an individual once infected spends in each of the compartments C, I, R. The sojourn time t C in compartment C is interpreted as incubation time, i.e. the delay between infection and the outbreak of the disease. The sojourn times t I is the duration of the disease (infected and infectious state) and t R indicates the immunity period after recovery. It appears natural to assume that these variables are not fixed constants, but individually fluctuating random variables drawn by specific distributions. With these assumptions we derive the SCIRS evolution equations where we focus on several pertinent situations of waiting times with existing and non-existing mean. For waiting time probability density functions (PDFs) with existing mean we derive explicit formulae for the endemic equilibrium and analyze its stability.
In the second part of the paper we combine our SCIRS model with a multiple random walker's approach which is implemented into computer simulations. We present a case study and give numerical evidence that a simple random walk approach ('Brownian motion') offers an appropriate description of the SCIRS dynamics with memory.

Four compartment model 2.1 Basic notions
We consider a constant population of indicate their numbers in the compartments S C I R at time t. We consider here first a continuous time model for Z ≫ 1 walkers with the compartment population fractions We allow for each walker the following transition pathway: Intuitively we infer that allowing for persistent oscillatory and (quasi-) periodic behavior it is necessary to have a closed (cyclic) pathway S → . . . → S of transitions. We assume that t C , t I , t R ≥ 0 are mutually independent random waiting times an individual spends (after infection) in the compartments C, I, R drawn from specific probability density functions (PDFs) which we prescribe by the kernels K C (t), K I (t), K R (t), respectively. These kernels have to be causal functions 1 with (read P as 'probability'), for instance indicates the probability that the incubation time t C ∈ [τ, τ + dτ ] which is clearly non-zero only for τ ≥ 0 as t C ≥ 0. In other words: Given a walker has entered compartment C at τ = 0, then K C (τ )dτ indicates the probability that this walker leaves compartment C during the infinitesimal interval [τ, τ + dτ ] (by a transition C → I). We introduce the infection rate A(t) = A(j(t), s(t), t) ≥ 0 (number of infections per walker and unit time or entry rate to the incubated compartment C at time t having of units sec −1 ). We point out that A(t) is not a known given function of time t but rather an implicit non-linear function of s(t), j(t), namely This quantity contains macroscopic information on the instantaneous number of collisions of individuals in the 'reactive compartments' S and I and needs to be specified by model assumptions. It is important to point out that A(t) does not depend on the fractions c(t), r(t) of the 'invisible' compartments C and R. We assume here the most simple nonlinear form A(t) = βj(t)s(t) where β > 0 is a constant independent of time where β −1 introduces a characteristic time scale. β is also a macroscopic measure for the infection probability in a collision of an S and an I walker. The microscopic picture which we will invoke subsequently and implement into computer simulations indeed is that the infection rate is driven by random collisions between infected and susceptible independent random walkers.
With these simple assumptions we can establish model equations governing the time evolution of the fractions s(t), c(t), j(t), r(t). We consider the dynamics starting at time t = 0 with an initial condition s(0) = 1 − j 0 , and j(0) = j 0 (many healthy and a few infected walkers) and no incubated or recovered walkers c(0) = r(0) = 0. Our (SCIRS)-model evolution equations have the following general structure The sum of these rates is vanishing due to s(t) + c(t) + j(t) + r(t) = 1. In order that an epidemics starts, it is necessary that the globally healthy state becomes unstable when the initial state is close to the healthy state with j 0 ∼ 1 Z → 0+, i.e. where for instance only one infected walker is present at t = 0 (among a large population Z). We analyze the stability of the healthy state subsequently in detail. In (1) we have employed the notation for convolutions of causal functions K(t), f (t). Be aware that convolutions commute and are associative which can be seen by the representation (t ≥ 0) of the multiple convolution where δ(. . .) indicates the Dirac δ-function and (3) is non-zero only for t ≥ 0 (causality of the K i (τ )). Introducing the Laplace transform of our kernels which we define aŝ where λ indicates the Laplace variable. K(λ) λ=0 = 1 indicates that the kernels are normalized PDFs. A further observation is worthy of mention. Any (multiple) convolution of PDFs gives again a PDF. This can be easily seen by integrating (3) over time to yield one. Therefore the kernels K C (t), (K C ⋆ K I )(t), (K C ⋆ K I ⋆ K R )(t) which appear in (1) all are (normalized and causal) PDFs representing, respectively the densities of the random variables t C , t C + t I , and t C + t I + t R . We can see in (1) that for each compartment where an individual remains for a random duration an additional convolution occurs describing the delayed transition out of the compartment. These delayed transitions introduce memory into the rates of individuals leaving the compartments. We will give a careful account for this issue later (subsequent representation (11)).
The interpretation of (1) is as follows. In the first line A(t) is the rate of infections (transitions S → C) at time t. The second term (A ⋆ K C ⋆ K I ⋆ K R )(t) is the rate of individuals losing their immunity (R → S) having undergone the full pathway of SCIRS transitions. In the second Eq. of (1) A(t) re-appears as the entry rate to state C. The convolution (A ⋆ K C )(t) is the rate of delayed transition C → I of individuals that fall ill at time t. This rate reappears in the third line of (1) as entry rate into compartment I. Then (A ⋆ K C ⋆ K I )(t) is the rate of recovery I → R reoccurring in third equation as entry rate into compartment R. Finally (A ⋆ K C ⋆ K I ⋆ K R )(t) is the rate of individuals losing their immunity (transition I → S) and re-appears in the first equation as entry rate into S. Let us consider for a moment the second equation in (1) more closely d dt For our convenience we introduce the 'survival probability' Φ C,I,R (t) for individuals in compartments C, I, R (also called 'persistence probability' or 'survival function'), e.g. [14] (and many others) indicating the probability that an individual which is entering compartment C, I, R, respectively at t ′ = 0 at time t ′ = t still is ('survives') in this compartment. (5) is capturing all realizations with t C,I,R > t where d dt Φ C,I,R (t) = −K C,I,R (t). We observe in this relation the initial condition Φ C,I,R (0) = 1 as a consequence of the normalization of the waiting time PDFs and is telling us that an individual is, with probability one, in a compartment at the instant when entering it. Further we have Φ C,I,R (t → ∞) → 0+, i.e. an individual 'survives' only a finite time in compartments C, I, R to enter eventually the susceptible state S. We hence can integrate (4) and rewrite as where initial condition c(0) = 0 is assumed. We can also verify this relation by its Laplace transform (subsequent second equation of (9)) and directly by differentiating this expression with respect to t and using the initial condition Φ C (0) = 1. We interpret (6) as follows: A(τ )dτ is the fraction of walkers entering C during [τ, τ + dτ ] and Φ C (t − τ ) is the (survival-) probability that this fraction still is in C after a delay of t C = t − τ , i.e. at time instant t where this expression sums up over the complete history of entries into C from 0 to time t. With these considerations we can rewrite our SCIRS Eqs. (1) in the equivalent integral form where the initial conditions s(0) = 1 − j 0 , c(0) = 0, j(0) = j 0 r(0) = 0 are assumed.

Memoryless case -exponentially distributed waiting times
It is worthy to consider exponential waiting time kernels K C,I, Substituting this into (1) and accounting for (7) takes us to the particular simple form without memory Indeed exponential densities stand out by the memoryless feature and the Markov property (see textbooks, e.g. [18] and many others for details), thus the transition rates in (8) depend only on the instantaneous state s(t), c(t), j(t), r(t) but not on the history of the evolution. Setting the transition rates on the left hand side to zero yields straight-forwardly the endemic state which we derive subsequently (Eqs. (21), (22)) for arbitrary waiting time kernels with existing means.

Arbitrary waiting time distributions with memory
To explore the general cases with memory consider now the equivalent representations (7) and (1) in the Laplace spacê where all Laplace transforms (LTs) depend onÂ(λ) of the (unknown) infection rate and are the LTs of the survival probabilities. We further have that as we deal with a constant population. We observe that in the first Eq. of (9) the term is the LT of the survival probability P[t C + t R + t I > t] (due to the PDF (K c ⋆ K I ⋆ k R )(t) of the random variable t C + t I + t R ), corresponding to individuals 'surviving' in one of the compartments C, I, R at time t. To shed more light on how the memory comes into play we introduce a 'memory operator' (For our convenience we employ here a slightly modified definition as in [14,26] Then readingK(λ) = M(λ)Φ(λ) as convolutions in the time domain we can rewrite our SCIRS Eqs. (1) (in the Caputo sense) of general (fractional) derivatives (see [29,15,16] and Appendix) as d dt These equations show that the transition rates at time t have a complete memory of their previous values. Alternatively we can write (11) in the Riemann-Liouville manner of general fractional derivatives (see Appendix). For exponential waiting time PDFs K exp (t) = ξe −ξt we haveM exp (λ) = ξ thus M exp (t) = ξδ(t) (i.e. is null for t > 0 indicating lack of memory) and we easily recover above memoryless SCIRS Eqs. (8).

Long waiting times -the time fractional case
Representation (11) is especially useful when we deal with fat-tailed waiting time PDFs without existing means such as the Mittag-Leffler waiting-time PDF (ML-PDF) which (see [28] and references therein for explicit formulas). The ML-PDF is a fractional generalization of the exponential PDF and has the LT [14,26] (and many others) where due to the fat tail the mean waiting time − d dλK β (λ) λ=0 → ∞ does not exist as β ∈ (0, 1) (occurrence of extremely long waiting times). For β = 1 (12) retrieves the LT of the exponential density. Expanding the LT for λ smallK β (λ) ∼ 1 − λ β ξ , Laplace inversion shows that ML-PDF has a fat power law tail corresponding to a long memory. For the ML-case we have for (10) where λ −β corresponds to the Riemann-Liouville fractional integral of degree β andM β (λ) = ξλ 1−β the Riemann -Liouville (R-L) fractional derivative of order 1 − β. Denoting D ν t the R-L fractional derivative where it is sufficient here to consider the range µ ∈ (0, 1) defined by [28] having LT λ µf (λ) and with the definition of the Caputo fractional derivative Consult the Appendix and [29] for the connections with general (fractional) derivatives. In the limit µ → 1− both Caputo and RL-fractional derivatives converge to the first order standard derivative where we can then write (11) when all waiting times t C,I,R are drawn from ML PDFs in terms of Caputo fractional derivatives as retrieving for β = 1 the Eqs. (8) of exponential waiting time PDFs (see (15), (16)). Indeed the time-fractional case with fat-tailed waiting time distributions is of utmost importance deserving further thorough investigation, see Appendix B for a brief account.

Endemic equilibrium for waiting time PDFs with existing mean
Here we confine ourselves to waiting time PDFs with existing means. The endemic state is defined as the long-time limit of the evolution. That is we seek a stationary (constant solution) [s(t), c(t), j(t), r(t)] → [S e , C e , J e , R e ] and A(t) → A e = βJ e S e which also is the asymptotic solution for t → ∞ if the endemic state is stable. This solution can be obtained from Eqs. (9) in the limit of small λ, where these equations then take (withŝ(λ) → S e /λ,..) the form which has to be considered for λ → 0 where the lowest order in λ determines the endemic equilibrium. Now witĥ where with .. we denote mean values and o(λ) stands for the Landau symbol. Thus we havê and as said we assume that the means exist. Therefore, Then we get straight-forwardly for the lowest orders in λ the endemic equilibrium as follows and with our assumption A e = βJ e S e we arrive at Consider now the third relation in (22) which is an implicit equation for J e leading to where we have introduced The quantity R 0 = β t I can be interpreted as the 'basic reproduction number' and is a crucial (control-) parameter. (23) has the roots Clearly the endemic equilibrium exists only for those roots with J e ∈ [0, 1] where the full endemic state then is determined by (22). For b = 0 (j 0 = 0) one root is zero corresponding to the healthy state S e = 1. The second root is J e (j 0 = 0) = 2a = R 0 −1 β T ∈ [0, 1) for R 0 = β t I > 1 which therefore is the condition of the existence of an endemic state for j 0 = 0. In Appendix B we show for waiting time PDFs with existing means that R 0 > 1 also is the condition that the healthy state is unstable thus an epidemic can start to spread. We also show there that the healthy state is always unstable if K I (t) has a fat tail (i.e. infinite mean). Indeed for Consider now j 0 > 0 and R 0 > 1. Then we have a < 1 as t R +t C j 0 2T < 1 2 and R 0 −1 2β T < 1 2 and b < 1 thus only the positive root is whereas the other (J e ) 2 = a − √ a 2 + b < 0. The endemic state exists if (J e ) 1 ≤ 1. As a, b are monotonously increasing functions of j 0 let us check the root for the maximum value j 0 = 1 where a(j 0 = 1) = 1 2 − 1 2βT = 1 2 (1 − ǫ) and b(j 0 = 1) = 1 β T = ǫ. Then we get i.e. the initial and endemic state coincide j 0 = J e (1) = 1 where we infer that j = 1 (for R 0 > 1) should be a stable endemic state. From the monotony of J e (j 0 ) (see (24)) we can see that i.e. in the entire range j 0 ∈ [0, 1] the endemic value is given by the positive root (26) together with (24). This behavior is shown in the plot of Figure 1 where we draw J e versus j 0 for different values of R 0 > 1. One can see that J e is monotonously increasing with R 0 when t R , t C , j 0 are kept constant.
We focus now on initial conditions s 0 = 1− and j 0 = 0+ (healthy state). The endemic equilibrium (22) then writes and exists solely for R 0 > 1 depending only from R 0 and the mean waiting times t C,I,R . We may consider the following limiting cases. 1) t R → ∞ (infinitely long immunity), t C , t I , β are kept constant: Then we have t R T → 1 and t C T , t I T → 0 and hence S e = 1 R 0 remains unchanged thus C e , J e → 0 and R e → R 0 −1 R 0 with S e + R e = 1 which corresponds to a fully healthy population (susceptible or immune) with C e + J e → 0, i.e. the disease in the long time limit dies out. 2) In the same way the limits t I → ∞ (long illness time) and and long incubation time t C → ∞, respectively, are straightforward, where the endemic values of respective compartments with infinite waiting times tend to R 0 −1 R 0 with unchanged S e = 1 R 0 . This corresponds to the fact that the individuals remain eventually trapped in the respective compartments (C or I) with infinite waiting times thus the cyclic transition pathway S → . . . → S is suppressed. We infer that in the limiting cases 1), 2) the endemic states are stable similar to the class of classical SIR models and do not exhibit oscillatory instabilities.

δ-distributed waiting times
An instructive case consists in the deterministic limit when the sojourn times in the compartments for all walkers are constant. Then the waiting time PDFs are Dirac δ-functions with LTsK C,I,R (λ) = e −λt C,I,R . Recall now the convolution of two δ-kernels which yields a new δ-kernel with a shifted peak at t 1 + t 2 . With this observation we have Note that A(t) = A[j(t), s(t)] is causal, i.e. all functions with negative time arguments are vanishing. For t < t C the infected individuals accumulate in the compartment C and for t ≥ t C some start to leave compartment C to enter I with the delayed rate A(t − t C ). Then for t ≥ t C + t I the first individuals are healed starting transitions I → R with rate A(t − t C − t I ).
Finally for t ≥ t C − t I − t R the transitions R → S due to individuals losing their immunity occur. From these observations we can infer that where Θ(τ ) denotes the Heaviside unit step function defined as Θ(τ ) = 1 for τ ≥ 0 and Θ(τ ) = 0 for τ < 0. We can see for sharp waiting times the Heaviside functions 'switch on and off' the respective transitions between compartments. It is worthy of mention that (31) is consistent with (7) when we take into account the survival probabilities i.e. Φ C,I,R (t) = 1 for t < t C,I,R ('survival' in C, I, R, respectively) and Φ C,I,R (t) = 0 for t > t C,I,R ('death', having left C, I, R, respectively). From Eqs. (31) we find the endemic equilibrium S e , C e , J e , R e representing a stationary solution of (30). Plugging the stationary value A e = βJ e S e into (31) we re-arrive at Eqs. (22) with t C,I,R = t C,I,R for Dirac δ-distributions.

Stability analysis of endemic and healthy states
Here we investigate the stability of the endemic equilibrium for healthy initial conditions s 0 = 1 (j 0 = 0) for δ-distributed waiting times. To this end we set where u, v, w, x are 'small' time independent constants. Clearly at least one of these equations is redundant as s + c + j + r = 1. Then it follows that A(t) ≈ A e + A 0 (u, w)e µt where we take into account only the linear orders in u, v, w, x. Therefore, A e = βJ e S e and A 0 (u, w) = β(uJ e + wS e ). Plugging this into (30) leads to the following system of equations µ which are uniquely determined by u and w. Therefore the solvability condition is determined uniquely by the first and third equation for s and j containing only the coefficients u and w. The complete determinant of the system (34) leads to We have µ = 0 as a threefold eigenvalue and with one non-zero eigenvalue µ 1 determined by and hence µ 1 = µ 1 (β, j 0 , t C , t I , t R ). However, we will see that we can reduce the set of pertinent parameters. For the outbreak of an epidemic, it is necessary that the healthy (initial) state s 0 = 1 becomes unstable. To explore this issue we consider s(t) = s 0 + ue µt and j(t) = j 0 + we µt with s 0 = 1, j 0 = 0. Then we get the solvability condition by replacing in (35) J e → j 0 = 0 and S e → s 0 = 1 in the formμ The healthy state becomes unstable if there is a positive solutionμ 0 > 0 of this equation. Clearly there is a positive root only if d dμ g(μ, R 0 )|μ =0 = R 0 > 1, otherwise the healthy state is stable. This justifies our interpretation of R 0 as the basic reproduction number where the healthy state is unstable for R 0 > 1 thus the epidemic starts to spread. In Appendix B we show that the condition for an outbreak R 0 = β t I > 1 holds true for waiting time PDFs with existing means. This can be seen more closely as the right-hand side g(μ, R 0 ) of (37) is a concave function ofμ (see Figure 2). We have forμ small the expansion g(μ, R 0 ) = R 0μ + O(μ 2 ) >μ only if R 0 > 1 and always g(μ, R 0 ) → 0 <μ as µ → ∞. Therefore exists a positive root µ 0 (R 0 ) > 0 only for R 0 > 1. This behavior is sown in Figure 2 where µ 0 (R 0 ) increases monotonically with R 0 enhancing the (non-persistent) exponential growth of j(t) at the outbreak of the epidemic. Returning to the stability analysis of the endemic state we re-scale (36) with (38) as At threshold of an oscillatory instability, the eigenvalue is purely imaginary. Pluggingμ = iω into (39) and separating real and imaginary parts yields the two conditions which need to be both simultaneously fulfilled for an oscillatory (Hopf) instability where we introduced the reduced control parameter r 0 = J e R 0 = (R 0 − 1)/t 2 . In Figure 3 (left frame) we plot the zero lines of f i in the r 0 , ω plane for t 1 = 1/2, t 2 = 10. They intersect at r 0 ≈ 0.063, ω = 0.35, corresponding to a critical R 0 ≈ 1.63. The right frame of Figure 3 shows numerical solutions of the full delay system (30) with the same parameters t 1 , t 2 . The system is solved using an Euler forward scheme with fixed time step δt = 10 −4 . The initial conditions are j 0 = 10 −4 , s 0 = 1 − j 0 , r 0 = c 0 = 0 close to the healthy state. We varied R 0 from 1.5 (subcritical, upper left) to 2.0 (lower right frame). Clearly for R 0 = 1.7 the stable focus turns into a limit cycle that becomes wider with increasing R 0 .

Simple random walk model
To explore the SCIRS phenomenology we combine this model with a multiple random walker's approach which we implemented into a Python code [30]. In this random walk Z walkers navigate independently on a periodic two dimensional lattice. Each walker performs at integer times t = 0, 1, 2, . . . instantaneous independent random steps to a closest neighbour lattice point (simple walk). The position of walker j (j = 1, . . . Z) can be described by the random variables with the random steps This simple multi-walkers motion is a microscopic model with scaling limits to (standard) Brownian motion [27]. We assume a N x = N y = N -periodic lattice with x (j) (t) = x (j) (t) mod N , y (j) (t) = y (j) (t) mod N for the position of each walker j. In order to connect the random walk with the epidemic dynamics we apply the following infection rule [11]: If a walker j ∈ I meets a walker k ∈ S on the same lattice point ('collision' of an infected-infectious and a susceptible walker), i.e.
then the susceptible walker k gets infected with probability P inf and undergoes an instantaneous transition S → C followed by the transition pathway as described previously. The infection probability P inf is constant for all walkers and time independent. The random paths of the walkers are not affected by transitions between the compartments or 'collisions' of walkers.
In the simulations at each integer time instant t we count the populations Z S,C,I,R (t) in the compartments where the total population Z = Z S (t) + Z C (t) + Z I (t) + Z R (t) remains constant over time. As in the macroscopic model we focus on the fractions We implement the mutually independent waiting times t C,I,R as random numbers drawn from specific Gamma PDFs K C,I,R (t). The Gamma-distribution provides sufficient flexibility to generate a wide range of possible behaviors such as sharp δ-peaks, broadly scattered waiting times with a maximum or a monotonically decreasing PDF. The Gamma (also called Erlang-) PDF writes where α is the so-called shape parameter and 1/ξ the time scale parameter. The Gamma PDF has a maximum for α > 1 at t m = α−1 ξ and is monotonically decreasing for α ≤ 1 and weakly singular at t = 0 for α < 1. For α = 1 we get the exponential PDF K ξ,1 (τ ) = ξe −ξτ . Some cases of Gamma PDFs are drawn in histograms of Figure 4. Useful is its Laplace transform from which we can easily retrievê = α ξ (mean waiting time), and V = t 2 − t 2 = α ξ 2 (variance). We use in the simulations to generate constant (sharp) waiting times the feature, e.g. [11] lim where the mean τ 0 = τ is kept constant in this limit and the variance is vanishing. This limit is easily seen in the Laplace space

Validation of the macroscopic SCIRS model and case study
In the simulations we remove unimportant fluctuations by recording the ensemble averaged 3 compartment populations s(t) , c(t) , r(t) , j(t) . We average numerically over a number of equivalent random walk realizations with identical parameters, waiting time distributions and observation times. Each realization employs different random numbers ('python seeds') for the waiting times and random walk. We perform a case study to validate the macroscopic t to the endemic values where in the left frame the time is not very large. In the right frame the time is increased thus all r C,I,R are very close to one indicating excellent accordance with Eq. (28). The slightly different absolute values are due to slightly different R 0 = 1/S e despite all parameters P inf , Z, N are identical. This confirms our observation in Figures 7 that R 0 depends on further characteristics of the waiting time PDFs such as for instance their variances.

Conclusions
We studied a macroscopic four compartment SCIRS model with memory effects introduced by random compartmental waiting times. We derived evolution equations for different arbitrary waiting time distributions and considered pertinent cases such as exponential (memoryless), Mittag Leffler (fat-tailed with long memory), and sharp (δ-distributed) waiting times. We highlighted connections with general fractional calculus and showed that the evolution equations are of general (fractional) type (Eqs. (11), (49)). For waiting time PDFs with existing mean we obtained exact formulas for the endemic equilibrium and conditions of its existence and identified the 'basic reproduction factor' R 0 = β t I as crucial parameter controlling whether or not a SCIRS epidemic starts to spread. We also found (Appendix B) that the healthy state always is unstable if t I has a fat-tailed PDF regardless of the distributions of t C,R . In this case R 0 and t I do not exist. We interpret this instability by the occurrence of very long infectious periods t I strongly boosting the epidemic spreading. Due to their general importance, the time-fractional cases call for further thorough investigations.
We compared the macroscopic SCIRS model with a random walk approach. Our computer experiments have given numerical evidence that Eqs. (28) are fulfilled for any of the implemented waiting time distributions (with initial conditions very close to the healthy state) whenever the endemic equilibrium exists. This is further confirmed by comparison of results produced by the random walk simulations and by direct numerical integration of Eqs. (30) for two cases of δ-distributed waiting times, see Figures 6 (left), 8 and 10 (left), 11, respectively. All these results give strong evidence that random walks offer an appropriate microscopic picture of the SCIRS dynamics. It is also found in the simulations that waiting time PDFs with different variances but otherwise identical parameters may influence R 0 .
For future research it would be desirable to relate the endemic value S e (i.e. R 0 ) (here determined numerically) with characteristics of the waiting time PDFs (such as their variances and others), infection probabilities in a collision of I and S walkers and of the random walk. A further challenge is the development of a microscopic random walk model for the oscillatory (Hopf) instabilities (condition for persistent oscillatory behavior).
Generally SCIRS type models with memory as introduced in the present paper may open a wide research field to better understand the phenomenology of real world epidemic dynamics.

Acknowledgement
T.G. gratefully acknowledges to have been hosted at the Institut Jean le Rond d'Alembert for the sake of the present study and development of the Python codes.

A General derivatives
We can alternatively represent Eqs. (11) in the Riemann-Liouville sense of general derivatives as which also is obtained straight-forwardly from the Laplace transformed representation (9) and employing definition (10). On the right-hand sides general derivatives come into play. The notion of 'general' (fractional) derivative was introduced by Kochubei [29] generalizing d dt y(t) as (in the Caputo sense) . The part D t · y(t) = t 0 k(t − τ )y(τ )dτ is a general derivative in the Riemann-Liouville sense for some admissible kernels k(τ ) (see [29] for an outline of this theory). The general derivative (50) has the LT The general derivatives contain the class of standard fractional derivatives of Caputo and Riemann-Liouville type and for a k(t) = δ(t) the standard first order derivative.

B Stability analysis for arbitrary waiting time densities
In Section 3 we have derived that the healthy state is unstable if R 0 = βt I > 1 for δ-distributed kernels. Here we consider arbitrary waiting time kernels first with finite means to generalize this result and at the end of this appendix we briefly look at the time fractional case. To this end we take into account that the PDFs K C,I,R (t) can be seen as a superposition (average over t C,I,R ) of the δ(t − t C,I,R )-kernel as follows and in this way we can average f (t C,I,R ) = ∞ 0 K C,I,R (τ )f (τ )dτ (53) for sufficiently good functions f thus e −µt C,I,R =K C,I,R (µ) (54) yielding the Laplace transform of the PDF. Therefore, averaging Eqs. (30) for δ-kernels over the waiting times t C,I,R brings us back to the general SCIRS Eqs. (1) which can be seen from where we used causality of A and K C . Hence we can generalize the case of δ-kernels (34) to any kernels by accounting for e −µ(t I +t R +t C ) =K C (µ)K I (µ)K R (µ) where we always use the mutually independence of the waiting times. Hence by averaging Eq.
(34) we get A 0 (u, w)[1 −K C (µ)K I (µ)K R (µ))] + µu = 0 We have again (as in (35)) µ = 0 as a threefold eigenvalue. The remaining non-zero eigenvalue is determined by We have to point out that in order to obtain these results we have relaxed causality in the exponential ansatz (33) and assumed that it is defined for all t ∈ R. The difference of causal and non-causal averaging can be seen by comparing where Θ(τ ) indicates the Heaviside unit step function and corresponds to causality. To obtain the condition (57) we have averaged as in the first line of (59) where the exponential is non-causal allowing exponentials with arguments t−t C < 0 to contribute. Both sides become asymptotically equal in the limit t → ∞. Hence above relation (57) can also be interpreted as the large time asymptotics of causal averaging. Let us focus now on the stability of the healthy state s 0 = 1. Then (58) yields (S e → s 0 = 1, As in the case of δ-distributed waiting times G(µ) is a concave function with G(µ) → 0 as µ → ∞ thus we can argue in the same way (see Eq. (37) and Figure 2). Hence there is an intersection µ 0 > 0 only if d dµ G(µ) µ=0 > 1. Expanding G(µ) (assuming existing mean t I ) gives Therefore, d dµ G(µ) µ=0 =: R 0 = β t I > 1 (62) is the condition for this instability of the healthy state and hence the interpretation as 'basic reproduction factor' R 0 makes sense. This holds for waiting time PDFs with existing means. On the other hand relation (60) holds for any waiting time PDF including fat-tailed ones which we briefly consider in the following.

Time fractional case
Assuming K I (t) is a fat-tailed PDF such as Mittag-Leffler withK I (µ) = ξ I ξ I +µ ν , ν ∈ (0, 1) and no matter whether or not the other kernels have existing means we get for G(µ) where o(µ ν ) indicates orders higher than ν. Thus (60) (independent of β ξ I ) always has a positive intersection since ν ∈ (0, 1). We conclude that for fat-tailed K I (t) the healthy state always is unstable regardless of the distributions of t C and t R . Physically this instability can be understood by the occurrence of very long infectious times t I strongly boosting the epidemic spreading (corresponding to the limit R 0 , t I → ∞).