An SIS epidemic model with individual variation

1.   Introduction

The SIS (Susceptible-Infectious-Susceptible) Model was introduced by Weiss and Dishon ^[1] to study infections in a closed population of $n$ individuals, infections that do not confer any long lasting immunity (gonorrhea, or the common cold, for example). If $Y(t)$ is the number of infectives at time $t$, then $(Y(t), \, t\geq 0)$ is a continuous-time Markov chain taking values in $\{0, 1, \dots, n\}$ with transitions

The large-$n$ behaviour of this model is well understood, because the proportion of infectives $Y(t)/n$ obeys the following law of large numbers, which is an application of Theorem 3 of ^[2]; to the best of my knowledge, this result has not been stated explicitly elsewhere.

Theorem 1. If $Y(0)/n \to y_0$ as $n\to\infty$, then $(Y(t)/n)$ converges in probability, uniformly over finite time intervals, to the solution of the ODE

where $\rho = \mu/\lambda$, namely

We see immediately that if $\lambda > \mu$ then (1.1) has a globally stable equilibrium at $y^* = 1-\rho$, while if $\lambda < \mu$ the disease-free state $0$ is globally stable. The two behaviours are illustrated in Figure 1: quasi stationarity (the infection becomes endemic) versus evanescence (the infection dies out). Notice how the sample paths of the process $(Y(t)/n)$ (blue) "track" the limiting deterministic path $y(t)$ (green).

Several elaborations of the basic SIS model have been proposed where the assumption of homogeneity among individuals is relaxed. Perhaps the most is extensive is the study of epidemics where transmission occurs within households ^{[3,4,5,6,7,8]}. Another ^[9] assumes that individuals encounter one another in a communal meeting place. We incorporate heterogeneity by assuming that individuals have their own, possibly different, infection and recovery rates. In particular, we assume that individual $i$ $(i = 1, \dots, n)$ has an exponentially distributed recovery period with mean $\mu_i^{-1}$ and a resistance level $\lambda_i^{-1}$ (both finite). Let $X^{_{(n)}}_i$ be $1$ or $0$ according to whether individual $i$ is infected or not, and let $X^{_{(n)}}_{\phantom{i}} = (X^{_{(n)}}_1\!, \dots, X^{_{(n)}}_n)$ be the state of the population. Suppose that $(X^{_{(n)}}_{\phantom{i}}(t), \, t\geq 0)$ is a continuous-time Markov chain taking values in $\{0, 1\}^n$ with transitions

where $\bar X^{_{(n)}}_{\phantom{i}} = \frac{1}{n} \sum_{j = 1}^n X^{_{(n)}}_j$ is the proportion of the population that is infected. So we assume here that each infected individual makes the same contribution to the infection potential of the group of infectives. The model in described in ^[9] has $\mu_i = 0$ and $\lambda_i$ to be interpreted as the proportion of time individual $i$ spends in the communal meeting place. Daley et al. ^[13] also consider the case where $\mu_i = 0$ for all $i$ and study the mean duration of the epidemic.

Notice that the disease free state ${\boldsymbol{0}} = (0, 0, \dots, 0)$ is the sole absorbing state and, since $(\lambda_i)$ and $(\mu_i)$ are strictly positive, the remaining states form a communicating class from which ${\boldsymbol{0}}$ is accessible, and indeed reached with probability $1$. None-the-less, the process may exhibit quasi stationarity. By investigating the large-$n$ behaviour of the model, we can determine conditions for such a quasi equilibrium to be achieved. We will see that accounting for variation amongst individuals can lead to different conclusions about equilibrium behaviour to those reached for the standard SIS model. Several examples are given to show a range of behaviour. For example, the model presented here may predict quasi stationarity in cases where the standard model predicts evanescence.

2.   Large-$n$ behaviour

Because the population is now heterogeneous, a quite different approach is needed. We will follow McVinish and Pollett ^[10], who considered a setup that also encapsulates metapopulation models, and other population models. Think of the individual characteristics $\theta_i: = (\lambda_i, \mu_i)$ as (random) points in a subset $S$ of $\mathbb{R}_+^2$, and define sequences of random measures $(\sigma^{_{(n)}}_{\phantom{i}})$ and random-measure-valued processes $(m^{_{^{(n)}}}_t\!, \, t\geq 0)$ by

where $B$ is any Borel subset of $S$. We will suppose that $\sigma^{_{(n)}}_{\phantom{i}} \stackrel{d}{\to} \sigma$ (weak convergence) for some non-random (probability) measure $\sigma$, and then hope to establish that $m^{_{^{(n)}}}_t \stackrel{d}{\to} m^{_{^{(n)}}}_{\phantom{t}}$, so that properties of $X^{_{^{(n)}}}_t$ can be approximated for large $n$. In particular, if the individual characteristics were chosen from some non-negative bivariate distribution, then $\sigma^{_{(n)}}_{\phantom{i}} \stackrel{d}{\to} \sigma$ would follow as a consequence of the law of large numbers.

The sequences $(\sigma^{_{(n)}}_{\phantom{i}})$ and $(m^{_{^{(n)}}}_t)$ may be defined, equivalently, by

