The effects of heterogeneity on stochastic cycles in epidemics

TL;DR

This paper analyzes how heterogeneity in susceptibility and infectivity influences stochastic epidemic cycles, deriving analytical spectra and examining synchronization effects in a modified SIR model.

Contribution

It introduces a reduced analytical framework for understanding the impact of heterogeneity on epidemic quasi-cycles and their spectral properties.

Findings

Cycle frequency depends on population averages of heterogeneity.

Amplitude of cycles is influenced by higher-order moments of heterogeneity.

Heterogeneity affects synchronization and phase lag between groups.

Abstract

Models of biological processes are often subject to different sources of noise. Developing an understanding of the combined effects of different types of uncertainty is an open challenge. In this paper, we study a variant of the susceptible-infective-recovered model of epidemic spread, which combines both agent-to-agent heterogeneity and intrinsic noise. We focus on epidemic cycles, driven by the stochasticity of infection and recovery events, and study in detail how heterogeneity in susceptibilities and propensities to pass on the disease affects these quasi-cycles. While the system can only be described by a large hierarchical set of equations in the transient regime, we derive a reduced closed set of equations for population-level quantities in the stationary regime. We analytically obtain the spectra of quasi-cycles in the linear-noise approximation. We find that the characteristic frequency of these cycles is typically determined by population averages of susceptibilities and infectivities, but that their amplitude depends on higher-order moments of the heterogeneity. We also investigate the synchronisation properties and phase lag between different groups of susceptible and infected individuals.

Figures

Equations

$$S_i\stackrel{\xi {\chi }_i{q}_a}{⟶}{I}_a$$

$$S_i+I_a\stackrel{{\beta }_a{\chi }_i{q}_b}{⟶}{I}_a+I_b$$

$$I_a\stackrel{\rho }{⟶}R$$

$$S_j\stackrel{p_i\kappa }{\to }{S}_i$$

$$\overline{\chi }=\sum _i{p}_i{\chi }_i,\text{and }\mathfrak{X}=\frac{1}{N}\sum _i{\chi }_i{n}_i.$$

$$\overline{\beta }=\sum _a{q}_a{\beta }_a\text{and }\mathcal{B}=\frac{1}{N}\sum _a{\beta }_a{m}_a,$$

$$\dot{x_i}$$

$${\dot{y}}_a$$

$$S=\sum _i{x}_i\text{and }I=\sum _a{y}_a.$$

$$\dot{S}$$

$$\dot{I}$$

$$\dot{\mathfrak{X}}$$

$$\dot{\mathcal{B}}$$

$${\dot{\mathfrak{X}}}_n=\kappa \overline{{\chi }^n}-\kappa {\mathfrak{X}}_n-(\xi +\mathcal{B}){\mathfrak{X}}_{n+1},$$

$$x_i^{\star }$$

$$y_a^{\star }$$

$$S^{\star }$$

$$I^{\star }$$

$${\mathfrak{X}}^{\star }$$

$${\mathcal{B}}^{\star }$$

$${\dot{\tilde{x}}}_i$$

$${\dot{\tilde{y}}}_a$$

$$\dot{\tilde{S}}$$

$$\dot{\tilde{I}}$$

$$\dot{\tilde{\mathfrak{X}}}$$

$$\dot{\tilde{\mathcal{B}}}$$

$${\hat{x}}_i$$

$${\hat{y}}_a$$

$$\hat{S}$$

$$\hat{I}$$

$$\hat{\mathfrak{X}}$$

$$\hat{\mathcal{B}}$$

$$A_i$$

$$C$$

$$D$$

$$E$$

$${\mathcal{P}}_{\mathcal{B}}(\omega )$$

$$|g{|}^2={\left[E+\overline{\beta }C\kappa \sum _i\frac{{\chi }_i^2{p}_i}{{\omega }^2+A_i^2}\right]}^2+{\omega }^2{\left[1-\overline{\beta }C\kappa \sum _i\frac{{\chi }_i^2{p}_i}{A_i\left({\omega }^2+A_i^2\right)}\right]}^2.$$

$${\omega }_d\approx \sqrt{\kappa \overline{\chi }\overline{\beta }}.$$

$${\mathcal{P}}_I\left({\omega }_d\right)$$

$${\mathcal{P}}_S\left({\omega }_d\right)$$

Full text

Introduction

It is now widely recognised that noise and uncertainty play an important role in modelling biological systems. Traditional approaches to modelling phenomena in biology[1] are often based on deterministic ordinary or partial differential equations, and do not aim to describe stochasticity. In order to capture epistemic uncertainty, static or dynamic noise variables are introduced in more modern mathematical biology. This randomness reflects the lack of detailed knowledge about phenomena at finer scales than described by the model at hand; any modelling approach necessarily operates at a set scale (e.g. cell, individual, or population), and does not capture in detail the processes at smaller scales. These are ‘emulated’ through effective randomness. Different types of such noise are frequently found in models of biological phenomena, including intrinsic demographic noise, extrinsic stochasticity, parameter uncertainty or heterogeneity between different types of interacting entities [2, 3]. Some of these random variables are static and do not evolve in time, others are described by dynamic time-dependent noise. Intrinsic noise, due to the stochastic dynamics of a system has lately been the focus of many studies (see for example [4, 5, 6]). Extrinsic or parametric noise, due to variations, heterogeneity or uncertainties in the parameters or the environment surrounding the process, has received similar attention (e.g. [7, 8]). To be able to adequately describe biological systems, however, it may be necessary to account for both these uncertainties which contribute to the noisy dynamics.

