Seasonal synchronization and unpredictability in epidemic models with waning immunity and healthcare thresholds

Categories: Disease & Virus

May 18, 2025

Limited treatment model

We model the dynamics of recurrent epidemics due to limited treatment by a system of differential equations, which represents a compartmental SIRS (Susceptible-Infectious-Removed-Susceptible) model with a nonlinear recovery rate. We omit the exposed compartment of infected and non-infectious individuals due to the simplicity of the model, since adding this compartment has no effect on the non-linear mechanisms we study. The following parameters can be defined:

Parameters:

$\beta$: The transmission rate, specifying the probability of disease transmission per time from infectious to susceptible individuals which includes the average number of contacts per person per time and the probability of transmission.
$\gamma$: The recovery rate, indicating the fraction of infected individuals recovering without treatment per unit of time.
$f(K, \eta , r, I)$: The function describing acceleration of removed rate. The function is defined by following variables and parameters: K – a population threshold (healthcare availability), $\eta$ – a maximal acceleration rate, r – a steepness by which the threshold K is reached, and I – the current infectious proportion. For more details see “Acceleration of removed rate” section.
$\nu$: The waning immunity rate, representing the rate at which recovered individuals lose immunity or exit the removed compartment.

Equations: The system of differential equations is as follows:

$$\begin{aligned} \begin{aligned} \frac{dI}{dt}&= \beta S I – \gamma f(K,\eta ,r,I) I, \\ \frac{dR}{dt}&= \gamma f(K,\eta ,r,I) I – \nu R, \\ S&= 1 – I – R. \end{aligned} \end{aligned}$$

(1)

The selection of system parameters in our study is based on a conceptual approach rather than direct experimental data. However, we made an effort to align this conceptual principle with realistic parameter values to enhance the understanding for researchers who work with epidemiological data. Our predictive model²³, which was also part of the European Centre for Disease Prevention and Control (ECDC) ensemble²⁴ with high predictive accuracy, demonstrated our experience with real data, where we thoroughly discussed the choice of the SIR(S) model and the methodology for parameter estimation. Consequently, the analysis presented in this paper is conducted within a parameter range that corresponds to a realistic epidemiological scenario.

Parameters are chosen to evoke the COVID-19 disease²³. The transmission rate, $\beta$, is set to small hundreds, which corresponds to the serial interval in days²⁵. It is heavily influenced by the contact rate; however, we assume an average long-term contact rate within the community, excluding short-term periods of non-pharmaceutical interventions, such as lockdowns and mandatory masking. Since both the contact rate and the serial intervals are stochastic, we will also demonstrate how this stochasticity interacts with seasonality. Specifically, we will explain four different mechanisms of unpredictability in “Unpredictability and stochastic differential equations” section. Health-seeking behavior is strongly associated with the severity of illness²⁶, while the severity of COVID-19 significantly depends on age and comorbidities. Therefore, we assume that the percentage of people seeking healthcare or isolating themselves is small, specifically, K is in small units of percentages. The recovery rate, $\gamma = 36$, indicates an approximately 10-day infectious period^27,28. The rate of waning immunity $\nu$ can be considered in small units or could be even higher for post-delta SARS-CoV-2 variants, which have a high percentage of reinfections and evidence of post-vaccination and post-infection protection decline²⁹. For simplicity, we omit mortality by the disease assuming this compartment would be very small with respect to the susceptible compartment.

Acceleration of removed rate

After infection with an infectious disease, individuals may experience a spectrum of symptoms, ranging from mild to severe (for COVID-19, see meta-analyses^30,31, clinical studies^32,33 or nation-wide data analyses^29,34). This variation is evident in the choices individuals make regarding seeking healthcare or entering isolation, with those experiencing severe symptoms more frequently opting to seek medical attention²⁶. If they take this course of action, they can either recover or, at the very least, be isolated, preventing further spread of the disease. It is worth noting that the proportion of the population at risk of developing severe illness is probably relatively small²⁶.

From a mathematical perspective, a higher removal rate $\gamma$ is associated with the proportion of people seeking healthcare or isolating. We follow a model², where acceleration was defined using a piece-wise constant function. We generalize this approach by using its standard continuous approximation. The function