for $h$ in $C_{_{^{\!b}}}\!(S)$, the class of bounded continuous functions that map $S$ to $\mathbb{R}$. For example if $(h\equiv 1)$ then $m^{_{^{(n)}}}_t\!(S) = \int m^{_{^{(n)}}}_t\!(\mathrm{d}\theta) = \frac{1}{n}\sum_{i = 1}^{n} X^{_{^{(n)}}}_{i, t}$, the proportion of population that is infected at time $t$.

Our main result is a consequence of Theorem 1 of ^[10].

Theorem 2. Suppose that $\sigma^{_{(n)}}_{\phantom{i}} \stackrel{d}{\to} \sigma$ and $m^{_{^{(n)}}}_{_{^{0}}} \stackrel{d}{\to} m_{_{^{0}}}$ for some non-random measures $\sigma$ and $m_{_{^{0}}}$. Then, the sequence of measure-valued processes $(m^{_{^{(n)}}}_t\!, \, t\geq 0)$ converges weakly to the unique solution $(m_t, \, t\geq 0)$ of

where (notation) $(h, m) = \int h(\theta) m(\mathrm{d}\theta)$, and

Lemma 5 of ^[10] implies that $m_t(B)\leq \sigma(B)$, for all Borel sets $B$, and in particular that $m_t$ is absolutely continuous with respect to $\sigma$. It follows that $m_t$ has a (uniquely determined $\sigma$-a.e.) Radon-Nikodym derivative $\phi_t$ $(\geq 0)$ with respect to $\sigma$: $m_t(B) = \int_B \phi_t(\theta) \sigma(\mathrm{d}\theta)$. And further it implies that $\phi_t\leq 1$.

Now we may differentiate both sides of (2.2) with respect to $\sigma$ to show that $\phi_t(\lambda, \mu)$ must satisfy

or equivalently,

Note that $\phi_t(\theta)$ can be interpreted as the limiting probability that an individual with characteristics $\theta = (\lambda, \mu)$ is infected at time $t$.

3.   Long-term behaviour

Equation (2.3) can be used to study the long-term $(t\to\infty)$ behaviour of our model. For example, any equilibrium $\phi_{{{\rm{eq}}}}$ of (2.3) must satisfy $0 = \lambda (1- \phi_{{{\rm{eq}}}}) \int \phi_{{{\rm{eq}}}}(\theta) \sigma(\mathrm{d}\theta) -\mu \phi_{{{\rm{eq}}}}$. On setting $\psi = \int \phi_{{{\rm{eq}}}}(\theta) \sigma(\mathrm{d}\theta)$, we see that

and so, on integrating this over $(\lambda, \mu)\in S$, we find that $\psi$ must solve the equation

The following result characterizes the large-$n$ equilibrium behaviour of our model in terms of $r_0: = R^{ {\, \prime}}(0+)$ ($\leq \infty$). It follows from Theorems 2, 3 and 4 of ^[10], noting that their Condition D$^{ {\, \prime}}$ holds in the present setting. Note also that since the integrand in (3.1) is less that $1$, and $\sigma$ is a probability measure, we can differentiate under the integral to obtain

which implies that $r_0 = \int (\lambda/\mu) \, \sigma(\mathrm{d}\lambda, \mathrm{d}\mu)$.

Theorem 3. (a) If $r_0\leq 1$, then $\psi = 0$ is the only fixed point of $R$, and $\phi_{{{\rm{eq}}}} = 0$ is globally stable, that is, for all $\phi_{_{^{0}}}$, $\phi_t\to 0$ on $S$. The latter entails $m_t(B)\to 0$, for all $B$, and hence the disease free state is globally stable.

(b) If $r_0 > 1$, then $R$ has two fixed points, $0$ and a positive fixed point $\psi_*$, and if $\left(m_0(S) = \right)$ $(\phi_{_{^{0}}}, \sigma) > 0$, then

The latter entails $m_t(B)\to m_*(B)$, for all $B$, where

in particular $m_*(S) = \psi_*$, thus implying quasi stationarity.

Example 1. Suppose that the individual characteristics are chosen independently from $\Gamma(a_\lambda, r_\lambda)$ and $\Gamma(a_\mu, r_\mu)$ distributions, respectively. Notation: $r_\lambda$ and $r_\mu$ are rate parameters. We have that

Figure 2 illustrates (a) quasi-stationary behaviour, and (b) evanescence. In (a), the green line is at $\psi_* = 0.8035$, very close to $0.8$ (dashed red), which is the equilibrium of the homogeneous model using $\lambda = 10$ and $\mu = 2$. Figure 3 also illustrates quasi-stationary behaviour, now with $\psi_* = 0.2435$, but the corresponding homogeneous model ($\lambda = 1.4$ and $\mu = 1.5$) predicts evanescence. This disparity between the two models happens when

Note that $\psi_*$ was evaluated by way of fixed point iteration of (3.1) using Matlab's ${\tt integral2}$, which evaluates double integrals numerically.