In the modelling of epidemics this is of particular importance. The infection process, driven by serendipitous contacts, is inherently stochastic, and heterogeneity in susceptibility to a disease or infectiousness of different individuals are known to exist and play a role in viral spread. For example, variation in host susceptibility and viral reproduction have been observed in [9], and behavioural, structural or contact differences between individuals are inevitable. However, the better part of the existing work focusing on heterogeneity of this type, does not explicitly seek to capture demographic noise. Instead one often assumes infinite populations and deterministic dynamics. This approach is often taken outside epidemics as well. Much existing work studies individual sources of uncertainty, heterogeneity and noise in isolation, but not their interacting together. A notable exception is the modelling of gene regulatory networks, in which the interaction of intrinsic and extrinsic noise is actively studied, see e.g. [10, 11, 12].

The effects of intrinsic noise have been recognised in recent years. In models with demographic processes, for example, intrinsic stochasticity has been seen to lead to sustained quasi-cycles [13, 14, 15, 16] in parameter regimes in which a deterministic model would converge to a stable fixed point. These quasi-cycles have been identified not only in models of epidemic spread, but also in other instances of population dynamics, including in genetic circuits, evolutionary systems and in game theory [17, 18, 19, 20]. Heterogeneity has been and is being considered in epidemics as well. Age structure is studied for example in [21, 22], seasonally changing infection rates in [23, 24], variation in infectivity and/or susceptibility are addressed in [25, 26, 27, 28, 29], spatial structure has been approached in [30, 31, 32, 33], and epidemics on static and dynamic networks are studied in [34, 35, 36, 7, 37, 38]. Heterogeneity has been found to generate outbreaks that propagate hierarchically [35, 39], grow faster than in homogeneous populations [36], and have a lower total number of infected individuals [40, 41].

Much of this work, whether describing a well-mixed population, a compartmented or structured one, is based on variants of the celebrated susceptible-infective-recovered (SIR) model. They can be described either by deterministic differential equations, or as a stochastic process involving a population of discrete individuals. In the former approach the population is effectively assumed to be infinite, so that the timing of stochastic infection, recovery or birth-death events ‘averages’ out, and smooth laws for the time evolution of the population are obtained. The latter approach explicitly captures the intrinsic randomness of infection, recovery and demographics. The population is taken to be finite, and its state discrete. The model evolves through discrete events (e.g. infections). In the simplest case this defines a Markovian random process, which often can be analysed further mathematically, at least to a good approximation. Starting from the master equation in a well-mixed population a set of stochastic differential equations can be derived in the limit of large, but finite populations [42]. These can then be studied further within the ‘linear-noise approximation’ (LNA) [43]. The mathematics are tractable and the corresponding theory is now well established. While remarkably powerful, this approach so far has mostly been used for well-mixed populations. The linear-noise approximation has also been applied to networked systems with contact heterogeneity (see e.g. [16, 44]), but progress is then much harder and often relies on further moment-closure approximations.

The aim of our work is to introduce agent-to-agent heterogeneity into the SIR dynamics in a finite well-mixed population. This provides a middle ground between homogeneous well-mixed models and an explicitly networked population. At the same time, we maintain tractability and are able to characterise stochastic effects in finite populations via the linear-noise approximation. This allows us to systematically investigate the combination of parameter heterogeneity and demographic noise. We divide the population of agents into $K$ different groups of susceptible individuals, where members of different groups have different susceptibilities. Similarly, in our model there are $M$ classes of infective individuals, with each class representing a different propensity to pass on the disease. This follows the lines of [29], but we explicitly focus on the combination of heterogeneity and intrinsic noise. Intrinsic stochasticity had not been included in [29].

Our paper is organised as follows: In Sec. 1 we describe our model in detail. As a baseline we then construct the deterministic rate equations in Sec. 2. They describe the deterministic dynamics in the limit of infinite populations, and are required to carry out the LNA. The most natural deterministic description will generally involve $K+M$ coupled non-linear equations (one for each subclass in the population). We discuss when and how these can be reduced to a smaller set of equations for aggregate quantities. In Sec. 3 we perform then the linear-noise approximation and use this approximation to characterise the fluctuations about deterministic fixed points. In particular we set up the theory to obtain the spectra of noise-driven quasi-cycles. Using this theory we then present our main results in Sec. 4, where we investigate in detail how the heterogeneity in the population affects the properties of stochastic outbreaks of the disease. Finally, in Sec. 5 we summarize our findings.

1 Model

We use an extension of the standard SIR model [45], in a population of fixed size $N$. Broadly, each individual can be of one of three types, susceptible (S), infective (I) or recovered (R). The spreading of the disease is described by infection events. These occur either through contact of a susceptible with an infective individual, as described below, or through spontaneous infection. Individuals recover at rate $\rho$, and they die at rate $\kappa$. The death rate is assumed to be independent of the disease status of an individual. To keep the number of individuals in the population constant, any death event is immediately followed by a birth of a new susceptible individual. This modelling assumption is made for simplicity and is commonly made (see e.g. [46, 13, 47]).

