We consider a generalized SEIS (susceptible, exposed, infectious, and susceptible) model where individuals are divided into three compartments: S (healthy and susceptible), E (infected but not just infectious, or exposed), and I (infectious). Finite waiting times in the compartments yield a system of delay-differential or memory equations and may exhibit oscillatory (Hopf) instabilities of the otherwise stationary endemic state, leading normally to regular oscillations in the form of an attractive limit cycle in the phase space spanned by the compartment rates. In the present paper, our aim is to demonstrate that in the dynamics of delayed SEIS models, persistent chaotic attractors can bifurcate from these limit cycles and become accessible if the nonlinear interaction terms fulfill certain basic requirements, which to our knowledge were not addressed in the literature so far. Computing the largest Lyapunov exponent, we show that chaotic behavior exists in a wide parameter range. Finally, we discuss a more general system and show that a sudden falloff of the infection rate with respect to increasing infection number may be responsible for the emergence of chaotic time evolution. Such a falloff can describe mitigation measures, such as wearing masques, individual isolation, or vaccination. The model may have a wide range of interdisciplinary applications beyond epidemic spreading, for instance, in the kinetics of certain chemical reactions.

Within the frame of population dynamics, compartmental low-dimensional schemes are widely discussed and applied in epidemiological models, such as SIR (Susceptible–Infected–Recovered) or its various extensions. In their standard forms, these models show instabilities and relaxational behavior in time with stable endemic equilibrium fixed points. More generalized variants include time delay terms (delay-differential equations) and may develop solutions in form of limit cycles, exhibiting periodic behavior of the infection rates. We demonstrate that a minimal compartment model, including a delayed infection, may show even chaotic trajectories in phase space if the infection mechanism fulfills certain simple conditions. A comparison with the dynamics of the Mackey–Glass system is presented as well.

Beginning with the book of Thomas Robert Malthus in 1798,1 the modeling of population dynamics has a long tradition. The nowadays called “Malthusian growth model” was later generalized by Richards2 and others. Predator–prey systems were studied from 1925 by Lotka and Volterra.3,4 The seminal work of Kermack and McKendrick5 in 1927 where the nowadays called “SIR model” was introduced can be seen as the initiation of epidemic models. Being a specialization of a predator–prey system, the original SIR equations have obtained many extensions and modifications. A huge number of publications exists based on this model, which is built on a low-dimensional system of ordinary differential equations; for a review, see Ref. 6.

SIR models refer to the class of “compartmental models” since they divide the individuals into several compartments depending on their state of health. In this way, SIR stands for an acronym from ( S = susceptible, I = infected, R = recovered). It turned out that the features of infectious diseases, such as measles, mumps, and rubella, could to a certain extent be captured by such simple models.

On the other hand, all these diseases come in waves with a more or less regular periodicity. The waves can be triggered by external reasons, such as seasons, but also by some internal mechanisms. Such an intrinsic behavior was already speculated in 1929 by Soper in a model for the time evolution of measles cases.7 

In the standard SIR models, the interplay between infected and susceptible individuals is inspired by an interaction of the predator–prey kind where the infection rate is expressed in the form of a simple bilinear term β 0 j ( t ) s ( t ). Here, j ( t ) and s ( t ) are the relative number of infected ( I) and susceptible ( S) individuals and β 0 is related to the constant probability of infection at each contact (infection rate). The predator (infected) “catches” the prey (susceptible) by infection. Predator–prey models of the classical SIR type are not able to describe sustained oscillations that originate from the instability of a fixed point, here the endemic equilibrium, but rather account for single outbreaks if the healthy state becomes unstable. In the long time limit, herd immunity and the endemic equilibrium as a stable fixed point is reached where the fractions of the population in the different compartments attain constant values. However, it turns out that there is a further class of predator–prey models with modified infection rates A ( s ( t ) , j ( t ) ) and delayed transitions exhibiting oscillatory instabilities and routes to chaos.