$$\begin{aligned} f(K, \eta , r, I) = \frac{\eta + \eta e^{-r (I-K)}}{\eta + e^{-r (I-K)}} \end{aligned}$$

(2)

governing the acceleration of the removal rate $\gamma$ is determined by the following variables and parameters: K, signifying the proportion of individuals seeking healthcare or isolation; $\eta$, representing the maximum acceleration rate; r, indicating the steepness with which the threshold K is attained; and I, denoting the current infectious proportion.

Fig. 1 illustrates the function f across various values of the parameter r. It is evident that, for higher values of r, the function tends towards a step function. This implies the existence of two potential removal rates $\gamma$. From that perspective, the function f is considered as an approximation of the step function. When analyzing the model using continuation software like MatCont^21,22, it is crucial for the function to be continuous. Also in reality, the function is likely not a step function at some given threshold K; rather, the threshold K may exhibit some variation, and a continuous function is a better fit for such variations. For the purpose of this paper we use $r=100$.

Limited treatment and seasonal transmission rate model

Following a previous research⁷, the transmission rate $\beta$ may exhibit variation over the course of a year, influenced by seasonal changes. This variation can be attributed to the seasons, where during the summer, the transmission rate tends to be lower due to increased room ventilation, reduced indoor occupancy, enhanced outdoor activities, and decreased contact rate resulting from summer holidays, etc. In contrast, during winter, individuals tend to stay indoors, and the colder weather makes them more susceptible to the disease, leading to an increase in the transmission rate.

Therefore $\beta (t)$ is time-dependent periodic function:

$$\begin{aligned} \beta (t)=\beta _0(1+\alpha \cos (2\pi T t)), \end{aligned}$$

(3)

where $\alpha \ge 0$ is the amplitude of the seasonal oscillation as percentage of $\beta _0$, $T=1$ is the period of function $\beta (t)$, and t is time with year-long scaling. Using time-dependent transmission rate, described by Eq. (3), we obtain system of following differential equations

$$\begin{aligned} \begin{aligned} \frac{dI}{dt}&= \beta _0 \left( 1 + \alpha \cos (2\pi T t)\right) S I – \gamma f(K,\eta ,r,I) I, \\ \frac{dR}{dt}&= \gamma f(K,\eta ,r,I) I – \nu R, \\ S&= 1 – I – R. \end{aligned} \end{aligned}$$

(4)

Embedding

We model the effect of seasonality with a harmonic oscillation driving the epidemic sub-system, which naturally displays a limit cycle (as seen in Eq. (4) and reference²), resulting in a non-autonomous system. We can utilize a previously described embedding method^35,36,37 for harmonically driven oscillators with ease to a generalized higher-dimensional autonomous system. Suppose the system can be expressed in the following form:

$$\begin{aligned} \dot{\textbf{x}} = \textbf{f}(\textbf{x}, \alpha \cos (2\pi T t), \pmb {\beta }), \end{aligned}$$

(5)

where $\textbf{f}:\mathbb {R}^{m+p+1} \rightarrow \mathbb {R}^m$ is a function smooth enough, and $\pmb {\beta } \in \mathbb {R}^p$ are given parameters. $\alpha \cos (2\pi T t)$ denotes a harmonic external force with amplitude $\alpha$ and frequency T. The system (5) can be embedded into an autonomous system:

$$\begin{aligned} \begin{aligned} \dot{\textbf{x}}&= \textbf{f}(\textbf{x}, u, \pmb {\beta }),\\ \dot{u}&= 2\pi T(A u – v – u(u^2 + v^2)),\\ \dot{v}&= 2\pi T(u + A v – v(u^2 + v^2)), \end{aligned} \end{aligned}$$

(6)

