A refractory density approach to a multi-scale SEIRS epidemic model

Abstract

We propose a novel multi-scale modeling framework for infectious disease spreading, borrowing ideas and modeling tools from the so-called Refractory Density (RD) approach. We introduce a microscopic model that describes the probability of infection for a single individual and the evolution of the disease within their body. From the individual-level description, we then present the corresponding population-level model of epidemic spreading on the mesoscopic and macroscopic scale. We conclude with numerical illustrations taking into account either a white Gaussian noise or an escape noise to showcase the potential of our approach in producing both transient and asymptotic complex dynamics as well as finite-size fluctuations consistently across multiple scales. A comparison with the epidemiology of coronaviruses is also given to corroborate the qualitative relevance of our new approach.

[inst1]organization=Centre Inria d’Université Côte d’Azur,addressline=2004 Rte des Lucioles, city=Biot, postcode=06410, country=France

[inst2]organization=Universite Claude Bernard Lyon 1, CNRS, Ecole Centrale de Lyon, INSA Lyon, Université Jean Monnet, ICJ, UMR5208, Inria,city=Villeurbanne, postcode=69622, country=France

[inst3]organization=Bernstein Center for Computational Neuroscience,addressline=Philippstr. 13, Haus 6, city=Berlin, postcode=10115, country=Germany [inst4]organization=Technische Universität Berlin Institute of Mathematics,addressline=Straße des 17. Juni 136, city=Berlin, postcode=10623, country=Germany

[inst5]organization=Department of Mathematical Sciences “G. L. Lagrange”, Politecnico di Torino,addressline=Corso Duca degli Abruzzi 24, city=Torino, postcode=10129, country=Italy

1 Introduction

The susceptibility of an individual to an infectious disease as well as the infectiousness of an infected individual considerably depends on the time of the person’s last infection, or “infection age”. Such dependence on the infection course and history can be modeled by age-structured population models. These models have been used since the seminal paper by Kermack and McKendrick [17] for a wide variety of mathematical models in epidemics, see e.g. [2, 3, 6, 13, 18]. Taking into account the time since the last infection (TSLI) is extremely convenient to introduce a temporal heterogeneity in the population concerning an ongoing epidemic, and to explicitly model a complex disease evolution with consequences both within and between hosts.

Population models structured by this TSLI are deterministic macroscopic models that relate to the idealized limit of an infinitely large population. Two important questions arise here: first, a fundamental question is how such macroscopic models are linked to “microscopic models”, where the TSLI or other internal variables are tracked for each individual in the population. Microscopic models can be more directly related to the mechanisms underlying the spreading of disease and to parameters measured in clinical studies. However, these models are too complex to understand the collective dynamics of an epidemic. Second, populations in reality (say, a local community relevant to the disease or an age group) are not infinitely large and thus exhibit fluctuations that reflect the infections of single individuals. Therefore, this raises the question of whether there exists an intermediate “mesoscopic model” that combines the simplicity of age-structured population models with the ability to capture fluctuations due to finite population size. Again, the parameters of such a mesoscopic model should be linked to the microscopic parameters. Answering these questions requires a multi-scale modeling framework for epidemic diseases.

We found our inspiration in the neuroscience field, where the age-structured population dynamics models have already been used in the form of the refractory density method [7, 8, 11, 12, 14, 24, 27, 29, 28]. Indeed, in that case, the basic units at the microscopic scale are the neurons, whose probability to fire a spike is largely determined by the time since the last spike as well as the spike input from other neurons in a large network of neurons. Taking into account the time since the last spike is important because neurons exhibit a period of absolute and relative refractoriness after a spike, where the firing probability is strongly reduced. At the macroscopic scale, the evolution of a neuronal population is governed by a system of first-order partial differential equations called refractory-density equations (RDE). The RDE tracks the density of the times since the last spikes within a continuous phase state, thereby avoiding discrete states such as spiking, refractoriness, and resting. Importantly, the RDE is a bottom-up mean-field model that is directly derived from the single neuron equations at the microscopic scale. Simulations using this technique yield precise solutions for both equilibrium states and transient non-equilibrium regimes.

In the last two decades, there have been several advancements in the refractory-density method that offered the possibility for significant progress in the multi-scale modeling of cortical networks in the brain. First, the method has been adapted for the important case where spiking events are triggered by the threshold crossing of a neuronal membrane potential that is driven by biologically plausible Gaussian white [7] or colored noise [8, 26], while the classical RDE is based on a phenomenological hazard-rate-based spike-generation mechanism (“escape noise”) [14]. Second, the RDE has been generalized to finite-size populations [25, 28], resulting in a stochastic RDE that accurately reproduces the finite-size fluctuations at the mesoscopic scale. This advancement offers a unique and consistent description of a neural circuit at three levels of granularity: micro-, meso- and macroscopic scales. Additionally, there has been a further advancement in the generalization of RDE to accommodate complex, multi-dimensional [7], and even two-compartmental [9] neurons, where the neuronal state is given by multiple variables rather than a single membrane potential. These advancements collectively allow the application of RDE for constructing meso- and macroscopic models based on a wide range of microscopic models.

The main concept of this paper is to adopt the recent generalizations of the refractory-density method from neuroscience to epidemiology. In particular, we suggest a novel multi scale modeling framework that connects three different scales (micro, meso, macro). Moreover, the theoretical framework applies to two different noise mechanisms at the microscopic scale (Gaussian white noise and escape noise). We show their distinct effects on the form of the age-structure epidemic population model at the meso- and macroscopic scales. To this end, we consider an infectious disease with a relatively short infection period, followed by temporary immunity and the possibility of multiple subsequent infections (e.g., flu or SARS-CoV-2). Our model can be naturally structured within a classical SEIRS (Susceptible – Exposed – Infected – Recovered – Susceptible) epidemic framework, where we continuously track the interplay between viral load and the immune system within each individual.