We introduce heterogeneity by dividing the groups of susceptibles and infectives into subclasses. We will write $S_{i}$ and $I_{a}$ for these, with $i=1,\dots,K$ and $a=1,\dots,M$. Individuals in subgroup $S_{i}$ have susceptibility $\chi_{i}$ to the disease, and infectives in class $I_{a}$ have infectiousness $\beta_{a}$, which describes the propensity of the infective to pass on the disease to susceptible individuals. We write $n_{i}$ for the number of individuals of type $S_{i}$, and $m_{a}$ for the number of individuals in class $I_{a}$.

The dynamics are illustrated in Fig. 1, and can be summarised in the following reaction scheme:

[TABLE]

where $\{p_{i}\}$ and $\{q_{a}\}$ represent the probabilities of being assigned a susceptibility $\chi_{i}$ or infectiousness $\beta_{a}$ at birth or upon infection, respectively. The first of these reactions describes spontaneous infection, converting an individual in class $S_{i}$ into an individual of type $I_{a}$. The per-capita rate of events of this type is $\xi\chi_{i}q_{a}$, where $\xi$ is an overall inverse time scale for spontaneous infection, $\chi_{i}$ is the susceptibility of $S_{i}$ to the disease, and $q_{a}$ is the probability that the newly infected individual is in class $I_{a}$. Similarly, the second reaction describes infection of an individual of type $S_{i}$ upon contact with an individual of type $I_{a}$. The newly infected individual is in class $I_{b}$. Events of this particular type occur with a rate proportional to $\beta_{a}$ (the propensity of $I_{a}$ to spread the disease), to $\chi_{i}$ (the susceptibility of $S_{i}$) and to $q_{b}$. The third reaction describes recovery, and the final three reactions are birth/death events. The newly born individual is assumed to be randomly placed into one of the classes $S_{i}$ ($i=1,\dots,K$), occurring with respective probability $p_{i}$. We note that our model does not describe potential correlations between the susceptibility of an individual and its infectivity after they become infected; our focus is on heterogeneity of susceptibility due to physiological factors, and not primarily due to contact patterns. Extensions to include correlations can however be constructed among similar lines.

The model defines a continuous-time Markov process, and can be simulated straightforwardly using for example the celebrated Gillespie algorithm [48]. The starting point for the analytical study of the model is the master equation. Our analysis below will be based on approximating the solution to this master equation by performing a system-size expansion [43] and linear-noise approximation, leading to a stochastic differential equation describing the dynamics in the limit of large, but finite population size.

In order to do this it is useful to first introduce

[TABLE]

The quantity $\overline{\chi}$ is the mean susceptibility of a newly born individual, whereas $N\mathfrak{X}$ describes the aggregate susceptibility of the population. Similarly we define

[TABLE]

where $\overline{\beta}$ represents the mean infectivity of a newly infected individual, and $N\mathcal{B}$ the total ‘infective power’ in the population. We note that $\overline{\chi}$ and $\overline{\beta}$ are fixed in time, and are properties of the distributions $\{p_{i},\chi_{i}\}$ and $\{q_{a},\beta_{i}\}$. The quantities $\mathfrak{X}$ and $\mathcal{B}$, on the other hand, are time-dependent and evolve as the composition of the population changes.

2 Deterministic analysis

2.1 Dynamics

In the limit of an infinite population the dynamics can be described by deterministic equations for the quantities $x_{i}=\lim_{N\to\infty}n_{i}/N$, $y_{a}=\lim_{N\to\infty}m_{a}/N$. They are given by

[TABLE]

These ordinary differential equations can be derived either by using direct mass-action kinetics, or from the lowest-order expressions in an expansion of the master equation in the inverse system size [43].

Ultimately we will mostly be interested in aggregate quantities, i.e. the total density of susceptibles or infectives in the population, irrespective of what subclass they belong to. We therefore introduce

[TABLE]

From Eqs. (4) we find

[TABLE]

This system is not closed due to the presence of $\mathfrak{X}$ and $\mathcal{B}$ on the right-hand side. These quantities in turn evolve in time according to

[TABLE]

which again does not close the set of equations, due to the presence of the term $\mathfrak{X}_{2}(t)\equiv\sum_{i}\chi_{i}^{2}x_{i}(t)$. Modulo normalisation and recalling that the $\{x_{i}\}$ are time-dependent, this object is recognised as the second moment of the distribution of susceptibilities among the group of susceptibles at time $t$. It cannot be determined from Eqs. (6) and (7) alone. Instead we find

[TABLE]

where we have introduced $\overline{\chi^{n}}=\sum_{i}p_{i}\chi_{i}^{n}$ and $\mathfrak{X}_{n}=\sum_{i}x_{i}\chi_{i}^{n}$. This indicates that the deterministic dynamics at the aggregate level is described by an infinite hierarchy of equations. This set of equations does not close in the transient regime. However, as we will see next, closure can be achieved assuming the system settles down to a fixed point in the long run.

2.2 Fixed point

We proceed by a brief analysis of the fixed points of the deterministic dynamics. We will label these by a star. They can be obtained by setting $\dot{x}_{i}=0$ and $\dot{y}_{a}=0$ in Eqs. (4), leading to

[TABLE]

Similarly, we find the fixed points of the aggregate quantities $S$, $I$, $\mathfrak{X}$ and $\mathcal{B}$ from Eqs. (6,7). After re-arranging and using Eqs. (9) we arrive at

[TABLE]

which is a closed set of equations, for a given set of model parameters $\{p_{i},\chi_{i},q_{a},\beta_{a}\}$.

