A Model for Reinfections and the Transition of Epidemics

Abstract

Reinfections of infected individuals during a viral epidemic contribute to the continuation of the infection for longer periods of time. In an epidemic, contagion starts with an infection wave, initially growing exponentially fast until it reaches a maximum number of infections, following which it wanes towards an equilibrium state of zero infections, assuming that no new variants have emerged. If reinfections are allowed, multiple such infection waves might occur, and the asymptotic equilibrium state is one in which infection rates are not negligible. This paper analyzes such situations by expanding the traditional SIR model to include two new dimensionless parameters, ε and θ, characterizing, respectively, the kinetics of reinfection and a delay time, after which reinfection commences. We find that depending on these parameter values, three different asymptotic regimes develop. For relatively small θ, two of the regimes are asymptotically stable steady states, approached either monotonically, at larger ε (corresponding to a stable node), or as waves of exponentially decaying amplitude and constant frequency, at smaller ε (corresponding to a spiral). For θ values larger than a critical, the asymptotic state is a periodic pattern of constant frequency. However, when ε is sufficiently small, the asymptotic state is a wave. We delineate these regimes and analyze the dependence of the corresponding population fractions (susceptible, infected and recovered) on the two parameters ε and θ and on the reproduction number R₀. The results provide insights into the evolution of contagion when reinfection and the waning of immunity are taken into consideration. A related byproduct is the finding that the conventional SIR model is singular at large times, hence the specific quantitative estimate for herd immunity it predicts will likely not materialize.

1. Introduction

The COVID-19 pandemic has created tremendous interest from a plethora of investigators to better understand the onset of infection waves, their propagation in time and space, and the effect of the multitude of its parameters, both physiological and behavioral. The literature abounds with such studies. Many rely on the standard SIR (Susceptible, Infected, and Recovered) model [1] and its variations. Stochastic analogs [2,3] and other extensions, e.g., those that interpret infection data using statistical techniques [4,5,6], have appeared. Here, we cited only very few examples of such studies, recognizing that many important related works have been published.

In a recent paper [7], we presented an approach to model the SIR framework by drawing analogies with chemical reaction engineering processes. Based on areal densities (namely population per unit area) rather than absolute population numbers, this approach allows for a more correct formulation characterized in terms of dimensionless parameters that express chemical and biological kinetics and also mobility and transport mechanisms (advection and diffusion). Interestingly, in the absence of spatial gradients, which is the underlying assumption of the standard SIR model, one finds that infection is governed by only one parameter, the so-called reproduction number, $R_0=\frac{L}{T}$, expressed in [7] in terms of physiological and behavioral properties. Here, $T$ is the average characteristic time parametrizing infection and depends on both physiological and behavioral parameters, while $L$ is the average characteristic time for the recovery of infected individuals (e.g., assumed 14 days, for the case of COVID-19). It was shown in [7] that in the SIR framework under the assumptions that $R_0$ remains constant, that no new variants emerge, and that there are no reinfections, the process asymptotically wanes towards an equilibrium state of zero infections. This state is identified as denoting the onset of herd immunity, since any new infections will not lead to a new infection wave, but to an immediate decay. In the standard SIR model this state is asymptotically approached in an exponential manner. At the same time, we also showed [8], that the approach to the state of herd immunity can be slower—for example, following a power law (rather than an exponential decay) in time, if $R_0$ is not constant, but increases at late times, following a power-law dependence on the infected fraction. This is indeed likely to occur, at the later stages of infection, given the societal tendency to relax behavior as the infection wanes. Under this assumption, the resulting herd immunity also differs from that predicted by the traditional SIR model of constant $R_0$.

The standard infection models, including those in Ref. [7], do not consider reinfection. When the latter is allowed, however, a steady state of zero infections will not be reached at large times, whether exponentially or algebraically. Instead, epidemics will evolve into different, endemic states, in which infections continue. The nature of these states, how they are reached, and how they depend on the various process parameters form the objectives of this paper. We examine whether, in the case of reinfections, there are multiple, rather than single, infection waves, what their frequency is, and what are the properties of the resulting fractions of susceptible, infected, and recovered individuals.

Reinfection has been discussed in the general literature before, using a number of different approaches. We will cite a few. Thieme and Yang [9] and Song et al. [10] proposed generic SIR models with reinfection occurring without delay. Safan et al. [11] concentrated on a two-stage model for pertussis. Ashraf [12] summarized aspects of a SIR model that involves birth rate, vaccination, and waning immunity. Okhuese [13] and McMahon and Robb [14] predicted reinfection rates for COVID-19 using a more elaborate SEIRU model. Camacho and Vargas-De-Leon [15] developed a model for disease-free equilibria, with the stability of these equilibria controlled by the reproductive number and the recurrence rate. The possibility of delay times before reinfection sets in has also been included. Finally, Cooke [16] proposed general differential-difference models for diseases, while Beretta and Takeushi [17] and Enatsu et al. [18] generalized an SIR model by also including a delay time and studied the stability of the resulting state. Mena-Lorka and Hethcote [19] used a more elaborate SIR model to describe disease equilibria but did not endeavor in reinfection processes.

The review of the literature shows that we do not have at present a definitive approach to how to account for reinfections. Some works have proposed delay-type models, which lead to differential-difference equations. Others have postulated generic integral generation terms that lead to integrodifferential equations. This paper sheds light on this issue by adopting the reaction process-like contagion model proposed in ref. [7], here extended to include reinfections through the introduction of a novel rate of generation of new susceptible individuals. The model includes two new parameters, $\epsilon$ and $\theta$, denoting respectively, the kinetics of reinfections, and a delay time for the onset of reinfections. In general, the resulting description is in terms of a single non-linear integrodifferential equation, which, however, in the limit of fast reinfection kinetics (large $\epsilon$), reduces to a differential-difference equation, thus unifying the previous two different approaches in the literature. The paper extends the problem formulation, delineates the asymptotic regimes obtained, and examines the dependence of the various properties and patterns of infected, recovered, and susceptible population fractions on the three parameters, $R_0$, $\epsilon$, and $\theta$.

2. Mathematical Formulation and Results

We follow the approach in ref. [7], where contagion and recovery were represented by two equivalent chemical reactions, namely:

$$\mathcal{S}+\mathcal{J}\to \mathcal{J}+\mathcal{J}$$

(1)