Unlike the classical ordinary differential equation (ODE)-based compartmental models of epidemics (see, e.g., [19]), we choose a continuous phase space, focusing on the time since infection and the corresponding evolution of the interaction between the virus and the immune system. In this space, the SEIRS stages can be defined either strictly or loosely; see section 2 for a detailed explanation of our definition. In reality, an individual does not experience strict boundaries between SEIRS stages, and our approach allows us to avoid such additional assumptions. Instead, we propose that two characteristics—- the time since the previous infection and the interplay between viral load and immune response -— sufficiently describe an individual’s state within a population. We believe that avoiding strictly defined compartments and introducing a continuous space of states is a key advantage of our model.

Another advantage of our approach is that the spread of infectious disease at the mesoscopic and macroscopic population level is derived from the equations that stand for the dynamics of infection, disease development, and recovery of a single individual. A similar approach governs the derivation of macroscopic equations starting from their microscopic counterparts in kinetic theory, see e.g. [10, 20]. Our construction differs however significantly, as we explain in sections 2 and 3. In particular, the present paper is, to the best of our knowledge, novel in that it combines such a multi-scale (microscopic, mesoscopic, macroscopic) model with an age (in the sense of “infection age”) structured component, exploiting tools and analytical techniques coming from the so-called Refractory Density approach.

The paper is structured as follows: in Section 2, we present the microscopic model describing the evolution of the disease within a single individual. In Section 3 we introduce the equations governing the macro-scale population, while the mesoscopic ones are given in Section 4. We illustrate our work in Section 5 with numerical simulations and compare them to real data. We finally discuss our work in Section 6.

Figure 1: Schematic of the micro- and macroscopic models described in sections 2 and 3. A, Microscopic model. The state variable $V_{i}$ of an individual $i$ evolves in time $t$ (bottom trace); its crossings of the threshold $V^{T}$ determine the onsets of disease, the infection time moments, when the drive $D(t)$ appear, which describes virus production and immune response (top). The time since the last infection $t^{*}$ is zeroed at the onset of the disease. B, Macroscopic model. At time $t$, individuals with different $t^{*}$ are distributed with the density $\rho$ in the $t^{*}$-space (bottom). With time, the density transports rightwards (as $t^{*}$ increases with $t$) and undergoes a return flux to the infection state $t^{*}=0$ according to the hazard function $H$, so the dashed arrows represent subsequent infections. The hazard depends on the mean state variable $U$, which is similar to $V_{i}$ in A but is extrapolated over the threshold.

2 Microscopic model of single individuals

We consider a population of $N\gg 1$ individuals, among which an infectious disease is spreading. We do not consider birth and death processes, and consequently, we assume the population to be constant (although this assumption is not critical for our framework). We assume that the risk of infection for an individual depends on two state variables $V_{i}$ and $t^{*}_{i}$ (see Table 1 for a summary of the variables and parameters of the model). The variable $V_{i}$ represents the interplay of the viral load and the immune response within individual $i$. We call this variable viral state. The dynamics of $V_{i}(t)$ are governed by the interaction of multiple processes (Fig.1A, bottom). First, due to the immune system, the number of viral particles relaxes with time constant $\tau$. Second, it grows proportionally to the fraction of infected individual $I(t)$, with a coefficient of proportionality $k$, which is assumed to be proportional to the basic reproduction number $R_{0}$ of the corresponding disease. This modeling assumption is qualitatively close to the underlying homogeneous mixing assumption of classical ODE models, whereby any individual is equally likely to meet with any other individual in the population. Third, it mimics the time course of the disease composed of a rapid increase of the viral load at disease onset and a subsequent decrease caused by the immune response. This time course is driven by the function $D(t_{i}^{*})$, started at each onset of the disease (Fig.1A, top). Note that, since $V_{i}$ represents the interplay between the virus and the immune system, it can become negative at the negative phase of $D(t_{i}^{*})$.

Variable	Notation	Units
$t$	Time	days
$t^{*}$	Time since last infection, or “infection age”	days
$V_{i}$	Interplay of viral load and immune response for individual $i$	v-units
$U$	Mean interplay of viral load and immune response	v-units
$\rho$	Density of individuals distributed in $t^{*}$-space	people/day
$\nu$	Population rate of new infections	people/day
$I$	Fraction of infectious individuals in the population	nondim.
$\tau$	Time scale of passive immune response	days
$\tau_{r}$	Time scale of disease progress and cure	days
$\tau_{I}$	Infectious period	days
$\tau_{R}$	Time from infection to loss of immunity (“Refractory period”)	days
$D$	Time course of the interplay for disease and treatment	v-units/day
$N$	Total number of people in the population	people
$N_{c}$	Cumulative cases in the population	people
$R_{0}$	Basic Reproduction Number	nondim.
$k$	Average individual rate of viral load	v-units/day
$H$	Hazard; $H\Delta t$ is the probability of infection in the interval $(t,t+\Delta t)$	people/day
$V^{T}$	Threshold of infection	v-units
$\sigma$	Noise amplitude	v-units

Table 1: Notation for the micro- and macroscopic models described in sections 2 and 3.