Other work8 considered a nonlinear infection rate according to
(1)
and obtain limit cycle solutions for certain parameters m = n 2 , α > 0. Tang et al.9 studied the case m = n = 2 and found a weak focus and the existence of two limit cycles. For m = n = 1, the dynamics is qualitatively the same than for the standard SIR model, and sustained oscillations cannot be observed. The denominator 1 + α j n accounts for mitigation measures against the epidemics that naturally increase with increasing j. Note that for the case n > m, the interaction A has a maximum at a certain infection number. For such non-monotonic behavior, it was shown in Ref. 10 that the dynamics in the long time limit approaches a stable fixed point as for the original SIR models.

In earlier work, we showed that a finite (long) immunity lifetime11 as well as delayed mitigation measures12 lead to a system of delay-differential or memory equations and may show oscillatory instabilities of the otherwise endemic state, leading to normally regular oscillations of the compartment rates. Memory terms were introduced in epidemiological models by many other researchers; for an overview, see Ref. 13. From the mathematical point of view, the presence of a delay term in an ordinary differential equation makes a low-dimensional system infinitely dimensional and may allow for the occurrence of periodic, quasiperiodic, or even chaotic behavior, rendering the dynamics much more complex.14–16 

A crucial element of our extended model studied in the present paper is the introduction of a class of infection rate functions exhibiting a sudden falloff with respect to the infection numbers, according to (1) for large enough n. We demonstrate that the complex interplay of these kinds of infection rates and delayed transitions between the compartments is the actual source of the chaotic dynamics. The well known Mackey–Glass equation17 refers to our model class, and we show in Sec. IV that the Mackey–Glass equation can be derived systematically from our model for the approximation of a constant susceptibility rate. However, our model generalizes the Mackey–Glass system in two respects. First, it highlights the interpretation as a compartmental model and considers continuously distributed compartmental waiting time distributions for the delayed transitions, and second, it admits a wide range of infection rate functions (beyond the one used in the Mackey–Glass model) with the above-mentioned falloff feature. To the best of our knowledge, there is so far no such approach in epidemic modeling.

Our aim is to develop a model as close as possible to the non-delayed standard models. For further sake of simplicity, we consider first only two compartments, namely, S and I, and assume that recovered individuals transfer directly back to the S compartment, leading to a SIS model. The crucial point is that the infection S to I is assumed to take a certain finite time. This can be modeled by introducing a new compartment E (exposed) where the individuals are already infected but not infectious. The sojourn time in this compartment corresponds to the latent time of the infection. The transition from E to I is defined as the moment where the individual becomes infectious. Symptoms may emerge even later. The time between infection and occurrence of symptoms is typically named incubation time. In general, the incubation time is lager than the latent time. Note that some authors do not distinguish between latent and incubation times. Throughout this paper, however, we shall use the above definition.

In the literature, the compartment containing latent individuals is also named “exposed.” Thus, we shall consider here a SEIS (susceptible, exposed, infectious, and susceptible) model with finite delay times for the transitions E I and I S, respectively.

Our model could apply to diseases with a pronounced latent period. The latent times can be distributed stochastically and drawn from a probability density function that should have a single maximum at a finite time (unimodal). Infectious diseases, such as HIV/AIDS [infectious incubation (I) up to years and non-infectious latent period (E) some weeks], tuberculosis [latent tuberculosis (E) some years], but also common colds or influenza where both latent and incubation times are on a much shorter time scale, come to mind.

This paper is organized as follows: In Sec. II, we introduce the SEIS model for general distributions of the waiting times in the different compartments. Equilibrium points and their stability are given and examined in more detail for δ-shaped kernel functions, leading to a set of two coupled delay-differential equations for the fractions of S and I. The critical conditions for a Hopf bifurcation are determined analytically. In Sec. III, numerical solutions of a delay system are computed and the emergence of a chaotic attractor is demonstrated; its largest Lyapunov exponent is found to be positive for large regions of the basic reproduction number. In Sec. IV, the analogy to the Mackey–Glass system is shown for the approximation of a constant susceptible population. A mechanism for the occurrence of chaotic attractors is proposed.