We highlight that while the transient dynamics of the system described in terms of the four macroscopic variables $S$, $I$, $\mathfrak{X}$ and $\mathcal{B}$ generates an infinite hierarchy of equations, potential fixed points can be uniquely described by a closed set of equations, assuming that the distribution of susceptibilities at birth and of the propensity of newly infected individuals to pass on the disease are known. In other words, the fixed point can be obtained in terms of the model parameters $\{q_{a},\beta_{a}\}$ and $\{p_{i},\chi_{i}\}$. While we cannot provide an analytical proof that the deterministic system will always converge to a fixed point, we note that, for the range of parameter used, we have not detected a single case in which numerically integrating Eqs. (4) did not lead to a fixed point. In this context it is useful to point out that, in a homogeneous model, any combination of susceptibility and infectivity within the range of parameters used here would lead to a basic reproductive number above unity. For such models it is known that stable fixed points are eventually reached [49].

3 Linear-noise approximation

We now proceed to analyse the effects of stochasticity in the model, with a particular focus on the interaction between heterogeneity of individuals in the population and the noise induced by the demographics of the finite system.

We illustrate these effects in Fig. 2, and show an example of both the deterministic time-evolution of the system (thick continuous lines) and a realization of an individual-based simulation (thin dashed lines); the latter illustrates the intrinsic stochasticity of the process. Even after the deterministic model has reached a fixed point, the individual-based model shows sustained oscillations around it. We will focus our attention on these stochasticity-driven periodic outbreaks in the remainder of this article. In particular we will study how the heterogeneity in the population affects their properties.

3.1 Stochastic Dynamics

In order to carry out an analysis of the stochastic dynamics, we write $n_{i}/N=x_{i}+\tilde{x}_{i}/\sqrt{N}$, and $m_{a}/N=y_{a}+\tilde{y}_{a}/\sqrt{N}$, where $x_{i}(t)$ and $y_{a}(t)$ are the solutions of the deterministic equations (4) and the quantities with a tilde describe the stochastic fluctuations about the deterministic trajectory. The above ansatz reflects the anticipation that these fluctuatons will have a relative magnitude of order $N^{-1/2}$. We then carry out an expansion in the inverse system size up to and including sub-leading order [43], and arrive at

[TABLE]

The $\{\eta_{i}\}$ and $\{\nu_{a}\}$ are Gaussian white noise variables, with variance and co-variance (across components) as described in more detail in the Supplement (see S1). Writing $\tilde{S}=\sum_{i}\tilde{x}_{i}$ and $\tilde{I}=\sum_{a}\tilde{y}_{a}$ we find the following dynamics of fluctuations at the aggregate level,

[TABLE]

As in the deterministic analysis, this set of equations for the transient dynamics is not closed. This is due to the term $\sum_{i}\chi_{i}^{2}\tilde{x_{i}}$ in the equation for $\dot{\tilde{\mathfrak{X}}}$. However, as in Section 2.2, we will show below that a closed set of equations for fluctuations in the stationary state can be derived.

3.2 Fluctuation around the deterministic fixed point

We here show that although Eqs. (3.1) are not closed, we can explore noise-induced oscillations around the deterministic fixed point. To this end we introduce the Fourier transforms (with respect to time) of the variables $\tilde{x}_{i}$ and $\tilde{y}_{a}$. We will denote these by $\widehat{x}_{i}$ and $\widehat{y}_{a}$. From the Langevin equations (3.1) we find, after re-arranging,

[TABLE]

The noise variables $\{\eta_{i}\}$ and $\{\nu_{a}\}$ are uncorrelated in time, and their variance and correlation across components can be expressed in terms of known quantities (see Eqs. (S1) in the Supplement). The variable $\omega$ is the conjugate of time under Fourier transform. Similarly, we find the following for the relevant aggregate quantities,

[TABLE]

where, for simplicity, we have introduced the notation

[TABLE]

Eqs. (3.2) constitute a closed set of equations for the Fourier transforms of the aggregate fluctuations $\tilde{S},\tilde{I},\tilde{\mathfrak{X}}$ and $\tilde{\mathcal{B}}$ in the stationary state. We thus make an observation similar to that in Section 2: although we cannot describe the evolution of fluctuations in the transient regime, we can derive a closed description of the statistics of fluctuations about deterministic fixed points within the linear-noise approximation.

3.3 Power Spectral Density

Eqs. (3.2) can be used describe the periodic cycles shown in Fig. 2; we will now proceed to analyse these in more detail. Specifically we will use the above results to compute the power spectral density (PSD) of fluctuations. This allows us to identify the characteristic frequency of noise-driven epidemic cycles, and to infer information about their amplitude.

The (average) power spectral density of a time series, $z(t)$, generated from the stochastic individual-based model is given by $\mathcal{P}_{z}(\omega)=\langle|\widehat{z}(\omega)|^{2}\rangle$, where $\left\langle{\cdots}\right\rangle$ stands for an average over realizations of the stochastic dynamics. The PSD can be computed analytically for all individual signals $x_{i}$, $y_{a}$, and for the aggregate variables $S$, $I$, $\mathfrak{X}$ and $\mathcal{B}$. The resulting expressions are lengthy; for completeness we provide them in the Supplement (see S2). As an illustration we here show the PSD of $\mathcal{B}$,

[TABLE]

with

[TABLE]