The state variable $t_{i}^{*}$ measures the time relative to the last infection time (onset of the disease) of individual $i$. These onsets define an ordered sequence $t_{1,i}<t_{2,i}<t_{3,i}<\dotsb$ of infection times for individual $i$. The continuous states of the epidemic in the $t^{*}$-space may be roughly interpreted in terms of discrete SEIRS stages (Fig.1A). To this end, we may partition the population with respect to an ongoing epidemic in the following compartments: $S$, individuals who are Susceptible to the disease; $E$, individuals who have been Exposed to the disease but are not yet infectious; $I$, Infected individuals who may spread the disease; and $R$, Recovered individuals who are temporarily immune from infection. According to these definitions, we can loosely define the stages in the $t^{*}$-axis. When $V_{i}(t)$ overcomes a threshold $V^{T}$ (in the case of Gaussian noise, see below) or an infection event is triggered by a probabilistic “soft” threshold (in the case of escape noise), the individual becomes infected and consequently belongs to compartment $I$. During the infectious phase, $V_{i}(t)$ rises above the threshold because of the viral production during disease and then decreases because of the immune response. The $I$ state lasts for a fixed duration $\tau_{I}>0$. After the time $\tau_{I}$ has elapsed since the onset of infection, the individual is effectively immune to the disease and consequently belongs to compartment $R$. The sojourn in this compartment lasts until the end of the recovery period $\tau_{R}>\tau_{I}$, measured from the onset of infection. After this time, when $V_{i}$ is low, the individual becomes Susceptible and then Exposed until the next infection when the individual becomes Infected and infectious again. Note that only the Infectious and Recovered states are strictly defined by the interval $0<t^{*}<\tau_{I}$ and $\tau_{I}\leq t^{*}<\tau_{R}$, respectively. However, the duration of these states does not directly affect the evolution of the individual’s state $V_{i}$, but instead determines the periods when the individual affects the population and is in an “absolute refractory period”, respectively.

Taken together, the dynamics for single individuals $i$, $i=1,\dotsc,N$, between two subsequent infection times $t_{k,i}$ and $t_{k+1,i}$ is modeled as

$$\frac{\text{d}V_{i}}{\text{d}t} =-\frac{V_{i}}{\tau}+k(t)I(t)+D(t_{i}^{*})+\sigma(I(t))\leavevmode\nobreak\ \xi_{i}(t),$$ (1)

$$\frac{dt_{i}^{*}}{dt} =1,$$ (2)

with the fraction of infected individuals

$$I(t) =\frac{N_{c}(t)-N_{c}(t-\tau_{I})}{N},$$ (3)

where $N_{c}(t)\in\mathbb{N}$ is a global counting variable that keeps track of the cumulative number of cases in the population up to time $t$. Clearly, this counting variable is constant when no new infection occurs. At the infection times $t_{k,i}$, however, the overall dynamics is supplemented with an additional jump condition (reset):

$$V_{i}(t_{k,i}^{+})=V^{T},\quad t_{i}^{*}(t_{k,i}^{+})=0,\quad N_{c}(t_{k,i}^{+})=N_{c}(t_{k,i}^{-})+1,$$ (4)

where $t_{k,i}^{-}$ and $t_{k,i}^{+}$ denote the left- and right-sided limits, respectively. In Eq. (3), we indicate with $\tau_{I}$ the average time spent in the infectious state or average duration of the infectious window. We note that Eq. (3) can be easily generalized to permit a more fine-grained, graded account of infectiousness depending on infection age $t^{*}$ (see Discussion, Eq. (16)). Furthermore,

$$D(t^{*}) =4a\exp\left(-\frac{t^{*}}{\tau_{r}}\right)-a\exp\left(-\frac{t^{*}}{4\tau_{r}}\right)$$ (5)

is a response function that drives the time course of $V_{i}$ during disease. The parameter $a$ is the severity of the disease, and $\tau_{r}$ is the time of disease progress and cure. The last term in Eq. (4) represents an independent Gaussian white noise with mean $\mathbb{E}[\xi(t)]=0$, auto-covariance function $\mathbb{E}[\xi(t)\xi(t^{\prime})]=\delta(t-t^{\prime})$ and strength $\sigma(I)$ ($\delta(\cdot)$ denotes the Dirac delta function). As explained below, this term is only present in our model variant with Gaussian noise and is absent in the variant with escape noise where we set $\sigma(I)\equiv 0$.

As initial conditions of the model, we assume that before time $t=0$ there were no cases, $N_{c}(t)=0$ for $t<0$ and that at time $t=0$ a fraction $I_{0}$ of the population gets infected, i.e. $N_{c}(0)=I_{0}N$. Thus the state variables for the $I_{0}N$ initially infected individuals are $V_{i}(0)=V^{T}$ and $t_{i}^{*}(0)=0$, while the state variables of all other (non-infected) individuals (the remaining $(1-I_{0})N$ individuals in the population) is initially $V_{i}(0)=0$ and $t_{i}^{*}(0)=+\infty$.

Importantly, it remains to define how the infection events are triggered. We consider two different variants to define these events, and hence the infection times $t_{k,i}$, corresponding to two different ways of including stochasticity in the model. We refer to these variants as the white- and escape-noise cases, respectively.

Model with escape noise

To introduce stochasticity that phenomenologically captures the influences of various sources of noise, we model the infection events of an individual $i$ by a probabilistic risk of becoming infected depending on the individual’s current state variables $V_{i}(t)$ and $t^{*}_{i}(t)$. Mathematically, we keep the dynamics of $V_{i}$ between infection events deterministic (i.e. setting $\sigma=0$) and generate infection events stochastically by a state-dependent hazard rate

$$H_{i}(t)=H(V_{i}(t),t_{i}^{*}(t)).$$ (6)

The hazard rate means that an infection of individual $i$ occurs in a small time step of length $\Delta t$ with conditional probability $H_{i}(t)\Delta t$ given the individual’s current state $(V_{i}(t),t^{*}_{i}(t))$. More precisely, we have the following conditional probability:

$$\text{Pr}\left\{\text{individual $i$ gets infected in }(t,t+\Delta t)\;|\;V_{i}(t),t_{i}^{*}(t^{-})\right\}=H_{i}(t)\Delta t+o(\Delta t)$$ (7)

as $\Delta t\to 0$. For a concrete choice of the hazard rate in simulations, we use a rectified power function with absolute refractory period $\tau_{R}$:

$$H(V,t^{*})=c\max(0,V/V^{T})^{m}\mathbbm{1}_{t^{*}\geq\tau_{R}},$$ (8)

where $m>0$. We illustrate the effect of varying $m$ in Figure 2.

Model with white Gaussian noise

An alternative way to model noise is to consider an additive white-noise drive in the dynamics of $V_{i}$. This “diffusive noise” is modeled by the term $\sigma(I(t))\leavevmode\nobreak\ \xi_{i}(t)$ in Eq. (1). The stochastic term may reflect both the fluctuations of the number of viral particles of a given type around a particular individual and the fluctuations of activity of the immune system of this individual. To prevent that in the absence of a viral infection in the population an infection occurs spontaneously by noise, we enforce the noise strength to be zero when $I(t)$ is zero by choosing the specific function $\sigma(I)=\hat{\sigma}\mathbbm{1}_{I>0}$. In the white-noise case, infection events are triggered when the viral state $V_{i}(t)$ crosses a certain threshold $V^{T}$, provided that at least an absolute refractory time $\tau_{R}$ has passed since the last infection. Hence, the infection times $t_{k,i}$ are defined by the condition

$$V(t_{k,i}^{-})\geq V^{T},\quad t^{*}(t_{k,i}^{-})\geq\tau_{R}.$$ (9)

Note that if this condition is fulfilled at some time $t_{k,i}^{-}$, the reset condition Eq. (4) ensures that the above threshold condition no longer holds at $t_{k,i}^{+}$, i.e. immediately after this time, as it should be.

The above equations conclude the definition of the microscopic model. Given the following application of the refractory-density (RD) method, we note that the case of white Gaussian noise can be approximately mapped to the case of escape noise [7, 26]. In this case, the noise term in Eq. (1) is omitted and the condition Eq. (9) is substituted by (7), where $H_{i}(t)$ is now given by the following hazard function (10). The hazard rate $H(t,t^{*})$ depends on the actual state $V_{i}$, its rate of change $\text{d}V_{i}/\text{d}t$, the threshold $V^{T}$, and the noise amplitude $\sigma$, i.e., $H(t,t^{*})=H(V_{i}(t),\dot{V}_{i}(t),t^{*}_{i}(t);\sigma(I(t)))$. In [7], this function was found as an approximate solution of a first-time passage problem based on the Kolmogorov-Fokker-Planck equation:

$$H(U,\dot{U},t^{*},\sigma) =\{A(T)+B(T,\dot{U},\sigma),\leavevmode\nobreak\ \leavevmode\nobreak\ \hbox{if}\leavevmode\nobreak\ \leavevmode\nobreak\ t^{*}>\tau_{R},\leavevmode\nobreak\ \leavevmode\nobreak\ \sigma>0;\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ 0,\leavevmode\nobreak\ \hbox{otherwise}\},$$ (10)

$$A(T) =\frac{1}{\tau}\leavevmode\nobreak\ \exp(0.0061-1.12T-0.257T^{2}-0.072T^{3}-0.0117T^{4}),$$

$$B(T,\dot{U},\sigma) =\frac{2}{\sigma\sqrt{\pi}}\bigg{[}\dot{U}\bigg{]}_{+}\frac{\exp(-T^{2})}{1+\hbox{erf}(T)},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ T=\frac{V^{T}-U}{\sigma}.$$

Here and in the following, $[x]_{+}=\max(0,x)$ denotes the rectified linear function. We refer to [26] for a slightly more elaborate approximation of diffusive noise by escape noise.

3 Macroscopic model of a population

For a large number of individuals $N$, the dynamics of the microscopic model can be either evaluated with a Monte-Carlo simulation or obtained much more efficiently by integrating corresponding macroscopic population equations. To this end, we apply the refractory-density (RD) approach [27] to the epidemiological model described at an individual level in section 2. According to this approach, the state variables of a single “particle”/individual are parameterized by the “age” $t^{*}$ (with the total derivative in time substituted by the sum of partial derivatives, i.e. $\text{d}/\text{d}t=\partial/\partial t+\partial/\partial t^{*}$), and the population is characterized by the density of particles in the one-dimensional, semi-infinite $t^{*}$-space, i.e. $\rho(t,t^{*})$ (Fig.1B). In our case, in line with other age-dependent models in mathematical epidemiology, the “age” is the time elapsed since infection. This means that at any time $t$, the number of individuals in the population who became infected between $t^{*}_{1}$ and $t^{*}_{2}$ time units ago are represented by

$$\int_{t^{*}_{1}}^{t^{*}_{2}}\rho(t,t^{*})\text{d}t^{*}.$$

At the macroscopic level, we define the rate of new infections

$$\nu(t)=\lim_{\Delta t\to 0}\lim_{N\to{}\infty}\frac{N_{c}(t+\Delta t)-N_{c}(t)}{N\Delta t}.$$ (11)

We also introduce a function $U(t,t^{*})$ as follows: for the model with escape noise, $U(t,t^{*})$ is the unique function that maps the infection age $t^{*}_{i}(t)$ to the state variable $V_{i}(t)$, i.e. $V_{i}(t)=U(t,t^{*}_{i}(t))$. For the model with Gaussian white noise, there is no such deterministic map. However, we can first map the model to the approximate escape-noise model, Eq. (10), for which we can define again the unique function $U(t,t^{*})$. Intuitively, this function can be interpreted as the variable $V_{i}(t)$ averaged across realizations of the Gaussian-white noise for individuals $i$ that had the time of last infection equal to $t^{*}$, i.e., the same $t^{*}$. The equation for $U$ then follows from the corresponding averaging of Eq. (1). The bottom traces in Fig.1B illustrate distributions of $U$ and $\rho$ in $t^{*}$-space.