In the present paper, we shall consider a generalized SEIS model without birth or death processes, leading to a constant total population of N individuals. According to their state of health, the individuals are divided into the three compartments S (healthy and susceptible), E (exposed, infected but not just infectious), and I (infectious), see Fig. 1, with the total population number N = S + E + I.

FIG. 1.

The SEIS model with a j-dependent basic reproduction number R ( j ) and finite (random) delay times τ E , τ I between compartments E, I and I, S, respectively.

FIG. 1.

The SEIS model with a j-dependent basic reproduction number R ( j ) and finite (random) delay times τ E , τ I between compartments E, I and I, S, respectively.

Close modal

Contrary to the standard SEIRS models, we neglect the compartment R, which means that we consider diseases with no or very short immunity phase. This allows the individuals to become infected again after the duration of illness. However, the crucial difference of our model is the finite time for the transition from S to I, biologically motivated by a rather long latent time, where the individuals stay in the compartment E. This transition time is chosen from a given probability distribution having a maximum at a finite time. The special case of an exponentially and monotonically decaying probability distribution would lead to the standard SEIS without memory.

We shall show below that this finite latent time together with a special form of the nonlinear infection rate as already outlined in Eq. (1) is responsible for the occurrence of limit cycles, period doublings, and, finally, the emergence of chaotic attractors. All these new features are not found in the standard SEIS model that simply shows the transition from an unstable fixed point (healthy state) to a stable one (endemic equilibrium) if the basic reproduction number exceeds one.

We define the fractions
with s + e + j = 1. The model reads
(2a)
(2b)
(2c)
where A ( t ) indicates the infection rate, which describes the nonlinear interaction at time t. The delay time τ E between the instant of infection (transition S E) and transition E I has the interpretation of a latent time. The latent time is followed by the sojourn time τ I in compartment I where the individual is infectious, meaning either showing symptoms, or, during the incubation period, with no symptoms. After spending the time τ I in compartment I, the individual recovers, undergoing a transition I S.
Considering the transitions occurring at time instant t gives an evocative interpretation of Eq. (2): A ( t ) indicates the infections (transitions S E) at instant t. The term A ( t τ E ) accounts for the past infections at t τ E for which at instant t, the latent time has elapsed and which undergoes a transition E I at instant t; see (2b) and (2c). Furthermore, the term A ( t τ E τ I ) comes from the infections that happened in the past at instant t τ E τ I for which at time t both, the latent time τ E and the illness time τ I have elapsed and, hence, are recovering at instant t with transition I S; see (2c) and (2a). In Eq. (2), stands for the means with respect to the independent random times τ E and τ I drawn from the normalized PDFs (memory functions) K E , K I and are computed from
We note that the system (2) is included in the more general case of a SEIRS cycle, examined in detail in a recent paper.16 
For the interaction A that models the infection rate, we consider the class
(3)
with a monotonically decreasing function,
(4)
with some positive α and integer n > 1 according to (1). This leads to a decrease of the effective basic reproduction number if j is increasing and vice versa, modeling mitigation measures, such as containment, masques, etc. Note that A has a maximum at
The parameters α and n can be used to adjust the infection rate to certain mitigation scenarios. Whereas α merely scales j so that for larger α, the measures become effective already for smaller infection numbers, and the parameter n turns out to be crucial for our following analysis. If n 1, the measures become effective rather abruptly if j exceeds 1 / α in the form of a pronounced sudden falloff of the infection rate; see Fig. 2.
FIG. 2.

The function f ( j ) for α = 100 and three different values of n. For larger n, a stepwise behavior can be seen that is crucial for the occurrence of a chaotic dynamics.

FIG. 2.