As detailed in the Supplement (see Sec. S2) the power spectra of $S,I$ and $\mathfrak{X}$ can be expressed in terms of that of $\mathcal{B}$; many of the characteristics of the spectra of $S,I$ and $\mathfrak{X}$ are shared with those of $\mathcal{B}$, or directly related to it. We note that the RHS of Eq. (16) is proportional to $1/|g|^{2}$, and the same is the case for the spectral densities of $\mathfrak{X},S$ and $I$ (see Eqs. (S2)); as a result, some of the key properties of the power spectra are determined by the behaviour of $|g|^{2}$, as discussed in more detail below.

3.4 Test Against Simulations

To illustrate the model and test our analytical results, we sampled possible heterogeneous populations. Specifically, the simulations shown in Fig. 3 are for populations with five susceptible and three infected subclasses. For each example, the probabilities $\{p_{i}\}$ and $\{q_{a}\}$ were drawn at random from a flat distribution over the simplexes $\sum_{i}p_{i}=1$ and $\sum_{a}q_{a}=1$. Susceptibilities and infectivities were assigned randomly in the intervals $0.5\leq\chi_{i}\leq 2.5$ and $0.3\leq\beta_{a}\leq 1.3$. Simulations are for $N=10^{6}$, and the rates for recovery, birth/death and immigration were set at \rho=0.07\leavevmode\nobreak\, $\kappa=5.5\times 10^{-5}$ and $\xi=5\times 10^{-6}$ respectively. The rates $\beta_{a}$, $\rho$, $\kappa$ and $\xi$ have units of $\mbox{days}^{-1}$, whereas $\chi_{i}$ is dimensionless. The chosen rates are representative of childhood diseases such as whooping cough, measles, rubella or chickenpox [50].

The resulting PSDs are shown in Fig. 3. The continuous thick lines show the analytical result, and dashed lines are obtained from simulations, as an average over realizations of the individual-based model. As can be seen from the figure, the predictions of Eqs. (S2) precisely match the results from simulations. In all figures, axes labelled ‘frequency’ show $f=\omega/2\pi$, and have units of $\mbox{days}^{-1}$.

4 Consequences of Heterogeneity

Having established an analytical description of quasi-cycles, we now use this theory to identify which properties of the distribution of $p_{i}$, $\chi_{i}$, $q_{a}$ and $\beta_{a}$ are most relevant for the characteristics of stochastic quasi-cycles in heterogeneous populations. Specifically, we study how heterogeneity in the population affects the dominant frequency of quasi-cycles, their amplitude and the sharpness of the spectra. We will then also discuss if and how the different subgroups synchronise during the epidemic cycles.

4.1 Dominant Cycle Frequency

Numerical inspection of the different terms in the analytical solution of the PSDs suggests that the dominating element is the factor $1/|g|^{2}$, as briefly indicated in Sec. 3.3. The frequency for which $|g|^{2}$ reaches its minimum roughly corresponds to the dominant cycle frequency, $\omega_{d}$, in the PSDs. The minimum of $|g|^{2}$ can be found by differentiation of the expression in Eq. (17). Assuming that $C\chi_{i}\gg\kappa$ we further approximate the location of this minimum. This assumption is valid if infection processes occur on a time scale which is much shorter than the life expectancy of an individual. Further, we assume that $\omega\gg A_{i}$, i.e. that a susceptible individual typically lives through several epidemic events before it becomes infected. Both approximations are intuitively plausible for childhood diseases, known to show periodic outbreaks [50]. Making these assumptions we find that the frequency for which $|g|^{2}$ is minimal can be approximated as

[TABLE]

This implies that the characteristic frequency is determined (mostly) by the mean susceptibility at birth and the mean infectivity at infection ($\overline{\chi}$ and $\overline{\beta}$) and the capacity of replenishment of the susceptible pool ($\kappa$).

The validity of our approach is confirmed in Fig. 4(a), where we test the approximation against simulations for a wide set of parameters. A perhaps more intuitive representation of our result can be found in Fig. 4(b), where we show the power spectra of several sample populations, each with different distributions of $\{p_{i},\chi_{i},q_{a},\beta_{a}\}$, but all with the same first moments $\overline{\chi}$ and $\overline{\beta}$. As seen in the figure, this produces spectra of different amplitudes but with the same characteristic frequency.

4.2 Amplitude of Stochastic Cycles

While we have found above that the dominant frequency of stochastic cycles is largely determined by the first moments $\overline{\chi}$ and $\overline{\beta}$, the results shown in Fig. 4(b) demonstrate that this is not the case for the amplitude of the spectra at the dominant frequency. To investigate this further we evaluate the analytic expressions for the PSDs in Eqs. (S2) at the approximation of $\omega_{d}$ in Eq. (18). Making the same assumptions as in Section 4.1, we find that the height of the peak in the power spectra can be approximated as

[TABLE]

We note the presence of the second moments $\overline{\chi^{2}}$ and $\overline{\beta^{2}}$, unlike in Eq. (18). This indicates that the spread of susceptibilities and infectivities is relevant to the size of the epidemic.

In Fig. 5 we plot results from the approximation in Eqs. (19) against the maximum amplitude of spectra obtained numerically from the full expression (within the LNA), see Eqs. (S2) in the Supplement. The data confirms that the approximation is valid for a wide range of parameters. While we find slight deviations at large amplitudes in the case of the infectives, the approximation is very robust for the susceptible population.

4.3 Sharpness of the Spectra