The RD model corresponding to Eqs. (1)-(4) consists of two transport equations for $\rho(t,t^{*})$ and $U(t,t^{*})$, and two integrals for $\nu(t)$ and $I(t)$:

$$\frac{\partial\rho}{\partial t}+\frac{\partial\rho}{\partial t^{*}} =-\rho\leavevmode\nobreak\ H(U(t,t^{*}),t^{*}),$$ (12a)

$$\frac{\partial U}{\partial t}+\frac{\partial U}{\partial t^{*}} =-\frac{U}{\tau}+k(t)I(t)+D(t^{*}),$$ (12b)

$$\nu(t) =\int_{0}^{\infty}\rho(t,t^{*})\leavevmode\nobreak\ H(U(t,t^{*}),t^{*})\,dt^{*}+I_{0}\delta(t),$$ rate of new cases (12c)

$$I(t) =\int_{t-\tau_{I}}^{t}\nu(s)\,ds,$$ fraction of infected people (12d)

The boundary condition for $\rho$, resulting from the conservation of people in the population, is

$$\rho(t,0)=\nu(t).$$ (13)

The boundary condition for $U$ corresponding to the reset in the microscopic model, Eq. (4), is $U(t,0)=V_{T}$. The initial condition for the density $\rho$ corresponding to an initial fraction of infected people $I_{0}$ at $t=0$ and to a situation, where individuals never encountered this viral infection before time $t=0$ (formally, $t^{*}_{i}(0^{-})=\infty$) is

$$\rho(0,t^{*})=\delta(t^{*})I_{0}+\delta(t^{*}-\infty)(1-I_{0}).$$

Furthermore, the corresponding initial condition for $U$ can be taken as the steady-state solution of Eq.(12b) (with $\partial U/\partial t=0$ and $I(t)=0$ for $t<0$):

$$U(0,t^{*})=V^{T}\exp{(-t^{*}/\tau)}+\int_{0}^{t^{*}}e^{-s/\tau}D(t^{*}-s)\,ds.$$

In the density equation (12a), the density diminishes because of the sink term $-\rho H$ but the total mass $\int_{0}^{\infty}\rho(t,t^{*})\,dt^{*}=1$ is conserved at all times because of the boundary condition, Eq. (13). This boundary condition acts as a source term at $t^{*}=0$ given by the rate of new cases $\nu(t)$ that exactly balances the integrated sink term, Eq. (12c). We note that this conservation law reflects the fact that we neglect mortality, and focus on the transmission mechanism.

Furthermore, in the density equation (12a), the risk of illness is evaluated by the hazard function $H$. In Eq. (12), we used the hazard function for the case of escape noise. For the case of white noise, the hazard function $H(U(t,t^{*}),t^{*})$ should be replaced by $H(U(t,t^{*}),\dot{U}(t,t^{*}),t^{*},\sigma(I(t)))$, where $H$ is given by Eq. (10) and $\dot{U}(t,t^{*})=(\partial_{t}+\partial_{t^{*}})U(t,t^{*})$ is given by the right-hand-side of Eq. (12b).

Once again, note that Eq. (12d) permits a simple generalization to graded infectiousness depending on the course of the infection (see Discussion, Eq. (16)).

4 Mesoscopic model of a finite-size population

At the intermediate mesoscopic scale, the size of (sub-)populations is finite. The finite size causes fluctuations in the fraction of infectious people $I(t)$, which may yield significant finite-size effects through the interaction of finite-size noise and nonlinear population dynamics. In the case of a finite size population consisting of $N\in\mathbb{N}$ individuals, $N\gg 1$, we apply the corresponding theory from [25, 28], which yields a stochastic generalization of the macroscopic dynamics given in the previous section. To this end, we introduce the pseudo-density $\rho(t,t^{*})$ in terms of the survivor function $S(t,t^{*})$ as

$$\rho(t,t^{*})=S(t,t^{*})\leavevmode\nobreak\ \nu_{N}(t-t^{*}).$$ (14a)

Here, $$S \nu_{N}$$ and are given as the solution of the following system of stochastic partial differential equations, generalizing Eq. (12)

$$\frac{\partial S}{\partial t}+\frac{\partial S}{\partial t^{*}} =-S\leavevmode\nobreak\ H(U),\qquad S(t,0)=1,$$ (14b)

$$\frac{\partial U}{\partial t}+\frac{\partial U}{\partial t^{*}} =-\frac{U}{\tau}+u(t)+k\leavevmode\nobreak\ I(t)+D(t^{*}),\qquad U(t,0)=V_{T},$$ (14c)

$$I(t) =\int_{t-\tau_{I}}^{t}\nu_{N}(s)\,ds,$$ (14d)

$$\nu(t) =\left[\int_{0}^{\infty}\rho(t,t^{*})\leavevmode\nobreak\ H(U(t,t^{*}))\,dt^{*}+\bar{H}(t)\left(1-\int_{0}^{\infty}\rho(t,t^{*})\,dt^{*}\right)\right]_{+},$$ (14e)

$$\nu_{N}(t) =\nu(t)+\sqrt{\frac{\nu(t)}{N}}\xi(t),$$ (14f)

with

$$\bar{H}(t)=\frac{\int_{0}^{\infty}H(U(t,t^{*}))(1-S(t,t^{*}))\rho(t,t^{*})\,dt^{*}}{\int_{0}^{\infty}(1-S(t,t^{*}))\rho(t,t^{*})\,dt^{*}},$$ (14g)

where $\xi(t)$ is the white Gaussian noise of unity amplitude. We remark that, for all $t,t^{*}$, $S(t,t^{*}),\in[0,1]$; hence, the function $\bar{H}(t)$ is non-negative for all $t$. Furthermore, we note that $\nu_{N}(t)$ must be regarded as an abstract quantity that mathematically only makes sense as a distribution, i.e. if it is integrated against some test function over time. Specifically, the empirical rate of new cases $\hat{\nu}_{N}$ measured over some finite time step $\Delta t$ (such that the expected number of new cases $N\nu(t)\Delta t\gg 1$) would be