The function f ( j ) for α = 100 and three different values of n. For larger n, a stepwise behavior can be seen that is crucial for the occurrence of a chaotic dynamics.

Close modal
Instead of studying the full system (2), we assume an exponential distribution for K I according to
(5)
which has the mean value τ I = 1 / γ. It can be easily seen by differentiation that for such a PDF, the relation
holds. Since e is not occurring explicitly, it is sufficient to consider Eqs. (2a) and (2c). Scaling time with t = t / γ and substituting (5) with j ( 0 ) = 0 yields the two coupled memory equations,
(6a)
(6b)
with R 0 = R 0 / γ , τ = τ γ. In the following, we shall leave all primes and time is measured in units of the time of recovery. The standard SIS model is recovered for α = 0 and K E ( τ ) = δ ( τ ). For K E ( τ ) = δ ( τ τ E ) with Dirac’s δ-function δ, one has a sharp delay time τ E for the latent period. This is a special case that will be treated in Sec. II C in more detail.
There are two equilibrium solutions or fixed points of (6). The first one is the state before the epidemic breaks out (healthy state) and reads
(7)
The other one is the endemic equilibrium found from
(8)
or
(9)
where we used j e q + s e q + e e q = 1 and e e q = A τ E .
We compute the stability of both equilibrium points applying linear stability analysis. Inserting
where s 0 , j 0 stand for the equilibrium point in (6) yielding a linear system that has the transcendental equation
(10)
as a solvability condition. Here,
(11a)
(11b)
and K ^ E as a Laplace transform of K E,
(12)
Note that R 0 is retrieved from R 0 = A j | s = 1 , j = 0. The state (7) yields a = 0 , b = R 0, and from (10),
The critical point is λ = 0. Expanding (12) yields
and in a linear order
Thus, the healthy state is stable for R 0 < 1 and becomes an unstable saddle node for R 0 > 1.
For the special case of a δ-shaped kernel, one has
and (6) turns into a system of delay-differential equations,
(13a)
(13b)
For τ E = 0 , K ^ E = 1, and the endemic equilibrium has the two eigenvalues
The first one belongs to a marginal (Goldstone) mode that describes a shift between j e q and s e q since for τ E = 0, the sum of the two equations (13) yields d t ( s + j ) = 0. The second eigenvalue is always negative, and the stability of the endemic state is proven without restrictions.
Since τ E > 0, a Hopf instability becomes possible and (10) can be solved with λ = i ω, leading to
(14)
and
(15)
Thus, the endemic equilibrium becomes unstable if τ E exceeds τ H with the critical frequency (14). Equation (15) has to be solved iteratively because a , b depend on τ. Figure 3 shows τ H, and the period T = 2 π / ω over R 0.
FIG. 3.

Critical delay time necessary for an oscillatory instability of the endemic equilibrium and the period of Hopf frequency (red).

FIG. 3.

Critical delay time necessary for an oscillatory instability of the endemic equilibrium and the period of Hopf frequency (red).

Close modal

We solved the nonlinear delay system by a standard fourth order Runge–Kutta method with fixed time step δ t = 10 4; for details, see Ref. 18.

For large values of n, j e q becomes an unstable focus if τ E exceeds a critical value. For even larger delay times, a chaotic attractor emerges. Figure 4 shows a time series for the periodic case with
and for the chaotic case with τ E = 2.
FIG. 4.

j ( t ) (black) and s ( t ) (red) in the periodic (left frame) and chaotic (right frame) regime.

FIG. 4.

j ( t ) (black) and s ( t ) (red) in the periodic (left frame) and chaotic (right frame) regime.

Close modal

Figure 5, left frame, shows the chaotic attractor at τ E = 2 where j ( t ) is plotted vs j ( t τ E ). In Fig. 5, right frame, the same attractor is plotted in the s j phase plane.