We now turn to the sharpness of the peak in the PSDs. The sharper the peak, the closer the stochastic outbreaks are to perfect cyclic behaviour. Conversely, cyclic behaviour is less distinct if the peak in the spectrum is shallow. This has been described before as the ‘coherence’ of the spectra [13]. As we will investigate a different notion of coherence in Sec. 4.4 and in order to avoid confusion, we will refer to the concentration of power near the peak of the spectrum as ‘sharpness’.

Following [13], we define the sharpness as the relative spectral power accumulated in an interval around the peak,

[TABLE]

We compute the sharpness numerically, using the expressions in Eqs. (S2). In order to evaluate the denominator in Eq. (20) we integrate up to an upper cutoff of $\omega_{max}=\pi/100\leavevmode\nobreak\ \mbox{days}^{-1}$. In the numerator we use $\Delta\omega=0.05\,\omega_{max}$. The choice of $\Delta\omega$ can be illustrated using Fig. 4(b), where the sharpness $\mathbb{S}$ of the peak roughly corresponds to the fraction of total power concentrated in the interval between frequencies of $0.0015$ and $0.002$ days*-1*.

In Fig. 6 we show the sharpness of spectra for $200$ random populations (as described in Fig 4). It is clear from the figure that there is a trend of increasing sharpness as the product of the mean susceptibility and infectivity at birth approaches unity (in the dimensions used here). The spread of the markers on the vertical axis indicates that there are significant effects of heterogeneity. It proves difficult, though, to find a functional dependence on higher moments of the distributions of susceptibilities and/or infectivities which would further collapse the data.

4.4 Synchronization between Subgroups

We have established so far that introducing heterogeneity leads to significant changes in the quasi-cycles of the aggregate numbers of susceptible and infective individuals. However, we have not yet said much about the dynamics of the individual subgroups. In Fig. 7 we show the same example of sustained oscillations as in the inset of Fig. 2, but instead of the total susceptible and infected population we now highlight the time evolution of each of the subgroups.

In the upper two panels, (a) and (b), we show time series of the number of individuals in each subgroup normalised by the total population size. More specifically, we show susceptible subclasses ($n_{i}/N$) in panel (a), and infective subclasses ($m_{a}/N$) in panel (b). For each of these, stochastic oscillations can be observed. These cycles are pronounced for the case of the infective subgroups, panel (b), and more shallow for the susceptibles, panel (a). This is to be expected, given that the total number of susceptibles is more than an order of magnitude larger than those of the infectives (see also Fig. 2). From Fig. 7 (a) and (b) it is clear that all subgroups undergo cycling of roughly the same frequency. This is confirmed by the power spectra in Fig. 8.

We note that these statements rely on expressing number of individuals in each class as a fraction of the total population, and not relative to the time-dependent total number of susceptibles or infectives respectively. We contrast the above with a representation in which we express the occupancy in each infective subgroup as a fraction of the infectives only, and similarly for the susceptibles. To this end we replot the simulation run shown in Fig. 7 (a) and (b), but now in terms of $n_{i}/(NS)$ and $m_{a}/(NI)$, respectively. The quantities $NS=\sum_{j}n_{j}$ and $NI=\sum_{b}m_{b}$ are the total number susceptible and infective individuals respectively, and they are time-dependent themselves. Results are shown in Fig. 7(c) and (d). Although the overall number of infectives, $NI$, undergoes the noise-driven cycles shown in Fig. 2, we find no discernible structure within the group of infectives; the time series $m_{a}/(NI)$ in Fig. 7(d) are essentially flat noisy lines. This is what one would expect, since the allocation to each subgroup, $I_{a}$, of infectives is random when an individual is newly infected, and the recovery rate is the same for all infective subgroups.

A more complex behaviour can be seen within the group of susceptibles. This group as a whole undergoes stochastic cycles (see Fig. 2), but an interesting structure is observed within the group of susceptibles as well. The time series $n_{i}/(NS)$ in Fig. 7(c) show cyclic behaviour, and – to a good approximation – any pair of these time series is either in phase with each other, or they have a phase difference of $\pm\pi$.

To explore the phase lag between the different time series we use the so-called complex coherence function [51]. This technique relies on computing the cross-spectrum $\left\langle{\widehat{x}_{i}(\omega)\widehat{x}_{j}^{*}(\omega)}\right\rangle$ between time series $x_{i}(t)$ and $x_{j}(t)$. The phase lag is then obtained as

[TABLE]

We stress that the subscript ∗ denotes complex conjugation, and is not to be confused with ⋆, used earlier to indicate fixed points of the deterministic dynamics. Eq. (21) returns a phase lag for each spectral component, $\omega$. Details can be found in the Supplement (see S3).