$$\hat{\nu}_{N}(t)=\left[\int_{t}^{t+\Delta t}\nu_{N}(s)\,ds\right]_{+}.$$ (15)

While in practice it is extremely rare that the integral takes negative values for large $N$, we added the rectification $[\cdot]_{+}$ to enforce a non-negative rate. The same statement also holds true for Eq. (14e).

Likewise, $\rho$ is not a normalized probability density [25, 28]. Therefore, it is no longer the strict interpretation of the empirical density of the times since the last infection but $\rho(t,t^{*})dt^{*}$ must be regarded as an abstract measure. In comparison with the macroscopic model, this system of equations is stochastic. Its solution corresponds to a single realization of the noise $\xi(t)$. When $N\xrightarrow{}\infty$, the solution converges to the solution of the macroscopic model. The precise numerical method employed for these simulations is described in [25] (see also the provided simulation code in Python).

5 Simulations

Model with escape noise

A first simulation of epidemic waves is shown in Fig. 2. The onset of an epidemic is modeled by a short pulse of magnitude $I_{0}$ of the rate of new infections $\nu(t)$ at time $t=0$ reflecting the appearance of contaminating patients. These conditions result in oscillations of the rate of cases $\nu(t)$ and the fraction of potentially infected population $I(t)$. These oscillations tend to relax. The dynamics of a single individual state $V_{i}(t)$ (Fig. 2A, next to bottom panel) follows the epidemic waves of $\nu(t)$ and $I(t)$. Each peak of $V_{i}(t)$ reflects the virus reproduction during the disease. The decreasing phase reflects the immune response, after which an individual becomes again susceptible and might get infected during the next contaminating flux $k(t)I(t)$. Each spike of the rate of cases $\nu(t)$ leads to a correspondingly shaped distribution of the density in the $t^{*}$-space.

A                                                                                      B

Figure 2: Escape noise in the form of (8). Simulation of the response to the rapid change of the rate of viral load $k$ from 0.6 to 2.4 and the short infecting pulse $I_{0}(t)$, at $t=100$. A, top to bottom: $I_{0}(t)$, $k(t)$, the fraction of infected individuals $I(t)$, the state variable $V_{i}(t)$ of a representative individual, and the rate of new cases $\nu(t)$. For $I(t)$ and $\nu(t)$, the solutions obtained with both macro- (black lines) and microscopic equations for $N=20000$ (blue lines) are shown. B, the distribution in the $t^{*}$-space of the mean state variable $U$, the source term $\rho H$, and the density $\rho$, in 4 time moments. Parameters: $m=1$, $c=0.015$ days^-1, $\tau=50$, $\tau_{I}=10$, $\tau_{r}=1$, $\tau_{R}=20$, $a=100$, and $V^{T}=2$.

The population waves move in the $t^{*}$-space from $t^{*}=0$ towards the region where $U$ is about to cross the threshold $V^{T}$ and form the peaks of $\rho H$, which is the distribution of falling sick individuals (Fig. 2B). Since the integral over $\rho H$ determines the rate of cases $\nu(t)$, the increase of the $\rho H$ determines the next peak of $\nu(t)$ and $I(t)$, and so on.

We assume that at $t=t_{1}=700$ the basic reproduction rate drops 4 fold, which might be due to any kind of containment measures. As a result, the epidemic wave generation stops through a significant decrease in the rate of cases $\nu(t)$ (Fig. 2A).

For the escape noise in the form of (8), the epidemic depends on the steepness $m$ of the hazard function (Fig. 3). With the increase of $m$ from 0.5 to 1, the epidemic oscillations increase. However, it is non-trivial that in the case of more steep $H$ with $m=2$ the epidemic shows only one peak. In this case, the rapid illness of the whole population results in a short splash of contamination $k(t)I(t)$ that occurs at the recovery phase of the population when no one is susceptible.

Figure 3: Escape noise in the form of (8) with $m=$ 0.5, 1 and 2. The remaining parameter values are as in Fig. 2.

Seasonality

One of the experimentally observed features of epidemics is its seasonal oscillations, as illustrated in Fig. 4 by data from [21]. In Fig. 4, we illustrate the case of seasonal change of $k$ through a step function alternating between $k=2.4$ and $k=1.2$ every 180 days. In response to a meander-like $k$, we observe a complex pattern of waves with various amplitudes.

A

B

Figure 4: Epidemiology of Seasonal Coronaviruses. A, Data from [21]; specifically, the monthly prevalence of seasonal coronaviruses (sCoVs) detected among patients with respiratory illness virologically tested in NHS Greater Glasgow and Clyde, Scotland, United Kingdom. B, The case of oscillating $k$. Escape noise in the form of (8). The remaining parameter values are as in Fig. 2.

Model with white Gaussian noise

When simulating the system with the white Gaussian noise (Fig. 5), we observed a similar behavior. The onset of an epidemic is again modeled by a short pulse of the rate of new infections $\nu(t)$ of magnitude $I_{0}$ at time $t=0$ reflecting the appearance of infectious individuals. At the same time, the viral load noise jumps ($\sigma$ changes from $0.5$ to $2$). These conditions result in sustained oscillations of the rate of cases $\nu(t)$ and the fraction of potentially infected population $I(t)$. Again, the dynamics of a single individual state $V_{i}(t)$ (Fig. 5A, next to bottom panel) follows the epidemic waves of $\nu(t)$ and $I(t)$. When at $t=t_{1}=900$ the basic reproduction rate $k$ drops due to assumed anti-epidemic measures, the epidemic wave generation stops.

A                                                                                      B

C                                                                   D