Example 2. Suppose that the individual characteristics are chosen from a bivariate Gamma distribution of the form given in Theorem 2 of ^[11]. The marginals are Gamma distributions with a common rate parameter, which we denote by $r$, and shape parameters $a_\lambda$ and $a_\mu$, respectively, and correlation $\rho = \sqrt{a_\lambda a_\mu}/c$, where $c > 0$. Using the representation Eq.3 on Page 768 of ^[11], we find that

The simplest case has $r = 1$, and $b_\lambda = b_\mu = 1$ (which necessitates $a_\lambda = a_\mu = c-1$, and $c > 1$), and $\sigma(\mathrm{d}\lambda, \mathrm{d}\mu) = f(\lambda, \mu)\, d\lambda\, \mathrm{d}\mu$, with

where $C^{-1} = \Gamma(c)/(c-1)^2$, and $\Gamma(a, x) = \int_x^\infty t^{a-1} \mathrm{e}^{-t}\, { \mathrm{d}}t$ is the complementary incomplete gamma function. Figure 4 illustrates quasi-stationary behaviour with $\psi_* = 0.1941$, noting the corresponding homogeneous model ($\lambda = \mu = 1.5$) predicts evanescence.

4.   Moments

We can get a handle on various statistics associated with the process by way of the moment generating function (MGF) of $m_t$:

which will be defined for $(a, b)$ in a region containing the origin. On setting $h(\theta)$ ($= h(\lambda, \mu)$) $= \exp(a\lambda + b\mu)$ in (2.2) we obtain an integral equation for $M_t$:

where $\Sigma(a, b) = \int \mathrm{e}^{a\lambda + b\mu} \sigma(d\lambda, d\mu)$ is the MGF of $\sigma$. This can be differentiated to obtain the PDE

where

On setting $({\partial}/{\partial t}) M_t(a, b) = 0$, we see that any equilibrium $M_{{{\rm{eq}}}}$ of (4.1) must satisfy

where $\psi = M_{{{\rm{eq}}}}(0, 0)$. This PDE is easily solved using the method of characteristics. We find that

Under the conditions of Theorem 3 (b), namely that $r_0 > 1$ and $\left(m_0(S) = \right)$ $(\phi_{_{^{0}}}, \sigma) > 0$, we will have that $M_t(a, b)\to M_{{{\rm{eq}}}}(a, b)$, $M_{{{\rm{eq}}}}$ being the MGF of $m_*$:

We will necessarily have $M_{{{\rm{eq}}}}(0, 0) = \psi_*$, and, since $R(\psi_*) = \psi_*$,

Example 3. Suppose that the individual characteristics are chosen from independent exponential distributions with rate parameters $l$ and $m$: $\sigma(\mathrm{d}\lambda, \mathrm{d}\mu) = l \mathrm{e}^{-l\lambda} \, m \mathrm{e}^{-m\mu} \, \mathrm{d}\lambda\, \mathrm{d}\mu$. With this choice we can perform some explicit calculations. For example,

where $\rho = l/m$. Therefore,

So, we have $R^{ {\, \prime}}(0+) = \infty$. This implies a stable equilibrium. In particular, $R$ has two fixed points, $0$ and a positive fixed point $\psi_*$. Again, a simple calculation yields

Notice that $M_{{{\rm{eq}}}}(0, 0) = R(\psi_*) = \psi_*$. For example,

a quantity which approximates the large-$n$ long-term behaviour of (refer to (2.1))

being the average recovery rate of infected individuals. This is illustrated in Figure 5, where the average recovery rate of infected individuals (cyan) is plotted along with the value of (4.3) (yellow).

5.   Conclusion

We have presented a model for an SIS infection that incorporates variation of individual characteristics. By imagining that these characteristics follow a point process on a subset of $\mathbb{R}_+^2$, we are able to exploit existing results on random measures, and random measure-valued processes, to obtain large-population limits, and then predict long-term behaviour. These results provide a means of studying various statistics associated with the process, including the proportion of individuals infected. Several examples demonstrate a range of behaviour that cannot be gleaned from the standard stochastic SIS model based on mean values of the characteristics. Explicit results on the moment generating function of the (large-$n$) limiting process allow us to estimate additional quantities such as the average infection rates and average recovery rates.

The methods described here do not provide explicit results on the quality of the approximations presented, and in particular how large $n$ would need to be for our approximations to be faithful. It might be possible to adapt the results in Chapter 4 of ^[12], based on their discrete-time counterparts ^[13], to provide explicit rates of convergence by attempting to bound weighted sums of differences:

where $p_t = (p_{i, t}, \ i = 1\dots n)$ is some suitable deterministic law of motion; this is the subject of ongoing research. A simple extension to the model presented would be to replace the upward transition rate in (1.2) by $a+\lambda_i \bar X^{_{(n)}}_{\phantom{i}}$, where $a > 0$, thus incorporating an external source of infection as in ^[14]. Now the disease-free state ${\boldsymbol{0}}$ would no longer be an absorbing state, but the methods employed here could be brought to bear, thus providing a means of estimating equilibrium behaviour.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare there is no conflict of interest.