Contrary to the standard SEIS model, it is remarkable that the variation of s is very small. This agrees well with the real world data, for instance, for COVID, where herd immunity was always far from being reached during or between the wave outbursts; see also Ref. 19.

FIG. 5.

Chaotic attractor for τ E = 2, left: j ( t ) over j ( t τ E ), right: j ( t ) over s ( t ).

FIG. 5.

Chaotic attractor for τ E = 2, left: j ( t ) over j ( t τ E ), right: j ( t ) over s ( t ).

Close modal
To demonstrate that the attractors found are really chaotic, we compute the largest Lyapunov exponent for two different values of τ E. First, a main trajectory s 0 ( t ) , j 0 ( t ) is computed iterating (6) by a Runge–Kutta method, then the system is linearized with respect to the main trajectory,
From t 0 = 4000, the linear delay system
(16a)
(16b)
is solved numerically by the same Runge–Kutta scheme. The chaotic time dependent functions a , b are obtained from (11) by replacing j 0 , s 0 by s 0 ( t ) , j 0 ( t ). Starting with a normalized random initial condition,
where ξ t η t are independent equally distributed random variables in [ 0.5,0.5], the L 2 norm
is computed each time after t reaches k τ E and ( u , v ) is normalized again; see, e.g., Ref. 18. This process is repeated until t e = 8000, and the largest Lyapunov exponent is found according to
The result is shown in Fig. 6.
FIG. 6.

Largest Lyapunov exponents for τ E = 2 (blue) and τ E = 3 (red) over R 0. Periodic windows can be clearly seen.

FIG. 6.

Largest Lyapunov exponents for τ E = 2 (blue) and τ E = 3 (red) over R 0. Periodic windows can be clearly seen.

Close modal
Next, we study the system (6) with a kernel function of a certain finite width. To this end, we take the Erlang function defined as
(17)
where Γ ( κ ) denotes the Γ-function and Γ ( n ) = ( n 1 ) ! for integer n > 0. The Erlang function is located around its maximum τ m = ( κ 1 ) / ξ and has a width 1 / ξ, resulting in a mean sojourn time in compartment E of τ E = κ / ξ; see Fig. 7. Thus, the two parameters ξ and κ can be used to adjust the distribution (17) to certain empirical data. The Erlang function approaches the δ-function δ ( τ τ E ) for ξ , κ but κ / ξ finite. Figure 8 shows some numerical solutions for different ξ and κ = 2 ξ, resulting in a fixed τ E = 2. Chaotic attractors are only obtained if ξ exceeds a certain value that depends also on τ E and R 0. This shows that the width of the distribution for the sojourn time in E must not exceed a certain critical value to obtain a chaotic evolution in phase space.
FIG. 7.

Erlang functions (17) for τ E = 2 and different ξ.

FIG. 7.

Erlang functions (17) for τ E = 2 and different ξ.

Close modal
FIG. 8.

Phase plots in the s j-plane for different ξ and τ E = 2 , R 0 = 2.5.

FIG. 8.

Phase plots in the s j-plane for different ξ and τ E = 2 , R 0 = 2.5.

Close modal
To demonstrate this more clearly, we computed the largest Lyapunov exponent also for this case. Now, Eq. (16) turns into a linear memory system,
(18a)
(18b)
The rest of the procedure remains the same as described in Sec. III B. The result for different ξ is shown in Fig. 9.
FIG. 9.

Largest Lyapunov exponents for τ E = 2 and the Erlang kernel (17) with ξ = 50 (gray), ξ = 40 (blue), and ξ = 30 (red).

FIG. 9.

Largest Lyapunov exponents for τ E = 2 and the Erlang kernel (17) with ξ = 50 (gray), ξ = 40 (blue), and ξ = 30 (red).