Figure 5: White noise. Simulation of the response to the rapid change of the rate of viral load $k$ from 3 to 1.2 and the noise amplitude $\sigma(t)$ from 0.5 to 2, at $t=t_{0}=0$. A, top to bottom: $k(t)$, $\sigma(t)$, the fraction of potentially infected population $I(t)$, the state variable $V_{i}(t)$ of a representative individual, and the rate of new cases $\nu(t)$. For $I(t)$ and $\nu(t)$, the solutions obtained with both macro (black lines) and microscopic equations for $N=20000$ (blue lines) are shown. B, the distribution in the $t^{*}$-space of the mean state variable $U$, the source term $\rho H$, and the density $\rho$, in 4 time moments. C, Response to periodically changing $k$ mimicking change of seasons ($V^{T}=4$). D, Response to periodically changing $k$ mimicking change of seasons ($V^{T}=3.7$). Parameters: $\tau=50$, $\tau_{I}=10$, $\tau_{r}=1$, $\tau_{R}=20$, $a=50$, $V^{T}=4$, and $\sigma=2\sqrt{\pi}$.

Mesoscopic model

In order to consider a finite number of individuals in a population, we apply the mesoscopic model described in section 4, based on Eqs.(14a)-(14g). The solutions are different for different $N$ (Fig.6). The solutions for small $N$ differ significantly for different realizations of the noise. The model in the limit $N\xrightarrow{}\infty$ is equivalent to the macroscopic model. Hence, the solutions for large $N$ converge for different realizations and approach the macroscopic model solution shown in Fig. 2.

Figure 6: Mesoscopic model (8) for different $N=200$, 2000, and 20000, and different realizations of noise. The parameters are as in Fig. 2. The solutions with $N\xrightarrow{}\infty$ converge to the macroscopic model solution shown in Fig. 2.

6 Discussion

In the present work, we have proposed the description of an epidemic that stems from a microscopic model of an individual and results in the meso- and macroscopic models of the entire population. The microscopic evolution is described by the dynamics of two state variables for each individual, the interplay between the viral load and immune response $V_{i}$, and the time since the last infection $t^{*}_{i}$. The corresponding mesoscopic and macroscopic models are then given by deterministic and stochastic refractory-density equations, respectively. We have shown the consistency between the models. The presented simulations remind of epidemic “waves”, for instance, the ones reported in [21]. One of our basic ideas to describe the epidemic in terms of the evolution of viral load depending on the time since the last infection is supported by the experimental data from [16], where the viral load was measured in patients with SARS-CoV-2 as a function of the time since infection. Our multi scale modeling framework allows us to take into account such microscopic parameters measured in individuals in order to constrain the parameters of the meso- or macro-scale models.

The approach and the model that we presented in this manuscript could naturally be extended and generalized. We conclude the manuscript with a few promising leads for future research, which we identified during the writing of this work.

6.1 Treatment

In simulations, we varied the individual rate of viral load $k(t)$, which accounted for the seasonality. This factor can also reflect such social measures as wearing masks and isolation. In contrast, medical treatment of patients is a function of the time since infection $t^{*}$, and it changes the course of the disease. Naturally, this factor can be taken into account through modification of the shape of the function $D$ entering the equation for the interplay of viral load and immune response, Eq. (1).

6.2 Modelling infectiousness depending on infection age

It has been shown that the infectiousness of SARS-CoV-2 depends on the time since the last infection [16]. This can be easily modeled in our framework by generalizing Eq. (3), Eq. (12d) and Eq. (14d) for the micro-, macro- and mesoscopic fractions of infected people to

$$I(t)=\frac{1}{N}\int_{0}^{t}\kappa(t-s)\,dN_{c}(s),\quad I(t)=\int_{0}^{t}\kappa(t-s)\nu(s)\,ds\quad\text{and}\quad I(t)=\int_{0}^{t}\kappa(t-s)\nu_{N}(s)\,ds,$$ (16)

respectively. The kernel $\kappa:\mathbb{R}^{+}\to[0,1]$ describes the infectiousness depending on the infection course. Our equations (3),  (12d) and Eq. (14d) are recovered as the special case $\kappa(t^{*})=\mathbbm{1}_{(0,\tau_{I})}(t^{*})$, allowing only two possible states – infectious and non-infectious. Our infectious period $\tau_{I}=10\leavevmode\nobreak\$days well correspond to the half-duration of cell culture infectivity measured in [16] for SARS-CoV-2.

6.3 Network of multiple populations

Certain types of heterogeneity (e.g. different social communities, age groups, etc.) and spatial structure (cities) can be treated by splitting a large heterogeneous population into smaller homogeneous subpopulations, again in analogy to RDE application in neuroscience, for instance, to simulate cortical neuronal populations [9]. Moreover, the methodology can be extended to interacting populations, such as in the scenario of epidemic transmission between countries, consideration of a network of communities, similar to the so-called multi-group approach (see e.g. [1, 23]). Since the multi-group approach naturally operates with smaller-size subpopulations, the proposed mesoscopic approach could be applied in this case.

6.4 Multiple internal variables

The proposed derivation of the population model from that for an individual applies not only to the 1-D individual model described by a single ODE for the viral state $V_{i}$, but also to multidimensional cases, in which one may be interested in including additional characteristics of individuals and/or of populations, and keep track of these as the disease spreads. In this way, additional internal variables can be introduced, describing, for instance, activation and/or inactivation of some processes in the immune system. This concept resembles the consideration of Hodgkin-Huxley-like neurons in neuroscience [7]. If an individual is characterized by $n$ ODEs, the population model comprises $n+1$ partial differential equations (PDEs) for the state variable and the density. These equations remain one-dimensional, with time since the infection $t^{*}$ is the sole independent variable, besides $t$. Consequently, this approach maintains computational efficiency even for multi-dimensional individuals. Leveraging this advantage, more intricate models can be developed.