The two-dimensional driving sub-system incorporated for $(u,v)^T$ is in the normal form of a super-critical Hopf bifurcation. When A is negative, it possesses a stable equilibrium at the origin, while for positive A, a stable limit cycle of the form $S = \{(u,v):u^2+v^2=A\}$ emerges. Consequently, when $A = \alpha ^2$, the dynamics of the system (6) on the invariant manifold $\mathbb {R}^m\times S$ become topologically equivalent to the system (5) dynamics. The synchronization of a driven oscillator (the epidemic endemic cycle) with driving oscillator (the seasonal cycle) in the multidimensional state space corresponds to the emergence of a unique limit cycle, while the non-synchronous oscillations of the two oscillators create a quasi-periodic trajectory on the torus or more complex dynamics. The asymptotic stability of cycle S (driving oscillator) ensures favorable numerical properties for the continuation of bifurcation manifolds^35,36,37, particularly in software packages like MatCont^21,22. Additionally, the autonomous nature of the system enables the computation of the Lyapunov spectrum, which will be discussed in the following “Analysis techniques” section.

In this case, we can embed model (4) into a new generalized system

$$\begin{aligned} \begin{aligned} \frac{dI}{dt}&= \beta _0 (1 + u) S I – \gamma f(K,\eta ,r,I) I, \\ \frac{dR}{dt}&= \gamma f(K,\eta ,r,I) I – \nu R, \\ \frac{du}{dt}&= 2\pi T(A u – v – u(u^2 + v^2)), \\ \frac{dv}{dt}&= 2\pi T(u + A v – v(u^2 + v^2)), \\ S&= 1 – I – R, \\ A&= \alpha ^2 \text { for } A \ge 0. \end{aligned} \end{aligned}$$

(7)

Analysis techniques

We provide bifurcation analysis of system (1) in MatCont^21,22, version 7p4. Bifurcation refers to a qualitative change in the dynamics of the studied system as a parameter or set of parameters varies. These changes often involve, for example, the emergence, destruction, or alteration of equilibrium points, periodic orbits, or other invariant sets^18,19. In this paper, we primarily focus on the emergence of cycles. We interpret the limit cycle dynamics as a periodic epidemic outbreaks¹.

It is crucial to note that bifurcation analysis, although computed numerically, fundamentally differs from simulation. Instead of observing the system’s behavior for specific parameter values, we use a continuation approach to systematically trace the boundaries where qualitative changes occur. This continuation process involves formulating the detection of bifurcation points as a system of algebraic equations, which are then solved using a Newton-like predictor-corrector method. This methodology allows us to accurately and efficiently map out bifurcation boundaries within the parameter space.

In epidemiology, researchers often analyze the basic reproduction number, which indicates whether a disease has the potential to spread. If this number exceeds one, the disease can spread. From the perspective of bifurcation theory, this phenomenon can be elucidated through the transcritical bifurcation. Transcritical bifurcation involves a switch in stability between two equilibrium states. In this scenario, it occurs between a trivial equilibrium (representing a disease-free state) and an endemic equilibrium. Transcritical bifurcation is not the only mechanism how the basic reproduction number exceeds one (see¹), but it is the only relevant one for our model (1).

As our primary aim is to study the emergence of cycles, we focus on one-parameter bifurcations starting with the Hopf bifurcation. In its super-critical version, it allows a stable limit cycle to arise from a stable equilibrium, that changes stability on the Hopf manifold. The sub-critical version has exactly opposite stability, it gives birth of an unstable limit cycle from unstable equilibrium. Along the Hopf manifold, we identified generalized Hopf points (GH). At these generalized points, the stability of the studied limit cycle changes, and a region of bistability between the stable limit cycle and the endemic equilibrium arises ¹⁹.

A second mechanism of a limit cycle emergence is the limit point of cycles bifurcation when a pair of limit cycles, one stable and one unstable, collide and disappear. The existence of this bifurcation is guaranteed at least in a neighborhood of the two parameter GH bifurcation¹⁹.

A third mechanism is the saddle-node infinite cycle (SNIC) bifurcation, which is characterized by limit cycle, which period becomes longer and longer, finally it splits on a homoclinic orbit. This bifurcation is guaranteed at least in a neighborhood of the two-parametric Bogdanov-Takens bifurcation¹⁹.