Close modal
Introducing the new variable x ( t ) = α j ( t ), Eq. (6a) turns into
To have an effective feedback, α j should be of order one, and since j is rather small, α 1. As a first approximation, we may then put
where s 0 is the initial value and, therefore, close to one. Then, from (6b), we find immediately
(19)
where time is scaled again with t = t ~ / R 0 s 0 and μ = 1 / R 0 s 0 < 1.
Equation (19) can be considered a Mackey–Glass equation with memory. It turns into the original Mackey–Glass equation for a δ-kernel K E = δ ( τ τ E ) and has the form (tildes removed)
(20)
It is well known since the remarkable paper of Mackey and Glass17 who derived it as a model for the control of blood cell concentration that this equation has a chaotic attractor for rather large n and if τ E exceeds a certain critical value. The attractor is born after an oscillatory instability of the non-trivial fixed point followed by period doubling. The bifurcation scenario and the resulting trajectories have much in common with our solutions shown in Fig. 5, left frame.
Instead of (20), let us examine the more general delay equation,
(21)
where f is a monotonically decreasing function with
and a certain more or less steep decrease close to x = 1. As a test function fulfilling all requirements, we may choose
(22)
where 1 / β denotes the width of the step and plays the role of the parameter n of the function f(j) defined in (4). Then, (21) has an unconditionally unstable fixed point at x = 0 and a second fixed point at x = 1 that is stable for τ < τ c and becomes an unstable focus at τ = τ c. A linear analysis yields
(23)
with the Hopf frequency ω c, Fig. 10. A limit cycle occurs if τ exceeds τ c, followed by several period doublings until eventually a chaotic attractor may emerge; see Fig. 11. Therefore, the steepness β of (22) plays a crucial role. The threshold (23) only exists if β > 2. Moreover, if f ( x ) is too flat, chaos never occurs or only for large delay times τ. The value x = 1 could be interpreted as a kind of separatrix. For x > 1, the delay term in (21) has no or only little effect and x decreases exponentially until the trajectory is back in the “delay region” x < 1. Then, the delay term determines again the dynamics and may throw x across the separatrix. This back and forth motion is sensible with respect to small perturbations and may be finally responsible for the chaotic nature of the trajectory.
FIG. 10.

Critical delay time necessary for an oscillatory instability according to (23) and the period of Hopf frequency (red). Note the similarity to Fig. 3.

FIG. 10.

Critical delay time necessary for an oscillatory instability according to (23) and the period of Hopf frequency (red). Note the similarity to Fig. 3.

Close modal
FIG. 11.

Numerical solutions of (21) for several values of β and τ, x ( t ) over x ( t τ). Chaos occurs earlier for larger values of β. For β < 2, the limit cycle ceases to exist, and the fixed point x = 1 remains stable for arbitrary τ.

FIG. 11.

Numerical solutions of (21) for several values of β and τ, x ( t ) over x ( t τ). Chaos occurs earlier for larger values of β. For β < 2, the limit cycle ceases to exist, and the fixed point x = 1 remains stable for arbitrary τ.

Close modal

Contrary to the standard SIR model and its various extensions, time-periodic trajectories occur in a natural way if the ordinary differential equations are not local in time but contain memory terms where the history of the evolution is included. We have shown that if these terms have a special form, chaotic solutions may exist that are not caused by fluctuations and, therefore, still deterministic. Starting from the most simple epidemic system containing only the two compartments for susceptible and infected individuals, we derived a minimal model by introducing a finite delay for the infection mechanism to become effective. Such a delay could be justified by an intermediate incubation state with finite waiting time, the latent time, where the individuals are infected but not yet infectious. We demonstrated that if the latent time, the recovery time, and the basic reproduction number have the appropriate relations, chaotic attractors emerge and persist around the endemic equilibrium.

The present model has a rather large potential of generalizations. One direction of interest may be the extension to multiple compartments in order to explore whether or under which conditions chaos still emerges. On the other hand, the exploration of infection rates exhibiting sudden falloffs with noisy slopes and other parameters may be of interest as well. In this way, the effect of natural fluctuations in the implementation of mitigation measures could be captured.