In comparison with the above-mentioned kinetic theory [20] that results in the Fokker-Planck equation describing the evolution of the population in the phase space of a state variable, the proposed approach considers the evolution in the phase space of the time since infection. Once the microscopic models are identical, these two macroscopic approaches would result in very similar numerical simulations, at least in the case of white Gaussian noise [27]. In the more general case of non-scalar state variables, the Fokker-Planck equation becomes multidimensional, whereas the transport equations of the refractory density approach remain 1-D PDEs that are effectively solvable.

6.5 Age-dependence

The epidemics to a different extent affect people of different ages [16]. Our approach allows the consideration of the population distributed by age $\hat{a}$ via an age-dependent hazard function $H(V,t^{*},\hat{a})$. In this case, the main variables have to be parameterized by $\hat{a}$ as $\rho=\rho(t,t^{*},\hat{a})$ and $U=U(t,t^{*},\hat{a})$. Assuming that the hazard function changes much slower with age $\hat{a}$ than with the time since the last infection (i.e. $|\partial H/\partial\hat{a}|\ll|\partial H/\partial t^{*}|$), the age $\hat{a}$ can be treated quasi-statically as a parameter and not as an independent dynamical variable. Under this plausible assumption, the equations of the macroscopic model remain to be 1-D transport equations. However, the fraction of infected individuals $I(t)$ would be an integral over $\hat{a}$. This parametrization would allow the consideration of different disease time courses, etc. The infectivity $\kappa(t^{*},\hat{a})$ can also be dependent on $\hat{a}$ because usually people interact stronger within their age groups.

6.6 Asymptomatic spread

We remark that, for our model, we are interested in the ability of an individual to spread the disease, rather than showing symptoms; hence, the infected and infectious population might be generalized to include asymptomatic individuals, thus allowing to adapt our construction to more complex compartmental models such as the ones presented in [4, 5, 15, 22, 23]. This may be done in multiple ways, depending on which characteristics of symptomatic and asymptomatic spread one wishes to capture in their model.

6.7 Heterogeneous severity

It could be of biological interest to assume a heterogeneous setting, e.g. the severity $a$ in Eq. (5) could vary among individuals depending on the vaccine status, stronger or weaker immune system, etc.; whereas we can assume that $\tau_{r}$ is approximately constant and depends only on the specific disease. However, this generalization would carry a non-negligible increase in the computational cost of our simulations, and make mean field approximations which we perform in section 3 and 4 considerably more difficult. One potential strategy to deal with such heterogeneity would be to split the population into several, approximately homogeneous subpopulations and use a multi-population version of our model as described above. Importantly, in this scenario, we could take full advantage of the mesoscopic model because it remains valid for small sub-populations. We leave such a generalization to heterogeneous systems as a promising outlook for future work.

6.8 Complexity and parameters

Despite more mathematical complexity in comparison with the classical SEIRS models (PDEs versus ODEs), the proposed model has a similar number of parameters. The SEIRS model has 5 parameters ($\beta$, $\sigma$, $\gamma$, and $\xi$, in traditional notations, and the time scale). Instead, in the case of omitted explicit course of the recovery $D(t)$, the proposed model has also only 5 parameters: $\tau$, $\tau_{I}$, $k$, $\sigma$, and $V^{T}$.

6.9 Inverse problem

The model may be used in future studies for solving reverse problems of finding the parameters from statistical data for a certain disease. For this purpose, one should have a sufficient set of experimental data that would reflect the effects of each of the parameters. For instance, the effect of $k$ can be studied with the effect of masks; the effect of $\tau_{I}$ with the effect of care; the effect of $V^{T}$ with the effect of provisional stimulation of the immune system of the entire population; the effect of $\sigma$ with a selection of more homogeneous subpopulation; the effect of $\tau$ with the selection of faster (or slower) recovering subpopulation, and so on.

6.10 Limitations

In this first work on the topic, we made several simplifying assumptions. One of the most important is that our model assumes a strictly conserved population: no removal, no mortality, closed community: no external influx, no out-flux. It would be of interest to consider, instead, a similar model in which the population is allowed to vary in time, for any of the reasons listed above.

Declaration of Competing Interest

The authors have no competing interest to declare.

Declaration of equal contribution

All authors contributed equally to this work.

Code availability

The code used to numerically solve the models presented in this manuscript is available at
https://github.com/schwalger/refracdens_epidem.

Appendix A Macroscopic vs microscopic models, white Gaussian vs. escape noise

To illustrate the agreement between microscopic and mesoscopic models as well as escape noise and Gaussian white noise, we consider a slightly simplified setup. Whereas the main model equation (1) includes the term $k(t)I(t)$, which is recurrently dependent on the activity of the entire population $\nu(t)$, here we consider a simplified problem with a step-wise input signal instead of the term $k(t)I(t)$. This case can be interpreted as a rough approximation of an epidemic problem, where the time course of the epidemic is assumed to be known and shaped as a step function. In this interpretation, we are interested in the probabilistic behavior of an individual, with $\nu(t)$ to be a probability for this individual to get ill. We simulate this response for the case of white noise and the escape noise in the form of Eq. (10) (Fig.7).

Figure 7: Solutions of a test with a stepwise term instead of $k(t)I(t)$ in the Eqs. (1) and (12d), obtained with the microscopic model ($N=2000$) with the white Gaussian noise and the escape noise described with Eq. (10), and the microscopic model. The value of $k(t)I(t)$ is 0 for $t\in[0,100)$ and $t>500$ and $0.4$ for $t\in[100,500]$.

From the methodical point of view, Fig. 7 shows the testing comparison of the macroscopic model with the microscopic models with white Gaussian and escape noise. The microscopic solutions converge to the macroscopic ones.