Interaction between delays and seasonality
Traditional epidemic models consider that individual processes occur at constant rates. That is, an infected individual has a constant probability per unit time of recovering from infection after contagion. This assumption certainly fails for almost all infectious diseases, in which the infection time usually follows a probability distribution more or less spread around a mean value. We show a general treatment for an SIRS model in which both the infected and the immune phases admit such a description. The general behavior of the system shows transitions between endemic and oscillating situations that could be relevant in many real scenarios. The interaction with the other main source of oscillations, seasonality, is also discussed.
This work was presented at the 9th AIMS Conference on Dynamical Systems, Differential Equations and Applications, Orlando, FL, USA, July 1 - 5, 2012.
Keywords: epidemics, delay equations
Many infectious agents are responsible of oscillating epidemics. In some cases these oscillations are almost periodic, appearing and disappearing with time in a well established way. Examples abound, even in classical textbooks [1, 2], of such measles, typhus or cholera epidemic waves. Some of these diseases have a cyclic natural history, in which a susceptible subject can become infected, then recovered and immune, and finally return to the susceptible state. However, the mere cyclic nature of the disease does not grant an oscillatory behavior of its epidemic, as can be exemplified by gonorrhea [3, 4]. On the other hand, many epidemics are strongly affected by seasonal factors, which could in principle produce a yearly variation in their prevalence. There are many instances, however, of epidemics that oscillate with longer periods: measless  and human parainfluenza  every two years, pertussis with a four year period , syphilis with 11 years , etc. Even for influenza-like illness, with its apparent yearly period, there exist some indications that its behavior has actually a period of two years .
How do these oscillations arise in a population, apparently synchronizing the infective state of many individuals? Indeed, it is known that external causes can produce oscillations in epidemic systems: seasonal or El Niño driving in Hantavirus, for example . Yet, inherently internal processes are also able to produce oscillations, such as a stochastic dynamics [8, 9], or a complex network of contacts . What are the roles and the interplay of external agents and the dynamic natural history of the disease?
The usual mean-field formulation of epidemic models assumes that the processes that change the states of the individuals proceed at constant rates. For an SIRS system we have three such processes: , susceptibles to infected through contagion, , infected to recovered, and , recovered to susceptible through loss of immunity. So, three rates characterize the dynamics. Let us call the probability of contagion per unit time in two-persons interactions. For the two other processes we can use characteristic times: is the duration of the infection, the infectious time (during which the agent stays in category , infected and infectious), and is the immune time, the duration of immunity in the state (see Fig. 1). The dynamics can be represented as follows:
where , and stand for the corresponding fractions of susceptible, infectious and recovered individuals in the population. Moreover, no demographic change has been included in so and the system is effectively two-dimensional.
As observed already by Anderson and May  this mathematically convenient treatment of the duration of the infection as a constant rate is rarely realistic. It is more common that recovery from infection takes place after some rather well defined time. A similar consideration can be made about the transition from recovered to susceptible. While the “constant rates” SIRS model shows, at most, damped oscillations towards an endemic state, it is known (see for example ) that sustained oscillations are possible in models with a fixed delay in the recovered to susceptible transition.
In a previous paper  we analized an SIRS model with distributed delays both in the infectious and the immune times. In our formulation the fixed times and the constant rates situations are recovered as extreme cases of the probability distribution of the delays. In other words, we can go from the SIRS model with deterministic delays to the classical constant rates SIRS, covering all the intermediate situations in between.
In the present contribution we analize the interplay of a delayed SIRS model with an external forcing that represents, for example, seasonal variations. We show that a reach dynamics arises, including non-seasonal periodicity of the epidemics of the kind observed in many real scenarios.
2 Delayed SIRS model
In  we developed a general SIRS system in which both the recovery and the loss of immunity processes are described by probability distribution functions of the corresponding times. Let be the probability (per unit time) of losing infectivity at time after having become infected at time . can be used as an integration kernel in a delayed equation for the infectious. The individuals that recover from infection at time are those that got infected at any time and then recover:
Correspondingly, let us represent the loss of immunity by a second kernel, . Now, the individuals that lose their immunity at time are those that got infected at earlier times, then recovered at some intermediate time with probability , and finally return to the susceptible state with probability :
which, given the constraint , are enough to describe the dynamics.
Appropriate choices of and in Eqs. (6) and (7) can reproduce different dynamical behaviors. For example, exponentially decaying distributions with characteristic times and respectively correspond to the constant rates model (local in time) represented by Eqs. (1)-(3). On the other hand, Dirac-delta distributions and represent a deterministic system in which each of the two processes, recovery and loss of immunity, last exactly and respectively:
where . Between these two extremes, other suitable and are able to represent more realistic situations, in which each process has a maximum at the corresponding characteristic time and some spread around it. One can use a convenient model for these in the form of gamma distributions adequately parametrized. They allow an analytic treatment of the stability of the fixed points. We used:
These distributions have mean and , respectively, for any value of the parameters and . Besides, they cover the range from exponentials (when or ) to Dirac delta distributions (when or ), with smooth bell-shaped functions for intermediate values of the parameters.
The model defined by Eqs. (8)-(9) must be supplemented with appropriate initial conditions. For an epidemic system, a reasonable choice is an outbreak of infection at the initial time: , , for example. Due to the nonlocality in time such conditions are insufficient, and extended initial conditions should be provided in order to have a well-posed problem. In more abstract contexts it is customary to provide arbitrary functions and in the interval . From an epidemiological point of view, however, it is more reasonable to provide just the initial conditions at , and the following complementary dynamics. In the interval : no loss of infectivity or immunity, just local contagion (only the first terms in (8)-(9)). And in : no loss of immunity (absence of the second term in (8)). Equations (6)-(7) can also be treated in this way.
These nonlocal SIRS epidemic models, at variance with the classic ones, display sustained oscillations. We showed in , by means of linear stability analysis and numerical computation of the solutions, that these oscillations appear as Hopf bifurcations in a variety of regions of parameter values. Observe, in Fig. 2, the line corresponding to , that represents the Hopf bifurcation in the fixed delays scenario (where and tend to Dirac deltas). The region of oscillations is above this curve. So, there is a transition in the dynamical behavior of the epidemic controlled by the basic reproduction rate of the infection, , as is usual in SIR and related models. And there is also a transition controlled by the ratio between the loss off immunity and infectious times, which is controlled by other factors of the disease and not by its reproduction rate (efficacy of vaccines, for example). When the kernels and are bell shaped instead of deltas, and there is some indeterminacy in the times of recovery and loss of immunity, the region of oscillations becomes enclosed by the curves shown in Fig. 2. There is a reentrance into the region of no oscillations in both directions of the control parameters. When the processes tend to constant rates (the kernels tend to exponentials, when ), the region of oscillations disappears towards the right in the figure. The kernels used in this example have the same shape, with , but the situation is the same with differently shaped and . Also, numerical solutions and stochastic simulations of the systems concur with these results, that correspond to a linear stability analysis.
The transition from no oscillations at constant rates, to oscillations at fixed delays, also happens as a Hopf bifurcation at a finite value of the shape parameters and . An example is shown in Fig. 3, where we plot the real part of the eigenvalue responsible for the destabilization of the fixed point as a function of the shape parameter . The inset shows three shapes of the kernel corresponding to increasing values of . It can be seen that the oscillations appear when this probability density is still very broad. this behavior represents a further transition in the dynamics of the epidemic system, this one controlled by the natural history of the disease.
3 The role of seasonality
The external forcing imposed by seasonality can be included in an epidemic model as a periodic oscillation of some of the parameters that control the dynamics. In our present model, defined by , and , the sensible choice is a varying contagion rate, . In the fixed delays equations this means:
with , oscillating around a mean value with an amplitude , and with a period that could be made equal to 1 year. In the region of oscillations it is expected an interplay of both effects: the oscillation inherent to the delayed dynamics, and that imposed externally through the varying contagion rate. This is, indeed, a situation that mimics real world epidemics in a much more realistic way than the simple model set by Eqs. (1)-(3).
Figure 4 shows the richness of behaviors that can be expected. The model is in the region of oscillations in the absence of external forcing, and defined by the parameters given in the caption. The plots show orbits (obtained by numerical computation) in the plane of susceptible-infected populations. Each panel corresponds to a different value of the period of oscillation of the contagion parameter . Some of the curves are closed and periodic, while others are quasiperiodic or seemingly chaotic. The orbits wind, naturally, due to the nonlocality in time, something that is impossible in the local models due to the uniqueness of solutions. As parameters change these windings change correspondingly, producing the observed range of solutions.
A compact representation of these transitions can be produced with a kind of Poincaré map, representing the time of return to each maximum of . In an epidemic system, this time is the recurrence of the outbreaks of the infection, and is generally well observed in the field. Figure 5 shows this for a model that represents a situation compatible with influenza (or ILI, influenza-like illness). The characteristic times are chosen as , , which, measured in days, correspond to such an epidemic scenario. is also compatible with influenza, while the amplitude of the oscillation has been left free and used as a control parameter. A simple recurrence of the epidemic is observed at low values of the amplitude, while a bifurcation cascade leads to more complex behaviors when . An extended region shows a period-2 behavior, which bears similarity to many cases observed in real systems. The time series of a period-2 solution is shown in Fig. 6. The successive outbreaks are of different magnitudes. The period is 2 years, as can be seen in the time between two large peaks, while the smaller peaks occur at a different moment of the year. The inset of Fig. 6 shows five consecutive years of cases of influenza-like illness in the South of Argentina , and it is easy to see that even years have their peaks around epidemic week 33 while odd years peak on week 25. The very large outbreak in 2009 corresponds to the epidemic of H1N1 virus (the “swine flu”).
We have studied generalized SIRS models with temporally spread infectious and immune processes. These models pose a better description of real epidemic systems than the usual mass-action equations at constant rates, providing a consistent treatment of the time delays inherent of infectious processes. Our model spans a whole range of scenarios in the delay kernels, going from constant rate equations to fixed delays with a single parameter. This allows to model infections that last a more or less fixed time, which is more realistic than constant rates. Even though the mathematical problem is more involved, a stability analysis of the stationary solutions can be made in many relevant cases. Besides this, numerical computation of the solutions, together with stochastic simulations of the processes, provide a reasonable set of tools for the analysis of these systems.
We have observed that the delayed models show oscillations in their solutions as a function of model parameters, such as the basic reproductive rate and characteristic times. Moreover, oscillations also appear as a function of the width of the characteristic times distributions, which represents the temporal indeterminacy of the corresponding processes. Finally, external forcing representing the seasonal variation of epidemiological parameters can be taken into account. In the case of an oscillating contagion rate, complex oscillations are observed that could be of relevance in the origin of some non-seasonal oscillations observed in real epidemic systems.
A cooperative agreement between the Brazilian Coordenação de Aperfeiçoamento de Pessoal de Nível Superior and the Argentinian Ministerio de Ciencia y Tecnología (grant CAPES-MINCyT 151/08-017/07) has supported our work. We also thank the Brazilian Conselho Nacional de Desenvolvimento Científico e Tecnológico (grant CNPq PROSUL-490440/2007), the Consejo Nacional de Investigaciones Científicas y Técnicas (PIP 112-200801-00076), and Universidad Nacional de Cuyo (06/C304).
-  R. M. Anderson and R. M. May, “Infectious Diseases of Humans: Dynamics and Control,” Oxford University Press, Oxford, 1991.
-  J. D. Murray, “Mathematical Biology”, Springer-Verlag, Berlin, 1993.
-  N. C. Grassly, C. Fraser and G. P. Garnett, Host immunity and synchronized epidemics of syphilis across the US, Nature, 433 (2005), 417–421.
-  B. Grenfell and O. Bjørnstad, Epidemic cycling and immunity, Nature, 433 (2005), 366–367.
-  A. M. Fry, A. T. Curns, K. Harbour, L. Hutwagner, R. C. Holman and L. J. Anderson, Seasonal trends of human parainfluenza viral infections: United States, 1990 2004, Clin. Inf. Dis., 43 (2006), 1016–1022.
-  M. N. Kuperman, personal communication, based on data from the Argentinian National Ministry of Health.
-  G. Abramson and V. M. Kenkre, Spatio-temporal patterns in the Hantavirus infection, Phys. Rev. E, 66 (2002), 011912.
-  S. Risau-Gusmán and G. Abramson, Bounding the quality of stochastic oscillations in population models, Eur. Phys. J. B, 60 (2007), 515–520.
-  J. P. Aparicio and H. G. Solari, Sustained oscillations in stochastic systems, Math. Biosciences, 169 (2001), 15–25.
-  M. N. Kuperman and G. Abramson, Small world effect in an epidemiological model, Phys. Rev. Lett., 86 (2001), 2909–2912.
-  H. W. Hethcote, H. W. Stech and P. Van Den Driessche, Nonlinear oscillations in epidemic models, SIAM J. Appl. Math., 40 (1981), 1–9.
-  S. Gonçalves, G. Abramson and M. F. C. Gomes, Oscillations in SIRS model with distributed delays, Eur. Phys. J. B, 81 (2011), 363- 371.