By introducing seasonal transmission rate model (4) the dynamics become more complex. We study an interaction of two periodic forces in the model (4), one given by seasonal transmission rate and the second by the existence of limit cycle dynamics in model (1). We identified there important dynamical scenarios:

1.

The model may undergo a cascade of period-doubling and end in chaotic regime.
2.

The model may exhibit dynamics on an invariant torus. Seasonally synchronized cycles in case of strong enough seasonality in transmission rate/quasiperiodicity if not.
3.

Due to bistability in model (1), we may encounter bistability also in model (4).

To investigate the phenomena described above, we employ the Poincaré sections³⁸ and the Lyapunov exponent spectrum^39,40, in addition to bifurcation analysis of system (7) in MatCont^21,22, version 7p4.

A Poincaré section provides a method to convert a continuous dynamical system into a discrete one. This process involves selecting a lower-dimensional surface, known as the Poincaré section, which intersects the trajectory of the continuous system. Equivalently, the selection of this section can be guided by criteria such as maximum and minimum values of certain variables along the trajectory, or a section at given times of the period. The points where the trajectory intersects this section are used to define a discrete sequence of states for the system.

Lyapunov exponents are essential mathematical measures that quantify the divergence or convergence of trajectories in dynamical systems, reflecting the system’s sensitivity to initial conditions. When two points are initially proximate in phase space, a positive Lyapunov exponent denotes exponential divergence in their separation, indicating chaotic dynamics, whereas a negative exponent implies convergence and stability. Understanding the mechanisms underlying the emergence and disappearance of chaotic dynamics in infectious disease outbreaks is crucial. Early investigations from the initial stages of the COVID-19 pandemic have pointed to the potential presence of such behavior⁴¹. In one-dimensional systems; there is generally a single Lyapunov exponent; however, higher-dimensional systems possess a Lyapunov spectrum, which comprises a collection of exponents indicative of stability over all directions. Lyapunov exponents offer valuable insights into a system’s stability and qualitative behavior, particularly near bifurcations that may lead to complex dynamics such as chaos; shifts in their sign or magnitude can signal critical transitions even when standard continuation methods are not feasible. To illustrate this, the following table 1 summarizes typical dynamical attractors along with their characteristic Lyapunov exponents and associated minimal embedding dimensions.

Table 1 Common types of dynamical attractors with their characteristic Lyapunov exponents and the minimum phase space dimension required for their occurrence.

To quantify the Lyapunov spectrum³⁹ for a nonlinear system, one generally starts by numerically simulating the system’s dynamics while concurrently monitoring the evolution of minor perturbations to the system’s state across time. This entails the integration of the original nonlinear system alongside its corresponding variational equations, which describe the linearised dynamics of perturbations around a trajectory. The rate of progression of these perturbations is determined by the system’s Jacobian matrix, computed along the trajectory at each time interval. To ensure numerical stability and avoid vector collapse, the Gram-Schmidt orthonormalization method (or QR decomposition) is often applied to the evolving perturbation vectors. With time, the average exponential growth or decay rates of these orthonormalised vectors converge to the Lyapunov exponents, resulting in the whole spectrum. We employ this technique to compute the spectrum by Govorukhin^40,42.

In the last section, we introduce stochasticity into the transmission rate in the model (7). We study the following set of stochastic differential equations:

$$\begin{aligned} \begin{aligned} \text {d}I&= \beta _0 (1 + u) S I \text {d}t – \gamma f(K,\eta ,r,I) I \text {d}t + \sigma \beta _0 S I \text {d}W\\ \text {d}R&= \gamma f(K,\eta ,r,I) I \text {d}t – \nu R \text {d}t, \\ \text {d}u&= 2\pi T(A u – v – u(u^2 + v^2)) \text {d}t, \\ \text {d}v&= 2\pi T(u + A v – v(u^2 + v^2)) \text {d}t, \\ S&= 1 – I – R, \\ A&= \alpha ^2 \text { for } A \ge 0, \end{aligned} \end{aligned}$$

(8)

where W denotes Wiener process and $\sigma$ measures transmission rate non-seasonal stochastic variability, and simulate them using Milstein method⁴³.

Source link