The phase lag between the different groups of susceptible individuals is shown in Fig. 9. The data in panel (a) corresponds to Fig. 7 (a). More precisely, in Fig. 9 (a) we pick the time series $n_{1}/N$ as a reference, and show the phase lag of all subgroups $n_{i}/N$ with respect to this reference time series. We find that the phase lag for frequencies around the dominant frequency in the power spectra is small, consistent with Fig. 7 (a); all time series $n_{i}/N$ oscillate (roughly) in phase with each other. In Fig. 9 (b) we repeat this procedure, but now taking the time series $n_{i}/(NS)$ as an input, corresponding to Fig. 7 (c). One then finds a rather different picture; the phase lag around the dominant frequency takes values either near zero, or close to $\pm\pi$. This indicates that the different classes of susceptible individuals fall into two groups. The time series in either group are in phase with each other, and in anti-phase with those in the respective other group. A closer inspection shows that these two groups are formed by the time series $i$ with $x_{i}^{\star}<S^{\star}/K$ and with $x_{i}^{\star}>S^{\star}/K$ respectively. This behaviour in turn can be understood intuitively by revisiting Eqs. (9). Assuming $\kappa\ll(\xi+\mathcal{B}^{\star})\chi_{i}$ for all $i$ (a valid approximation for the cases analysed here), we find $x_{i}^{\star}\propto 1/\chi_{i}$, indicating that the more susceptible classes are less populated at the deterministic fixed point than the less susceptible ones. During the increasing leg of a stochastic cycle, we expect the number of newly infected individuals among class $i$ to be proportional to $x_{i}^{\star}\chi_{i}$, suggesting that all susceptible classes are depleted in equal absolute numbers. This in turn means that subclasses with $x_{i}^{\star}>S^{\star}/K$ will represent an even larger fraction of the susceptible population as the total susceptible population decreases, while the subclasses with $x_{i}^{\star}<S^{\star}/K$ will represent a smaller fraction. This is what is observed in Fig. 7 (c).

5 Conclusions

In summary, we have explored the SIR model in finite populations, including demographic processes and allowed for agent-to-agent heterogeneity in both the susceptibility to a disease and the capacity to spread the disease. This system combines the effects of intrinsic demographic stochasticity (due to random infection, recovery and birth-death events), with quenched heterogeneity. The focus of our paper is to characterise the interplay between these two types of stochasticity, and to investigate how the heterogeneity between individuals affects quasi-cycles driven by intrinsic noise. Our analysis relies on the system-size expansion, which allows us to compute the properties of these cycles analytically in the linear-noise approximation.

Our principal results can be summarised as follows: (i) In the deterministic limit of infinite populations no closed set of equations for macroscopic quantities can be found in the transient regime. Fixed points for aggregate quantities of this deterministic dynamics can however be fully determined from a set of closed equations for the total susceptible ($S^{\star}$) and infected ($I^{\star}$) population, and weighted averages of the susceptibility ($\mathfrak{X}^{\star}$) and infectivity ($\mathcal{B}^{\star}$). (ii) Similarly, the Langevin equations in the linear-noise approximation do not close easily at the aggregate level, but a closed set of equations for the spectra of fluctuations in $S,I,\mathfrak{X}$ and $\mathcal{B}$ about the deterministic fixed point can be found in the stationary state. These can be used to analytically describe the stochastic oscillations about the fixed point. (iii) Within reasonable assumptions, the characteristic frequency of the noise-driven oscillations is determined mostly by the mean susceptibility and infectivity at birth or infection ($\overline{\chi}$ and $\overline{\beta}$). However, the amplitude of the oscillations and the sharpness of peaks in the power spectra will generally depend on the higher moments of the distribution of susceptibilities and infectivities,in particular also on the agent-to-agent heterogeneity. (iv) Finally, the number of individuals in the different subclasses of infectives and susceptibles undergo stochastic cycles as well. If expressed in relation to the total population, these time series are synchronised and in phase. Normalized against the time-dependent total number of infectives, however, the different infective subclasses show no discernible oscillatory behaviour. Using a similar normalization within the susceptible population, we find that different subclasses are syncronized and either in phase with each other or have a phase difference of $\pm\pi$. These results are confirmed analytically. Regardless of the normalization, we find that the periodic outbreaks do not follow a hierarchical infection process, and all subgroups have similar absolute depletion/increase in absolute numbers. This is in contrast to what has been reported in single outbreak studies [35, 39]. However, it is important to note that in this existing work the outbreak is tracked in an initial transient period. Our results are valid after this period, at a deterministic fixed point, where the susceptible population is distributed in inverse proportion to their susceptibility (as explained above); this is a scenario different to the one studied in [35, 39].

We think our results can be relevant for future work in several ways. First, our work contributes to the ongoing discussion about when and how a model with heterogeneity can be replaced or approximated by a homogeneous model. In previous studies, heterogeneous models were compared to homogeneous models with susceptibility equivalent to the arithmetic [52] or harmonic mean [41] of the susceptibilities in the different groups. More recently, the focus has been placed on equivalent basic reproduction numbers ($R_{0}$) [53]. In the heterogeneous model this requires estimating $R_{0}$ based on, for example, the outbreak size, and therefore the comparison is not straightforward. Here we have shown that all models within the class we have looked at and with equal values of $\overline{\chi}\overline{\beta}$ generate periodic outbreaks with the same dominating frequency. This characteristic frequency can be used to define a unique homogeneous model to which models of varying degrees of heterogeneity can be compared. Furthermore, the dependence of the spectra of oscillations on both the first and higher moments of the distribution of heterogeneity might provide an avenue towards estimating how heterogeneous a population is from the observation of epidemic cycles. Finally, the formalism we have developed is versatile and can be applied to study quasi-cycles in other areas in which heterogeneity might be relevant, for example in predator-prey dynamics or evolution [17, 54, 55, 19, 56, 57]. Our findings indicate that the frequency of quasi-cycles can, to a good approximation, be obtained from the first moment of the distribution of heterogeneous agent properties, but that their amplitude depends on higher moments of the disorder. We expect similar behaviour in other heterogeneous systems with noise-driven cycles.

Acknowledgements