$$\mathcal{J}\to \mathcal{R}$$

(2)

where $\mathcal{S}$, $\mathcal{J}$ and $\mathcal{R}$ denote susceptible, infected, and recovered individuals. To capture reinfections, we will propose an additional, equivalent chemical reaction that converts recovered populations to susceptible ones, namely:

$$\mathcal{R}\to \mathcal{S}$$

(3)

thus, leading to new susceptible individuals for infection. Reaction (3) will be assumed to be first order, with kinetics that account for various biological and behavioral factors. For the sake of generality, we will also allow the onset of the reaction to occur only after a characteristic delay in time has elapsed. Then, a recovered individual will have a finite probability of being reinfected, possibly multiple times. As in ref. [7], we will make significant simplifying assumptions, namely that there are no perished individuals, while we do not account for vaccinations that can confer additional immunity. These are made in order to focus on intrinsic, key aspects of reinfection. In short, we make two key assumptions: (1) A recovered individual can be reinfected only after a (dimensionless) time $\theta$, has elapsed. This time is the same for each recovered individual. (2) The rate by which recovered individuals are reinfected is proportional to the recovered population fraction—namely, a recovered individual has the same probability of being reinfected as any other infected one. Additional improvements can and should be added later as desired.

We proceed in the absence of spatial gradients (corresponding to what is also known as a “batch reactor” setting [7]) and a constant reproduction number $R_0$. The addition of chemical reaction (3) leads to a generalized SIR model containing an additional source term. Appendix A shows that this is described by a convolution integral. This leads to the following set of dimensionless equations for the three population fractions (susceptible, infected, and recovered):

$$s^{\prime }\left(t\right)=-R_{0}si+\epsilon {{\displaystyle \int }}_0^{t-\theta }i\left(\tau \right)\exp \left(-\epsilon \left(t-\theta -\tau \right)\right)d\tau H\left(t-\theta \right)$$

(4)

$$i^{\prime }\left(t\right)=R_{0}si-i$$

(5)

$$r^{\prime }\left(t\right)=i-\epsilon {{\displaystyle \int }}_0^{t-\theta }i\left(\tau \right)\exp \left(-\epsilon \left(t-\theta -\tau \right)\right)d\tau H\left(t-\theta \right)$$

(6)

Here, superscript indicates time derivative, and we introduced the two new dimensionless parameters, $\epsilon$ and $\theta$. In the absence of reinfections ($\epsilon =0$), Equations (4)–(6) revert to the standard SIR model. Note that the source term on the RHS of (4) and (6) is a specific convolution integral, which will be discussed further below. Parameter $\epsilon$ is the dimensionless kinetic rate of reaction (3), describing the loss of immunity of infected individuals. It is normalized with the rate of recovery as defined in ref. [7]: namely, if $I$ is the average characteristic time for the loss of immunity, then $\epsilon =\frac{L}{I}$, where $L$ was defined above as the average characteristic time for the recovery of infected individuals (e.g., assumed 14 days, for the case of COVID-19). In general, we expect small values of $\epsilon$, the case of large $\epsilon$ corresponding to the rather unrealistic case of the instant loss of immunity for all recovered individuals (because of its significant mathematical interest, however, this case will be also explored as well). Parameter $\theta \ge 0$ is a dimensionless time delay, normalized with $L$, after which reinfection commences. In our notation, it measures the number of time intervals, of a dimensional duration of about half a month, for COVID-19, which must elapse before reinfection is allowed. In the simulations below, we will typically focus on the range $0<\theta <24$ (the latter corresponding to about one-year delay). Finally, it is important also to note that Equation (9) does not include any of the two reinfection parameters.

The solution of (4)–(5) in the numerical simulations that will follow is subject to the initial conditions:

$$s\left(0\right)=0.99;i\left(0\right)=0.01; r\left(0\right)=0$$

(7)

while variables $s$, $r$, and $i$ obey the closure relation:

$$s+i+r=1$$

(8)

We must note that in the following results, both numerical and asymptotic, the solutions of (4)–(5) satisfy the physical constraints, $0\le s\le 1$, $0\le r\le 1$ and $0\le i\le 1$. Proving the validity of these conditions, in general, will require elaborate mathematical proofs which will not be attempted in this paper.

We now focus on deriving the gain and loss terms in (4)–(5). Consider the set $\mathcal{R}$ of recovered individuals. It consists of two different subsets, ${\mathcal{R}}_1$ and ${\mathcal{R}}_2$. Subset ${\mathcal{R}}_1$ contains infected individuals recovered at time $\tau$, such that $t-\theta \le \tau$, which are not yet subject to losing immunity (see also schematic in Appendix C). Conversely, the members of set ${\mathcal{R}}_2$ are individuals recovered at time $0<\tau \le t-\theta$. Thus, they have started and continue losing immunity. Subset ${\mathcal{R}}_2$ is empty if $t<\theta$, hence the need of the Heaviside step function $H\left(t-\theta \right)$, defined as the following:

$$H\left(x\right)=\{\begin{array}{c}1, x>0\\ 0, x<0\end{array}$$

Next, we consider a differential tranche of members in ${\mathcal{R}}_2$, created during time $\Delta \tau$. Its size is equal to $i\left(\tau \right)\Delta \tau$. Assuming first-order reaction kinetics for (3), this tranche loses members at the rate $\epsilon i\left(\tau \right)\Delta \tau$. At time $t$, therefore, its magnitude has been reduced to $\epsilon i\left(\tau \right)\Delta \tau \exp \left(-\epsilon \left(t-\theta -\tau \right)\right)$. Hence, the total loss at time $t$ is the integral over all $\tau$, namely $\epsilon {{\displaystyle \int }}_0^{t-\theta }i\left(\tau \right)\exp \left(-\epsilon \left(t-\theta -\tau \right)\right)d\tau H\left(t-\theta \right)$. We conclude that this is indeed the term on the RHS of (4) and (5).

Equations (4)–(6) will be solved subject to the typical initial conditions $i\left(0\right)=i_0, s\left(0\right)\equiv s_0=1-i_0, r\left(0\right)=0$, where $i_0\ll 1$ (see also (7)). Before we proceed further, we must make a couple of additional observations: First, Equations (4)–(6) are actually a system of only two differential equations, e.g., (4)–(5), expressed in terms of two variables, $s$ and $i$, since variable $r$ follows directly from the closure relation (8). In fact, one can simplify even further by solving Equation (5) for $s$, and substituting in (6). Then, one obtains the single non-linear, integrodifferential equation:

$$\frac{b}{\left(b+1\right)}{\left(\frac{i^{\prime }}{i}\right)}^{\prime }+i^{\prime }+i=\epsilon {{\displaystyle \int }}_0^{t-\theta }i\left(\tau \right)\exp \left(-\epsilon \left(t-\theta -\tau \right)\right)d\tau H\left(t-\theta \right)$$

(9)

Here, we defined $b\equiv \frac{1}{R_0-1}$, and assumed $R_0>1$ (which is required for the spreading of infection). Because $b$ is the inverse of $R_0-1$, it is expected to vary between a value of about 0.25 (for an aggressive contagion with $R_0=5$) to values much larger than 1 (when contagion is weak, $R_0\approx 1$).

Equation (5) generalizes the SIR model to the case of reinfections. Its solution and its dependence on the various parameters will be the main focus of this paper. We must note that when $\epsilon =0$, which is the SIR case, Equation (9) reverts to the standard SIR equation. This can be solved exactly with the solution [7]:

$$i\left(t\right)=1-r-s_0\exp \left(-R_{0}r\right); \mathrm{where} t={{\displaystyle \int }}_0^r\frac{du}{1-u-s_0\exp \left(-R_{0}u\right)}$$

(10)

For future reference we will also write (10) as

$$s\approx \exp \left(-R_{0}r\right)$$

(11)

where we assumed that the onset of the infection is at $s_0\approx 1$ (see (7)). From Equation (9) we find that in the limit corresponding to $i\approx 0$, the final value $r_{\infty }$, which solves $1-r_{\infty }-s_0\exp \left(-R_0{r}_{\infty }\right)=0$ is approached exponentially fast in time, with an exponent that is equal to $1-R_0{s}_{\infty }>0$.

2.1. Some Special Results

It is now worthwhile to first present some special results, the derivation of which is detailed in Appendix A.

i.

A general expression for r(t)

Because of its linear nature, Equation (6) can be used to also provide an expression of r(t) in terms of $i\left(t\right)$. The result is (Appendix A):

$$r\left(t\right)={{\displaystyle \int }}_0^{t}i\left(\tau \right)\left\{1-H\left(t-\tau -\theta \right)\left[1-\exp (-\epsilon \left(t-\tau -\theta \right)\right]\right\}dt$$

(12)

As expected, $0\le r\left(t\right)\le 1$.

ii.

The θ = 0 limit

When $\theta =0$, one anticipates that the loss term in (6) will reduce to $\epsilon r\left(t\right)$, since there is no delay in time. This is indeed the case, as shown in Appendix A. Now, the relevant equations become the following:

$$s^{\prime }\left(t\right)=-R_{0}si+\epsilon r\left(t\right)$$

(13)

$$i^{\prime }\left(t\right)=R_{0}si-i$$

(14)

$$r^{\prime }\left(t\right)=i-\epsilon r\left(t\right)$$

(15)

This important limit will be explored in considerable detail. In this case, one can also use the following alternative equation to (9):

$$b{v}^{″}-{v^{\prime }}^2-\left(1+2\epsilon \right)v{v}^{\prime }+\epsilon \left(\frac{1}{1+\epsilon }+b\right)v^{\prime }+\epsilon v\left(1-v\left(1+\epsilon \right)\right)=0$$

(16)

where we defined the following:

$$v=\frac{1}{\left(1+\epsilon \right)}-\left(b+1\right)r$$

(17)

This change of variables was motivated by the fact that $v(\infty )=0$, thus allowing for an easier asymptotic analysis.

iii.

The large ε limit

Even though unrealistic, it is worth examining the limit $\epsilon \gg 1$, corresponding to fast kinetics. In this limit, and after a delay time of $\theta$ is reached, all recovered individuals are instantly susceptible to infection. Now, the gain-loss terms in Equations (5), (6) and (9) reduce to (see Appendix A):

$$\underset{\epsilon \to \infty }{\lim }\epsilon {{\displaystyle \int }}_0^{t-\theta }i\left(\tau \right)\exp \left(-\epsilon \left(t-\theta -\tau \right)\right)d\tau H\left(t-\theta \right)=i\left(t-\theta \right)H\left(t-\theta \right)$$

(18)

and the corresponding equations read as follows:

$$s^{\prime }\left(t\right) = -R_{0}si+i\left(t-\theta \right)H\left(t-\theta \right)$$

(19)

$$i^{\prime }\left(t\right)=R_{0}si-i$$

(20)

$$r^{\prime }\left(t\right)=i-i\left(t-\theta \right)H\left(t-\theta \right)$$

(21)

and

$$\frac{b}{\left(b+1\right)}{\left(\frac{i^{\prime }}{i}\right)}^{\prime }+i^{\prime }+i=i\left(t-\theta \right)H\left(t-\theta \right)$$

(22)

Here, we find that the integral term reduces to a difference term, and the integrodifferential equations reduce to differential-difference equations. We noted in the introduction that a number of models in biological population dynamics have used a variety of models, some based on integrodifferential equations and some on differential-difference equations. Here, we show that the latter is the limit of the integral term in the case of fast kinetics, thus unifying the two different approaches.

iv.

Equilibrium states

When reinfections are present, the above formulation admits non-trivial equilibrium state solutions (denoted by subscript $\infty$). Assuming that these states are asymptotically stable (conditions for which will be developed below), their values are (Appendix A):

$$s_{\infty }=\frac{1}{R_0}=\frac{b}{\left(b+1\right)}$$

(23)

$$i_{\infty }=\frac{\epsilon }{\left(b+1\right)\left(1+\epsilon \left(1+\theta \right)\right)}=\frac{\left(R_0-1\right)\epsilon }{R_0\left(1+\epsilon \left(1+\theta \right)\right)}$$

(24)

$$r_{\infty }=\frac{\left(1+\epsilon \theta \right)}{\left(b+1\right)\left(1+\epsilon \left(1+\theta \right)\right)}=\frac{\left(R_0-1\right)\left(1+\epsilon \theta \right)}{R_0\left(1+\epsilon \left(1+\theta \right)\right)}$$

(25)

One concludes that allowing for reinfections ($\epsilon \ne 0$) leads at large times (assuming that such states are stable) to a constant, non-zero fraction in the population of infected individuals, $i_{\infty }$, which increases as the intensity of infection ($R_0$) increases. Moreover, note the important result that even for infinitesimally small values of $\epsilon$, the asymptotic values of $s_{\infty }$ and $r_{\infty }$ are not those predicted by the traditional SIR model (Ref. [7], and Equation (10)). The corollary is that the traditional SIR model is unstable, becoming singular at large times, because even a small value of $\epsilon$ will produce discontinuities in the final values $s_{\infty }$ and $r_{\infty }$. Since these values define the state of herd immunity, it is clear that the latter cannot be predicted using the SIR model and (9). We will comment on these features in more detail in the subsequent sections.

2.2. Numerical Results

Next, we consider the numerical solution of (4)–(6) and examine the sensitivity of the results to the values of the three parameters $\epsilon ,\theta , b$. These define a three-dimensional space, one boundary of which, namely the plane $\epsilon =0$, $\theta \gg 1$, corresponds to the conventional SIR model, which is well studied. For additional useful insights, we will focus on the other two boundaries of the $\left\{\epsilon , \theta \right\}$ parameter space, namely, the planes $\theta =0$, and $\epsilon \gg 1$. These correspond to Equations (16) and (22), respectively. We expect that they will provide prototypical regimes, valid for all other parameter values as well.

1.

The case θ = 0

In the absence of a delay in time, $\theta =0$, the problem reduces to Equations (13)–(15), which does not contain an integral source term. Numerical simulations for four different values of $\epsilon$ ($\epsilon =$ 10, 1, 0.1, and 0.01), spanning a range of behaviors, are shown in Figure 1. We observe the following: At large $\epsilon$ (Figure 1a) the variation of all three populations is monotonic, and after about an $O\left(1\right)$ time interval, they asymptotically approach a stable state (which will be shown to be a stable node), which from (23) reads as follows: $s_{\infty }=\frac{1}{R_0}=\frac{b}{\left(b+1\right)}$, $i_{\infty }=\frac{R_0-1}{R_0}=\frac{1}{\left(b+1\right)}$, $r_{\infty }\approx \frac{1}{\left(b+1\right)\epsilon }\ll 1$. Asymptotically, the fraction of the recovered population is small and of order ${\epsilon }^{-1}$, and that there is a constant, non-zero fraction of infected individuals, that increases with decreasing $b$ (increasing $R_0$), while the fraction of susceptible individuals $s_{\infty }$ decreases with increasing $R_0$. In this limit, closed-form results are also possible. At the large $\epsilon$ limit one obtains the analytical result

$$i\left(t\right)\approx \frac{1}{\left(b+1\right)+\left\{\frac{1}{i_0}-\left(b+1\right)\right\}\exp \left(-\frac{t}{b}\right)} ; r\approx \frac{i}{\epsilon } ; 1-s\approx \epsilon r$$

(26)

suggesting an asymptotic decay at the rate $b^{-1}$.

As $\epsilon$ decreases, the behavior changes, initially to a slight, and progressively to a stronger, undulation (Figure 1b–d). A number of distinct waves for the susceptible fraction appear, the amplitude of which decreases with time. We will show that such oscillatory behavior corresponds to a spiral (a “damped oscillator” [20]). It is worth examining Figure 1c,d in more detail. Here, over a length of time of about 2, in our dimensionless notation, the variables closely follow the SIR model of Ref. [7] (where $\epsilon =0$) which is described by the SIR limit (16), namely:

$$b{v}_0^{″}-v_0^{\prime }{}^2-v_0{v}_0^{\prime }=0$$

(27)

the integration of which gives the explicit SIR result in (10). While approximating well the early part of the curves (for example see Figure 1c), this result fails at larger times, since its asymptotic limit, reached exponentially fast, is not the correct limit $v_{\infty }=0$, but rather the solution $v_{\infty ,0}>0$ of the algebraic equation:

$$v_{\infty ,0}=1+b\ln \frac{z_{\infty ,0}+b}{\left(1+b\right)s_0}$$

(28)

As a result, and at about $t=2$ in Figure 1c, the number of susceptible individuals starts increasing, due to the loss of immunity of some of the recovered population, thus leading to an increase in infections, hence to the appearance of a secondary infection wave. This wave eventually wanes, but because of additional reinfections, infection continues, with the emergence of secondary waves, now of a lower amplitude, and so on. A sequence of such waves, with constant frequency, but with an amplitude that appears to decrease exponentially in time, follows. Figure 1d is a clear demonstration of this behavior.

To better analyze these results, consider the solution of (16) at large times. By neglecting higher-order terms and linearizing (namely, by conducting a linear stability analysis), we have

$$v^{″}+\left(m_0+\epsilon \right)v^{\prime }+m_0\left(1+\epsilon \right)v\approx 0$$

(29)

the solution of which is

$$v=\exp (-{\lambda }_{0}t)f\left(t\right)$$

(30)

Here, $f\left(t\right)$ satisfies

$$f^{″}+\frac{\Delta }{4}f=0$$

(31)

where we have defined

$${\lambda }_0\equiv \frac{\left(\epsilon +m_0\right)}{2}>0, m_0\equiv \frac{\epsilon }{b\left(1+\epsilon \right)}>0, \Delta \left(\epsilon ,b\right)\equiv 4m_0-{\left(\epsilon -m_0\right)}^2$$

(32)

When $\Delta <0$, namely for sufficiently large $\epsilon$, such that ${\left(\epsilon -m_0\right)}^2>4m_0$, the solution of (31) is a monotonic approach to a stable solution (a stable node) (Figure 1a). In the opposite case, ${\left(\epsilon -m_0\right)}^2<4m_0$, it is a sinusoidal wave (a spiral) (Figure 1b–d), with frequency ${\omega }_0=\frac{1}{2}\sqrt{4m_0-{\left(\epsilon -m_0\right)}^2}$. Here, within two arbitrary constants $C$ and $\phi$, we have

$$f=C\cos \left({\omega }_{0}t+\phi \right)$$

(33)

which, when inserted in (29), provides the damped oscillator solution. The frequency ${\omega }_0$ increases with decreasing $b$, namely with larger $R_0$. It also increases with $\epsilon$, until it reaches a maximum, above which it decreases with increasing $\epsilon$, and vanishes when $\Delta =0$, beyond which the asymptotic approach is no longer oscillatory.

Figure 2 shows a plot of variable $f$ as a function of the rescaled time ${\omega }_{0}t$ for the latter case. The approach to the asymptotic solution is relatively fast, roughly corresponding to $t>4{\lambda }_0{}^{-1}$. For example, for $\epsilon =$ 1, 0.1 and 0.01, and $b=0.25$, it is of the order of 2.5, 17, and 160, respectively (and in dimensional notation for COVID-19 equal to 1.25 months, 8.5 months, and 6 years, respectively). The latter value is of course unrealistic, given the various other circumstances that will emerge during that period.

If we express the solution of (29) as $v~\exp \left(at\right)$, we find the characteristic equation

$$a^2+\left(m_0+\epsilon \right)a+m_0\left(1+\epsilon \right)=0$$

(34)

for the generally complex number $a$. By further denoting $a=x+iy$, we find the following: (i) If $\Delta <0$, then $x=-{\lambda }_0\pm \sqrt{-\Delta }<0$ and $y=0$, while (ii) if $\Delta >0$, then $x=-{\lambda }_0<0$, and $y=\sqrt{\Delta }$. Figure 1a corresponds to the first case ($\Delta =$ −25.95), while Figure 1b–d, correspond to the second case ($\Delta =$ 7, 1.93, 0.16).

The demarcation of the region where the approach to the asymptotic state changes from that of an exponential decay to a “damped oscillator” is shown in Figure 3. The delineating curve is $\Delta =0$, namely $\epsilon =m_0-2\sqrt{m_0}$. (There is also another branch where $\Delta =0$, but which is not shown in Figure 3 because it corresponds to the physically unrealistic results of very large $R_0$).

We conclude that in the absence of a delay in time, two stable asymptotic regimes emerge when reinfections are allowed, both leading to an asymptotically constant state. The traditional SIR model cannot capture the long-time behavior, even when $\epsilon$ is infinitesimally small. Indeed, while the SIR model is followed for small times, (e.g., less than about 2 in Figure 1c), a different regime sets in when time is larger, leading to a different asymptotic state. The conclusion is that the SIR model is singular at times of order ${\epsilon }^{-1}$.

For future use it is useful to also plot the results in the phase plane $\left(s, r\right)$. Figure 4 shows results for the case $b=1$ ($R_0=2)$. Plotted in the figure are also two limiting curves, one corresponding to the SIR model (Equation (11)) and another corresponding to the case where $s\left(t\right)$ reaches an extremum (maximum or minimum). This can be readily shown to be the curve $r=R_0\frac{s\left(1-s\right)}{\epsilon +s{R}_0}$. When the solution trajectory intersects this curve, its derivative vanishes and its direction changes. The transition from a stable node to a spiral, namely to a wave behavior, as the value of $\epsilon$ decreases is clear in Figure 4.

2.

The limit of large ε

Consider, next, a different boundary in the parameter space, the plane defined by $\epsilon \gg 1$. This is the case where reinfection occurs instantly, after a delay time $\theta$ has elapsed. Figure 5 shows numerical results corresponding to four different values of $\theta$, selected such that $\theta <{\theta }_m$, ${\theta }_m<\theta <{\theta }_c$, and $\theta >{\theta }_c$, respectively, where the two critical values ${\theta }_m$ and ${\theta }_c$ depend on $b$ (and for Figure 5, they are equal to ${\theta }_m\left(1\right)=0.96$ and ${\theta }_c\left(1\right)=7.8$, respectively). We observe the following:

In Figure 5a,b, where, $\theta <{\theta }_m$ and $\theta <{\theta }_c$, respectively, the solution approaches at the large $\epsilon$ limit the equilibrium state (23)–(25)

$$s_{\infty }=\frac{1}{R_0}=\frac{b}{\left(b+1\right)}$$

(35)

$$i_{\infty }=\frac{1}{\left(b+1\right)\left(1+\theta \right)}=\frac{\left(R_0-1\right)}{R_0\left(1+\theta \right)}$$

(36)

$$r_{\infty }=\frac{\theta }{\left(b+1\right)\left(1+\theta \right)}=\frac{\left(R_0-1\right)\theta }{R_0\left(1+\theta \right)}$$

(37)

The approach is different depending on the value of $\theta$. When $\theta <{\theta }_m$ (Figure 5a), there is a monotonic exponential decay (corresponding to a stable node); when ${\theta }_m<\theta <{\theta }_c$ (Figure 5b), the solutions have the features of a damped oscillator (spiral); and when $\theta >{\theta }_c$ (Figure 5c,d), the asymptotic equilibrium is no longer stable, and the variables are approaching a periodic structure pattern of constant amplitude, with a period that depends on $\theta$ and $b$. This transition is consistent with a Hopf bifurcation [21], expected to arise when the rate of growth of the amplitude at large times becomes positive, and suggests the possibility of a limit cycle behavior. We note that in the case of large $\theta$, there are periods during which the infected fraction is infinitesimally small. We believe that this characteristic of any reinfection process, since after a time $\theta$ has elapsed, recovered individuals can become susceptible leading the system away from conditions of herd immunity, hence the onset of new infection waves. Of course, all this is under the assumption that all physiological and biological processes remain the same (no vaccinations, no changes in the biology of the recovered individuals, etc.).

To quantify the results shown in Figure 5, we consider the linear stability of Equation (22) around the equilibrium state (36). Assuming a small perturbation $\varpi$ of $i$ around $i_{\infty }$ and linearizing at large times we find

$$\frac{1}{m_{\infty }}{\varpi }^{″}+{\varpi }^{\prime }+\varpi =\varpi \left(t-\theta \right)$$

(38)

where we defined

$$m\equiv \frac{\epsilon }{b\left(1+\epsilon \left(1+\theta \right)\right)} \mathrm{and} m_{\infty }\equiv \underset{\epsilon \to \infty }{\lim }m\left(\theta ,\epsilon \right)=\frac{1}{b\left(1+\theta \right)}$$

(39)

The solution of (38) (which is the counterpart of (29) in the $\epsilon \gg 1$ limit) is the exponential

$$\varpi ~\exp \left(at\right)$$

(40)

where $a$ is a complex number. Substituting in (38) we find

$$\frac{1}{m_{\infty }}{a}^2+a+1-\exp \left(-a\theta \right)=0$$

(41)

and by further taking $a=x+iy$ we obtain the system of the two algebraic equations

$$\frac{1}{m_{\infty }}(x^2-y^2)+1+x=\exp \left(-x\theta \right)\cos \left(y\theta \right)$$

(42)

$$\frac{2}{m_{\infty }}xy+y=-\exp \left(-x\theta \right)\sin \left(y\theta \right)$$

(43)

where without loss, we can take $y>0$. Equations (42) and (43) do not have a simple explicit solution. However, one can show (Appendix B) that there exist two critical values ${\theta }_m\left(b\right)$ and ${\theta }_c\left(b\right)$ delineating three different regimes: (i) a stable node, if $0<\theta <{\theta }_m\left(b\right)$; (ii) a spiral, if ${\theta }_m\left(b\right)<\theta <{\theta }_c\left(b\right)$; and, (iii) a periodic state, arising from a Hopf bifurcation, if $\theta >{\theta }_c\left(b\right)$. As expected, the limit $\theta =0$ in (42)–(43) coincides with the corresponding limit of (30)–(31) in the case of large $\epsilon$. For future use, it is also worth presenting the results in the phase plane $\left(s, r\right)$ (Figure 6), which shows clearly the three asymptotic states. This figure will be used in the discussion section that follows.

For completeness we note the analytical result in the case of small $\theta$ (see Appendix B)

$$i\left(t\right)=\frac{1}{\left(b+1\right)\left(1+\theta \right)+\left\{\frac{1}{i_0}-\left(b+1\right)\left(1+\theta \right)\right\}\exp \left(-\frac{t}{b}\right)} ; s\left(t\right)\approx 1-\left(1+\theta \right)i\left(t\right) ; s\left(t\right)\approx 1-\frac{\left(1+\theta \right)}{\theta }r\left(t\right)$$

(44)

suggesting an asymptotic decay at the rate $b^{-1}$. Equation (44) reproduces the numerical result of Figure 5a, while it is also consistent with the $\theta =0$ case discussed previously at the large $\epsilon$ limit.

Figure 7 shows the dependence of ${\theta }_m\left(b\right)$ and ${\theta }_c\left(b\right)$ on the parameter $b$ (hence on the reproduction number $R_0$), obtained from the solution of (42)–(43) (Appendix B). As the intensity of the epidemic increases (larger $R_0$), it is easier for the asymptotic state to become a periodic pattern of constant amplitude and frequency. We anticipate somewhat similar behavior for finite or smaller $\epsilon$, which is the next case to be discussed.

3.

The general case

Having explored the three limiting boundaries in the ($\epsilon , \theta$) parameter space, we can now consider the more general problem, where $\epsilon$ is finite and $\theta$ is non-zero. We expect that, in general, there will emerge three different regimes: Two stable asymptotic states, for sufficiently small values of $\theta$, the nature of which could be either a stable node or a spiral, and an unstable equilibrium state at sufficiently large $\theta >{\theta }_c\left(b,\epsilon \right)$, which will lead to an oscillatory periodic state (a limit cycle). Numerical results are shown in Figure 8 and Figure 9 for a fixed $b$, for different values of $\theta$, and for two different values of $\epsilon$, respectively. In Figure 8 $\epsilon$ is chosen so that it falls in the upper region of Figure 3 (where $\Delta <0$) where the asymptotic state for $\theta =0$ is a stable node. Figure 9, on the other hand, corresponds to a smaller value of $\epsilon$, chosen to be in the lower region in Figure 2, (where $\Delta >0$) where the asymptotic state for $\theta =0$ is a spiral.

Figure 8a,b ($\theta =0.5$ and 5) suggest that the solution approaches a stable equilibrium, monotonically in Figure 8a, and as a damped oscillator in Figure 8b. The other figure, Figure 8c,d (where $\theta =$ 10 and 20) indicate that the equilibrium state is unstable and becomes a periodic solution. This is consistent with the limit of large $\epsilon$, discussed before, and where the order of appearance of the asymptotic regimes as $\theta$ increases from zero were monotonic, a damped oscillator or a periodic state, respectively. Two critical values ${\theta }_m\left(b,\epsilon \right)$ and ${\theta }_c\left(b,\epsilon \right)$ separate these asymptotic regimes.

On the other hand, Figure 9 shows that there are only two regimes, a damped oscillator, for sufficiently small $\theta$ (Figure 9a) and a periodic solution, at larger $\theta$ (Figure 9b). Again, a critical value ${\theta }_c\left(b,\epsilon \right)$ separates these two regimes. Here, the sequence of the asymptotic regimes as $\theta$ increases from zero are a damped oscillator and a periodic state, respectively.

To confirm these results, we conducted a linear stability analysis. Appendix C shows that now the rate of growth $a$ satisfies the algebraic equation:

$$\frac{1}{m}{a}^2+a+1-\frac{\epsilon }{a+\epsilon }\exp \left(-a\theta \right)=0$$

(45)

As expected, this reduces to the corresponding limits, namely Equation (32) for the case of $\theta =0$, and Equation (41) for the case of large $\epsilon$, respectively, With the usual representation $a=x+iy$, we further find (Appendix C).

$$\left(x^2-y^2+mx+m\right)\left(x+\epsilon \right)-\left(2x+m\right)y^2-\epsilon m\exp \left(-x\theta \right)\cos \left(y\theta \right)=0$$

(46)

$$\left[\left(x^2-y^2+mx+m\right)y+\left(2x+m\right)y\left(x+\epsilon \right)\right]+\epsilon m\exp \left(-x\theta \right)\sin \left(y\theta \right)=0$$

(47)

The solution of (46) and (47) is similar to the previous case of large $\epsilon$ (Figure 7) only if $\epsilon$ is larger than a ctitical value, $\epsilon >{\epsilon }_c\left(b\right)$, at which point ${\theta }_m=0$. For such values of $\epsilon$, there are two critical values, ${\theta }_m\left(\epsilon ,b\right)$ and ${\theta }_c\left(\epsilon ,b\right)$, that demarcate the three regions where there is monotonic exponential decay, a damped oscillator or a constant oscillation (Figure 10). However, when $\epsilon$ is smaller than the critical value, $\epsilon <{\epsilon }_c\left(b\right)$, (such that $\Delta >0$), a critical value ${\theta }_m>0$ does not exist, and there is only one critical value ${\theta }_c\left(\epsilon ,b\right)$ that separates a regime of a damped oscillator from that of a constant oscillation. The value ${\epsilon }_c\left(b\right)$ is the same as the one plotted in Figure 3, since it demarcates the region at which $\Delta =0$ and $\theta =0$. As shown in Figure 3, ${\epsilon }_c\left(b\right)$ decreases as $b$ increases or $R_0$ decreases. Finally, it is also not difficult to show that the asymptotic behavior of ${\theta }_c$ when $\epsilon \ll 1$ is given by the expression

$${\theta }_c\approx \frac{\sqrt{2\left(b+1\right)}}{\sqrt{\epsilon }}$$

(48)

which is consistent with Figure 10b. For completeness, the results of Figure 8 are also plotted in the phase space plane $\left(s, r\right)$ (Figure 11). The three asymptotic states are clearly portrayed. Similar phase portraits can also be obtained for the simulations shown in Figure 9.

3. Discussion

To provide a physical interpretation of the results obtained, it is useful to analyze the behavior in the phase space $\left(s, r\right)$, as portrayed in Figure 4, Figure 6, and Figure 11. In all cases, the solution trajectories lie within the right triangle bounded by the two coordinates $s, r$ (which are varying in (0,1)). As expected from the closure relation (8), the dotted green line (the hypotenuse) is the case of zero infections, and any line parallel to it denotes a fixed value of $i$. In the Figures, the color red denotes the trajectory corresponding to the SIR model (where $\epsilon =\theta =0$), given by Equation (11). For a better understanding of the results, it is worth focusing first on Figure 4, where $\theta =0$, and subsequently on Figure 6, where $\epsilon \gg 1$.

Consider, first, Figure 4. In all these cases the equilibrium state cannot be one in the absence of infections, (namely $i_{\infty }=0$). Otherwise, any recovered individual will become instantly susceptible, thus leading to an increase in $s$, at which point infection will recommence. Instead, the equilibrium state must correspond to $s_{\infty }=\frac{1}{R_0}$ (where the time derivative of $i$ also vanishes). In the first two panels (Figure 4a,b), where the asymptotic solution is a stable node, the solution trajectory is different from that of the SIR model and approaches the different equilibrium state $s_{\infty }=\frac{1}{R_0}=0.5$. In the limit of large $\epsilon$, the solution is a straight line, following Equation (26), as approximately shown in Figure 4a. Because of the absence of time delay and because $\epsilon$ is sufficiently large, reinfection is instant and fast, thus starting to interfere with the SIR process immediately. Thus, the infection process is reinforced from its onset, as new supplies of susceptible individuals are provided fast from those just recovered. Hence, one expects a fast and exponentially asymptotic approach to the equilibrium point.

In Figure 4a,b the infected fraction increases monotonically (as also shown in Equation (26)), while the susceptible fraction decreases monotonically as it approaches $s_{\infty }$. At equilibrium, there is a constant and fast conversion of recovered populations to infected ones, and vice versa, with infection and reinfection rates balancing each other, such that the system is at an equilibrium. Any deviation from equilibrium is rapidly restored. The characteristic time for this approach decreases as $R_0-1$ increases. This dependence on $R_0$ is notable, as the latter incorporates social behavior aspects, such as spatial distancing, face coverings, etc., increasing as the latter decrease.

When the asymptotic attractor is a spiral (Figure 4c,d), where $\epsilon$ is relatively small, infection first follows the SIR behavior described in (11). Because of reinfection, however, the rate of decrease of the susceptible population slows down, as some of the recovered individuals lose immunity to join the susceptible population. The maximum in the fraction of infected individuals occurs when $s$ approaches the value $s_{\infty }=\frac{1}{R_0}$ (see Equation (5)), which in the phase plane is a horizontal line parallel to the $r$-axis. This endows the solution trajectory with a counter-clockwise motion, since the second derivative of $i$ at that point is negative, as can be readily shown by taking the derivative of (5) and evaluating at $s_{\infty }=\frac{1}{R_0}$. As time proceeds, and because of reinfection, the trajectory will intersect at some time the curve $r=R_0\frac{s\left(1-s\right)}{\epsilon +s{R}_0}$, shown in light blue in Figure 4. This will signal that $s$ has reached a minimum (the RHS of (13) vanishes at that point). The counter-clockwise motion continues with decreasing infection populations, until $s$ increases to the value $s_{\infty }=\frac{1}{R_0}$, where the fraction of infected individuals now reaches a minimum, following which infection rates increase again. At a sufficiently large value of $s$ the rates of infection exceed the rates of replenishment, due to reinfections, the trajectory intersects again the curve $r=R_0\frac{s\left(1-s\right)}{\epsilon +s{R}_0}$, but at a different point, and a new infection wave sets in. The difference with the first wave is that because of the smaller overall susceptible fraction (which is close to $\frac{1}{R_0}$), the infection rate is correspondingly smaller, hence the trajectory is closer to the equilibrium state, which is given by $s_{\infty }=\frac{b}{\left(b+1\right)}$, $i_{\infty }=\frac{\epsilon }{\left(b+1\right)\left(1+\epsilon \right)}$, $r_{\infty }=\frac{1}{\left(b+1\right)\left(1+\epsilon \right)}$. Eventually, the system starts spiralling towards this state, which it approaches with a decreasing amplitude. At equilibrium, there will be progressively smaller infection waves, driven by the conversion of those recovered to susceptible, which are progressively attenuated due to the relatively low rates of infection and the fact that close to the equilibrium state $s=\frac{1}{R_0}$, the net rate of generation of infected individuals is zero. In Figure 6 and Figure 11, the stable node and spiral behaviors are very similar to Figure 4. In Figure 6a,b, the value of the delay time $\theta$ is relatively small and thus interferes relatively fast with the SIR model dynamics, thus the trajectory starts deviating from that of the SIR model relatively soon. A similar behavior occurs in Figure 11a,b.

When the delay time $\theta$ is sufficiently large, for example as in Figure 6c,d, the trajectory does not depart from the SIR model until it almost reaches its asymptotic limit. When $t$ begins exceeding $\theta$, Equation (21) becomes the following:

$$r^{\prime }\left(t\right)=i-i\left(t-\theta \right)$$

(49)

Following that time, the infection fraction $i$ will continue decreasing. The trajectory approaches the hypotenuse, while $i\left(t-\theta \right)$ will continue increasing, due to the time lag resulting from the delay time. At some time, the RHS of (49) will vanish, and the trajectory will commence reversing its direction ($r^{\prime }\left(t\right)<0$), consistent with the onset of reinfections. As long as $R_{0}s<1$, the fraction of infected individuals will keep decreasing (see Equation (5)), implying an even closer approach of the trajectory to the hypotenuse, until $s$ reaches the value $\frac{1}{R_0}$, at which point the trajectory will start curving away from it. Because of the corresponding increase in $i$, as time increases, there will be a time at which the RHS of (49) will become positive again, thus the trajectory will start curving towards increasing values of $r$, and the cycle will recommence. Thus, the periodic trajectory obtained at sufficiently large values of the delay time $\theta$ consists of the following: one segment, almost coinciding with the SIR curve, and another segment being very close to the hypotenuse, where the infection fraction is negligible (Figure 6d and Figure 10d). The characteristics of the latter segment can be readily obtained by using Equations (19) and (21) and by taking $i\ll 1$. Then,

$$s^{\prime }\left(t\right)\approx i\left(t-\theta \right)$$

(50)

$$r^{\prime }\left(t\right)\approx -i\left(t-\theta \right)$$

(51)

Hence,

$$s\left(t\right)\approx 1-r\left(t\right)\approx {{\displaystyle \int }}_{\theta }^{t}i\left(\tau -\theta \right)d\tau +s_{\infty ,SIR}={{\displaystyle \int }}_0^{t-\theta }i\left(\tau \right)d\tau +s_{\infty ,SIR}$$

(52)

where, the SIR limit, $s_{\infty ,SIR}$ satisfies Equation (10) in the limit $i\to 0$, namely:

$$s_{\infty ,SIR}=\exp \left(-R_0\left(1-s_{\infty ,SIR}\right)\right)$$

(53)

Using the expressions for the SIR limit in (52), we find

$$s\left(t\right)\approx r_{SIR}\left(t-\theta \right)+s_{\infty ,SIR}$$

(54)

where $r_{SIR}\left(t\right)$ is obtained by (10), namely:

$$t={{\displaystyle \int }}_0^{r_{SIR}}\frac{du}{1-u-s_0\exp \left(-R_{0}u\right)}$$

(55)

At large times we have,

$$s\left(t\right)\to r_{\infty ,SIR}+s_{\infty ,SIR}=1$$

(56)

The corresponding pattern has a period equal to $\theta$.

In the general cases depicted in Figure 8 and Figure 9, the behavior can be explained following a similar reasoning (see also Figure 11). The analysis is more complex at smaller values of $\theta$, and for more general values of $\epsilon$, since analytical results are not easily obtained. It is important to note, however, that when $\epsilon$ is sufficiently small, $\epsilon <{\epsilon }_c\left(b\right)$, which as Figure 3 shows is likely to be the case for most realistic situations, the predominant pattern is a spiral. For a periodic pattern to emerge requires large delay times, to exceed ${\theta }_c$, as also shown in Equation (48) and Figure 10. One concludes that in the general case, where $\epsilon$ is not large, the asymptotic behavior is likely to be a wave with a decaying amplitude and constant frequency, approaching asymptotically a constant state.

4. Conclusions

In this paper we considered the extension of the traditional SIR model, as further refined in ref. [7], to capture the possibility of reinfections. We proceeded by assuming that infected individuals are reinfected at a specified kinetic rate, and that they do not acquire permanent immunity. While this is a likely unrealistic assumption it does provide an interesting long-time behavior. By appropriately capturing the rate of reinfection in terms of a convolution integral, we examined the nature of the asymptotic state as a function of three dimensionless parameters—the normalized reinfection kinetic rate $\epsilon$, expected to be a small positive number; a delay time $\theta$ before which immunity is lost and reinfection commences; and the inverse of the distance of the reproduction number from unity, $b=\frac{1}{(R_0-1)}$, which appears to be the appropriate relevant parameter. We found that there are three possible asymptotic states. For sufficiently small $\theta$, the process approaches a stable attractor, either in the form of a stable node, at sufficiently large kinetic rates $\epsilon$, or in the form of a spiral, when $\epsilon$ is sufficiently small. This wave behavior is encountered at larger $b$ (namely at milder forms of epidemics). Such behavior will set in regardless of the value of the kinetic constant of reinfection, no matter how small that might be. A wave behavior will set in when $\epsilon$ is not large, assuming that $\theta$ is not too small. This is likely to be the practical case. Then, the asymptotic behavior will be in the form of waves of decreasing amplitude and in the approach to a steady state, in which the fraction of infected individuals is a non-zero constant. We explored its amplitude and frequency. When the delay time for the onset of the loss of immunity is sufficiently large, the asymptotic state is not stable, and the process will settle in a periodic pattern, which consists of a segment that parallels the standard SIR model, and of another segment, where the infection fraction is negligible. This pattern has a period proportional to $\theta$.

The model we used is based on a number of assumptions, some of which might be directly relaxed, and some that require significant additional work. One such assumption is that all recovered individuals are subject to reinfection. Clearly, this is not accurate, as fraction $r$ also includes individuals who perish. This can be rectified if we assume that a constant fraction $p$ of the population $r$ perishes. Then, the above remains valid subject to the substitution $\epsilon \to p\epsilon$. Under the assumptions considered, it captures a model of the transition of a pandemic to an endemic state. Importantly, we find that the herd immunity values predicted by the traditional SIR model are unlikely to materialize, since even in the presence of a small reinfection rate, the system will approach the value $1-s_{\infty }=1-\frac{1}{R_0}$, rather than the value in Equation (31) implied in the standard SIR model. The former value is also predicted when the large-time decay is algebraic (as in [8]).

We should also note that the present model assumes that conditions exist such that the simpler SIR model (in the absence of reinfections) can be applicable, which implies compartmental areas with negligible spatial gradients. If such conditions do not exist, for example, if mobility and diffusion are important parameters, then one could use the general model presented in reference [7] and expand models (3)–(6) to also include spatial gradients due to diffusion or other mobility processes. Of course, in any data comparisons, one would also need to account for vaccinations, etc., which were not included in our paper. This can be undertaken in further studies.