The authors have no conflicts to disclose.

Michael Bestehorn: Conceptualization (equal); Formal analysis (equal); Software (equal). Thomas M. Michelitsch: Conceptualization (equal); Formal analysis (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
T. R.
Malthus
, “An essay on the principle of population,” in Oxford World’s Classics (Oxford University Press, Oxford, 2008).
2.
F. J.
Richards
, “
A flexible growth function for empirical use
,”
J. Exp. Bot.
10
,
290
(
1959
).
3.
A. J.
Lotka
, “
Contribution to the theory of periodic reaction
,”
J. Phys. Chem.
14
,
271
(
1910
).
4.
V.
Volterra
, “
Variazioni e fluttuazioni del numero d’individui in specie animali conviventi
,”
Mem. Acad. Lincei Roma.
2
,
31
113
(
1926
).
5.
W. O.
Kermack
and
A. G.
McKendrick
, “
A contribution to the mathematical theory of epidemics
,”
Proc. R. Soc. A
115
,
700
(
1927
).
6.
R. M.
Anderson
and
R. M.
May
,
Infectious Diseases in Humans
(
Oxford University Press
,
Oxford
,
1992
).
7.
H. E.
Soper
, “
The interpretation of periodicity in disease prevalence
,”
J. R. Stat. Soc.
92
,
34
(
1929
).
8.
W.
Liu
,
A.
Levin
, and
Y.
Iwasa
, “
Influence of nonlinear incidence rate upon the behavior of SIRS epidemiological models
,”
J. Math. Biol.
23
,
187
(
1986
).
9.
Y.
Tang
,
D.
Huang
,
S.
Ruan
, and
W.
Zhang
, “
Coexistence of limit cycles and homoclinic loops in a SIRS model with a nonlinear incidence rate
,”
SIAM J. Appl. Math.
69
,
621
(
2008
).
10.
D.
Xiao
and
S.
Ruan
, “
Global analysis of an epidemic model with nonmonotonic incidence rate
,”
Math. Biosci.
208
,
419
(
2007
).
11.
M.
Bestehorn
,
T. M.
Michelitsch
,
B. A.
Collet
,
A. P.
Riascos
, and
A. F.
Nowakowski
, “
Simple model of epidemic dynamics with memory effects
,”
Phys. Rev. E
105
,
024205
(
2022
).
12.
M.
Bestehorn
and
T. M.
Michelitsch
, “
Oscillating behavior of a compartmental model with retarded noisy dynamic infection rate
,”
Int. J. Bifurc. Chaos Appl. Sci. Eng.
33
,
2350056
(
2023
).
13.
F. A.
Rihan
,
Delay Differential Equations and Applications to Biology
(
Springer Nature
,
Singapore
,
2021
).
14.
G. E.
Hutchinson
, “
Circular causal systems in ecology
,”
Ann. N.Y. Acad. Sci.
50
,
221
(
1948
).
15.
M.
Bestehorn
,
E. V.
Grigorieva
, and
S. A.
Kaschenko
, “
Spatio-temporal structures in a model with delay and diffusion
,”
Phys. Rev. E
70
,
026202
(
2004
).
16.
T.
Granger
,
T. M.
Michelitsch
,
M.
Bestehorn
,
A. P.
Riascos
, and
B. A.
Collet
, “
Four-compartment epidemic model with retarded transition rates
,”
Phys. Rev. E
107
,
044207
(
2023
).
17.
D.
Mackey
and
L.
Glass
, “
Oscillations and chaos in physiological control systems
,”
Science
197
,
28
(
1977
).
18.
M.
Bestehorn
,
Computational Physics
(
De Gruyter
,
Berlin
,
2023
).
19.
L.
Gostiaux
,
W. J. T.
Bos
, and
J.-P.
Bertoglio
, “
Periodic epidemic outbursts explained by local saturation of clusters
,”
Phys. Rev. E
107
,
L012201
(
2023
).