FHA thanks Consejo Nacional de Ciencia y Tecnología (CONACyT, Mexico) for support. TG acknowledges funding by the Engineering and Physical Sciences Research Council (EPSRC, UK) under grant number EP/K037145/1. TG would like to thank the Group of Nonlinear Physics, University of Santiago de Compostela, Spain for hospitality.

Author contributions statement

FHA and TG conceived and designed the study. FHA carried out the analytical calculations and computer simulations. FHA and TG interpreted results. FHA and TG wrote and reviewed the manuscript.

Additional information

The authors have no competing financial interests.

**The effects of heterogeneity

on stochastic cycles in epidemics

 

Supplemental Material

** Francisco Herrerías-Azcué and Tobias Galla

francisco.herreriasazcue@postgrad.manchester.ac.uk, tobias.galla@manchester.ac.uk

Theoretical Physics, School of Physics and Astronomy,

The University of Manchester, Manchester M13 9PL, United Kingdom

S1 Linear-noise approximation

Carrying out the system-size expansion for the model with heterogeneity is tedious, but straightforward and follows the lines of [43]. The final outcome is the linear-noise approximation in Eqs. (3.1). The variables $\eta_{i}$ and $\nu_{a}$, represent Gaussian noise, with no correlation in time, but with potential correlation between the different noise variables at equal time. These noise variables can be decomposed as

[TABLE]

where, broadly speaking, each term on the right-hand side represents one possible type of event in the microscopic model. For example, $u_{ia}$ relates to spontaneous infection of a susceptible individual of type $S_{i}$, resulting in a newly infective of type $I_{a}$. Similarly, $v_{iab}$ represents an event in which an individual of type $S_{i}$ is infected by an individual of type $I_{a}$, and the newly infected is of type $I_{b}$. The variable $w_{a}$ relates to a recovery event of an individual of type $I_{a}$, death of susceptible $S_{i}$ and simultaneous birth of susceptible $S_{k}$ is reflected by $x_{ik}$; death of an individual of type $I_{a}$ and simultaenous birth of susceptible $S_{i}$ is described by $y_{ai}$, and finally death of a recovered individual and simultaneous birth of susceptible $S_{i}$, by $z_{i}$. The signs on the right-hand-side in Eqs. (S1) reflect the fact that each of these events may either increase or reduce the number of individuals of type $S_{i}$ and $I_{a}$, respectively.

Each of the noise variables on the right-hand-side of Eqs. (S1) are uncorrelated in time, and they have no cross-correlations. Within the LNA their variances are set by the corresponding reaction rates at the deterministic fixed point, i.e. we have

[TABLE]

Using the shorthand introduced in Eqs. (3.2), we then find

[TABLE]

which are needed for the computation of the PSDs.

S2 Calculation of power spectra

We start from the result in Eqs. (3.2) in Sect. 3.2:

[TABLE]

As an illustration let us now compute the power spectrum of $\widehat{\mathcal{B}}$. To keep equations manageable, we define

[TABLE]

and so we write the Fourier transform of $\tilde{\mathcal{B}}$ as

[TABLE]

where $f_{i}$, $\beta_{a}$, $\widehat{\eta}_{i}$ and $\widehat{\nu}_{a}$ are all functions of $\omega$. We then find

[TABLE]

The notation $*$ denotes complex conjugation. Substituting the noise correlators from Eqs. (S1),

[TABLE]

and, using Eq. (S6), we find

[TABLE]

which is the PSD of $\mathcal{B}$, as also reported in Eq. (16) in the main text.

Following the same process, we can compute the PSD for the remaining quantities, $\mathfrak{X}$, $I$ and $S$. We do not report all details, but only the final results

[TABLE]

The power spectra of fluctuations for the individual subgroups of infectives and susceptibles are found as

[TABLE]

S3 Phase Lag

In order to explore the the phase lag we use the so-called complex coherence function, $\mathcal{CCF}_{ij}$, between subgroups $i$ and $j$, defined as

[TABLE]

where $\widehat{x}_{i}$ and $\mathcal{P}$ are functions of $\omega$.

For $i\neq j$ this is in general a complex-valued function (of $\omega$). The argument of $\mathcal{CCF}_{ij}$, given by

[TABLE]

is known as the phase spectrum; it describes the phase-lag between the time series $x_{i}(t)$ and $x_{j}(t)$ [58].

The cross spectra of the population in the susceptible classes normalized with respect to the total population ($x_{i}=n_{i}/N$) is given by

[TABLE]

This can be written as

[TABLE]

where we introduced the notation

[TABLE]

From these we obtain the phase lag as

[TABLE]

which yields the theoretical lines in Fig. 9a.

To explore the phase lag between the susceptible subgroups when normalized by the total susceptible population ($x_{i}^{\prime}=n_{i}/NS$), we first need to compute the cross-spectra of the renormalized signals $\mathcal{P}_{x_{i}^{\prime}x_{j}^{\prime}}(\omega)$. As in Section 3.1, we start from the ansatz

[TABLE]

We then have

[TABLE]

and so (after expanding in $1/\sqrt{N}$)

[TABLE]

In Fourier space this turns into

[TABLE]

For the cross spectra we then find

[TABLE]

which can be rewritten as

[TABLE]

We have introduced the notation $Y_{ij\,R}=\mbox{Re}\left[Y_{ij}\right]$ and $Y_{ij\,I}=\mbox{Im}\left[Y_{ij}\right]$ with

[TABLE]

From these, we can find the phase-lag as

[TABLE]

This expression was used to obtain the analytical predictions shown in Fig. 9b.