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

Anton Chizhov, Laurent Pujo-Menjouet, Tilo Schwalger, Mattia Sensi

► To cite this version:

Anton Chizhov, Laurent Pujo-Menjouet, Tilo Schwalger, Mattia Sensi. A refractory density approach to a multi-scale SEIRS epidemic model. Infectious Disease Modelling, 2025, 10 (3), pp.787-801. 10.1016/j.idm.2025.03.004 . hal-05000123

HAL Id: hal-05000123

https://hal.science/hal-05000123v1

Submitted on 24 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# A refractory density approach to a multi-scale SEIRS epidemic model

<figure><img src="image_2.png" /></figure>

Anton Chizhov<sup>a, b, c</sup>, Laurent Pujo-Menjouet<sup>d, *</sup>, Tilo Schwalger<sup>e, f</sup>, Mattia Sensi<sup>b, g</sup>

## ARTICLE INFO

**Article history:**

Received 26 July 2024

Received in revised form 7 March 2025

Accepted 10 March 2025

Available online 13 March 2025

Handling Editor: Dr Yijun Lou

**Keywords:**

Refractory density

Partial differential equations

Age-structured model

Time since last infection

Finite-size fluctuations

## 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.

© 2025 The Authors. Publishing services by Elsevier B.V. on behalf of KeAi Communications Co. Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

## 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 (Kermack & McKendrick, 1927) for a wide variety of mathematical models in epidemics; see e.g. (Iannelli & Milner, 2017; Alfaro-Murillo et al., 2019; Anderson, 1988; Brauer and Watmough, 2009, d'Onofrio et al., 2023; Li et al., 2019). Taking into account the time since the last infection (TSLI) is extremely convenient to introduce a temporal

<footnote><sup>a</sup> Institute for Theoretical Physics, University of Bremen, Bibliothekstr. 1, Bremen, 28359, Germany</footnote>

<footnote><sup>b</sup> Centre Inria d'Université Côte d'Azur, 2004 Rte des Lucioles, Biot, 06410, France</footnote>

<footnote><sup>c</sup> Ioffe Institute, 26 Polyteknicheskaya str., St.Petersburg, 194021, Russia</footnote>

<footnote><sup>d</sup> Universite Claude Bernard Lyon 1, CNRS, Centrale Lyon, INSA Lyon, Université Jean Monnet, ICJ, UMR5208, Inria, Villeurbanne, 69622, France</footnote>

<footnote><sup>e</sup> Bernstein Center for Computational Neuroscience, Philippstr. 13, Haus 6, Berlin, 10115, Germany</footnote>

<footnote><sup>f</sup> Technische Universität Berlin Institute of Mathematics, Straße des 17. Juni 136, Berlin, 10623, Germany</footnote>

<footnote><sup>g</sup> Department of Mathematical Sciences "G. L. Lagrange", Politecnico di Torino, Corso Duca degli Abruzzi 24, Torino, 10129, Italy</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: pujo@math.univ-lyon1.fr (L. Pujo-Menjouet).</footnote>

<footnote>Peer review under the responsibility of KeAi Communications Co., Ltd.</footnote>
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 (Chizhov and Graham, 2007, 2008; Dumont et al., 2024; Dumont & Tarniceriu, 2024; Gerstner & Kistler, 2002; Pakdaman et al., 2009; Schwalger & Chizhov, 2019; Schwalger et al., 2017; Tarniceriu, 2020). 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 (Chizhov & Graham, 2007) or colored noise (Chizhov & Graham, 2008; Schwalger, 2021), while the classical RDE is based on a phenomenological hazard-rate-based spike-generation mechanism ("escape noise") (Gerstner & Kistler, 2002). Second, the RDE has been generalized to finite-size populations (Schmutz et al., 2023; Schwalger et al., 2017), 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 (Chizhov & Graham, 2007), and even two-compartmental (Chizhov & Graham, 2021) 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. (Li & Wang, 2002)), 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. (Della Marca et al., 2022; Loy & Tosin, 2021). 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.

## 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*. We suggest that the $V_i$-variable might somehow be related to the measurable concentrations difference between viral RNA and immune cells. Here, we refer to the units of this variable as “v-units”. 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 a time constant $\tau$. Second, it grows proportionally to the fraction of infectious 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^*)$.

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} < \dots$ 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^I$ (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 \le t^* < \tau_R$, respectively. However, the duration of these states does not directly

<table><caption>**Table 1**
Notation for the micro- and macroscopic models described in sections 2 and 3.</caption><thead><tr><th>Variable</th><th>Notation</th><th>Units</th></tr></thead><tbody><tr><td>t</td><td>Time</td><td>days</td></tr><tr><td>t*</td><td>Time since last infection, or “infection age”</td><td>days</td></tr><tr><td>V<sub>i</sub></td><td>Interplay of viral load and immune response for individual i</td><td>v-units</td></tr><tr><td>U</td><td>Mean interplay of viral load and immune response</td><td>v-units</td></tr><tr><td>ρ</td><td>Density of individuals distributed in t*-space</td><td>people/day</td></tr><tr><td>ν</td><td>Population rate of new infections</td><td>people/day</td></tr><tr><td>I</td><td>Fraction of infectious individuals in the population</td><td>nondim.</td></tr><tr><td>τ</td><td>Time scale of passive immune response</td><td>days</td></tr><tr><td>τ<sub>R</sub></td><td>Time scale of disease progress and cure</td><td>days</td></tr><tr><td>τ<sub>I</sub></td><td>Infectious period</td><td>days</td></tr><tr><td>τ<sub>R</sub></td><td>Time from infection to loss of immunity (“Refractory period”)</td><td>days</td></tr><tr><td>D</td><td>Time course of the interplay for disease and treatment</td><td>v-units/day</td></tr><tr><td>N</td><td>Total number of people in the population</td><td>people</td></tr><tr><td>N<sub>c</sub></td><td>Cumulative cases in the population</td><td>people</td></tr><tr><td>R<sub>0</sub></td><td>Basic Reproduction Number</td><td>nondim.</td></tr><tr><td>k</td><td>Average individual rate of viral load</td><td>v-units/day</td></tr><tr><td>H</td><td>Hazard; HΔt is the probability of infection in the interval (t, t + Δt)</td><td>people/day</td></tr><tr><td>V<sup>T</sup></td><td>Threshold of infectiousness</td><td>v-units</td></tr><tr><td>σ</td><td>Noise amplitude</td><td>v-units</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>**Fig. 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$, defined as the interplay between viral load and immune response, 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. Note that the SEIRS stages are loosely defined rather than strictly compartmentalized, and $V_i$ does not have a strictly defined zero level. **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.</figcaption></figure>

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, ..., N$, between two subsequent infection times $t_{k,i}$ and $t_{k+1,i}$ is modeled as
$$ 
\begin{aligned}
\frac{dV_i}{dt} &= -\frac{V_i}{\tau} + k(t)I(t) + D(t_i^*) + \sigma(I(t)) \xi_i(t), \\
\frac{dt_i^*}{dt} &= 1,
\end{aligned}
\quad (1) $$

with the fraction of infectious individuals

$$ I(t) = \frac{N_c(t) - N_c(t - \tau_l)}{N}, \quad (2) $$

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, \quad (3) $$

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

$$ D(t^*) = 4a\exp\left(-\frac{t^*}{\tau_r}\right) - a\exp\left(-\frac{t^*}{4\tau_r}\right) \quad (4) $$

is a response function that drives the time course of $V_i$ during the 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. (3) represents an independent Gaussian white noise with mean $E[\xi(t)] = 0$, auto-covariance function $E[\xi(t)\xi(t')] = \delta(t-t')$ 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) = 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 infectious, 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)). $$

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:

$$ \Pr\{\text{individual } i \text{ gets infected in } (t, t + \Delta t) \mid V_i(t), t_i^*(t^-)\} = H_i(t)\Delta t + o(\Delta t) \quad (5) $$

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

$$ H(V, t^*) = c \max(0, V/V^T)^m \mathbb{1}_{t^* \ge \tau_R}, \quad (6) $$

where $m > 0$. We illustrate the effect of varying $m$ in Fig. 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)) \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 the 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) = \bar{\sigma} \mathbb{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$,
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Escape noise in the form of (6). 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*₀(*t*), at *t* = 100. **A**, top to bottom: *I*₀(*t*), *k*(*t*), the fraction of infectious individuals *I*(*t*), the state variable *V*ᵢ(*t*) of a representative individual, and the rate of new cases *v*(*t*). For *I*(*t*) and *v*(*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 *ρ*H, and the density *ρ*, in 4 time moments. Parameters: *m* = 1, *c* = 0.015days⁻¹, *τ* = 50, *τ*ᵢ = 10, *τ*ᵣ = 1, *τ*<sub>R</sub> = 20, *a* = 100, and *V*<sup>*T*</sup> = 2.</figcaption></figure>

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}. \tag{7}
$$

Note that if this condition is fulfilled at some time $t_{k,i}^-$, the reset condition Eq. (3) 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
(Chizhov & Graham, 2007; Schwalger, 2021). In this case, the noise term in Eq. (1) is omitted, and the condition Eq. (7) is
substituted by (5), where $H_i(t)$ is now given by the following hazard function (8). The hazard rate $H(t, t^*)$ depends on the actual
state $V_i$, its rate of change $dV_i/dt$, the threshold $V^I$, and the noise amplitude $\sigma$, i.e., $H(t, t^*) = H(V_i(t), \hat{V}_i(t), t_i^*(t); \sigma(I(t)))$.
(Chizhov & Graham, 2007), this function was found as an approximate solution of a first-time passage problem based on the
Kolmogorov-Fokker-Planck equation:

$$
\begin{align}
H(U, \dot{U}, t^*, \sigma) &= \{A(T) + B(T, \dot{U}, \sigma), \text{ if } t^* > \tau_R, \sigma > 0; 0, \text{ otherwise}\}, \nonumber \\
A(T) &= \frac{1}{\tau} \exp(0.0061 - 1.12T - 0.257T^2 - 0.072T^3 - 0.0117T^4), \tag{8} \\
B(T, \dot{U}, \sigma) &= \frac{2}{\sigma\sqrt{\pi}} [\dot{U}]_{+} \frac{\exp(-T^2)}{1 + \operatorname{erf}(T)}, \quad T = \frac{V^T - U}{\sigma}. \nonumber
\end{align}
$$

Here and in the following, [x]_+ = max(0, x) denotes the rectified linear function. We refer to (Schwalger, 2021) 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 (Schwalger & Chizhov, 2019) 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. d/dt = ∂/∂t + ∂/∂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^*) dt^*.$$

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

$$v(t) = \lim_{\Delta t \to 0} \lim_{N \to \infty} \frac{N_c(t + \Delta t) - N_c(t)}{N \Delta t}.$$

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. (8), for which we can again define 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)–(3) consists of two transport equations for $\rho(t, t^*)$ and $U(t, t^*)$, and two integrals for $v(t)$ and $I(t)$:

$$\frac{\partial \rho}{\partial t} + \frac{\partial \rho}{\partial t^*} = -\rho H(U(t, t^*), t^*), \qquad (9a)$$

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

$$v(t) = \int_{0}^{\infty} \rho(t, t^*) H(U(t, t^*), t^*) dt^* + I_0 \delta(t), \quad \text{rate of new cases} \qquad (9c)$$

$$I(t) = \int_{t-\tau_l}^{t} v(s) ds, \quad \text{fraction of infectious people} \qquad (9d)$$

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

$$\rho(t, 0) = v(t). \qquad (10)$$

The boundary condition for $U$ corresponding to the reset in the microscopic model, Eq. (3), is $U(t, 0) = V_T$. The initial condition for the density $\rho$ corresponding to an initial fraction of infectious 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. (9b) (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 (9a), 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. (10). This boundary condition acts as a source term at $t^* = 0$ given by the rate of new cases $v(t)$ that exactly balances the integrated sink term, Eq. (9c). We note that this conservation law reflects the fact that we neglect mortality, and focus on the transmission mechanism.

Furthermore, in the density equation (9a), the risk of illness is evaluated by the hazard function $H$. In Eq. (9), 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. (8) and $\dot{U}(t, t^*) = (\partial_t + \partial_{t^*})U(t, t^*)$ is given by the right-hand-side of Eq. (9b).

Once again, note that Eq. (9d) permits a simple generalization to graded infectiousness depending on the course of the infection (see Discussion, Eq. (12)).
## 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 (Schmutz et al., 2023; Schwalger et al., 2017), 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^*) \nu_N(t - t^*). \qquad (11a) $$

Here, $S$ and $\nu_N$ are given as the solution of the following system of stochastic partial differential equations, generalizing Eq. (9)

$$ \begin{aligned} \frac{\partial S}{\partial t} + \frac{\partial S}{\partial t^*} &= -S H(U), & S(t, 0) &= 1, \\ \frac{\partial U}{\partial t} + \frac{\partial U}{\partial t^*} &= -\frac{U}{\tau} + u(t) + k I(t) + D(t^*), & U(t, 0) &= V_T, \end{aligned} \qquad (11b) \\ I(t) = \int_{t-\tau_l}^{t} \nu_N(s) \, ds, \qquad (11c) $$

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

$$ \nu_N(t) = \nu(t) + \sqrt{\frac{\nu(t)}{N}} \xi(t), $$

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^*}, \qquad (11d) $$

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]_+. $$

While in practice, the integral rarely 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. (11c).

Likewise, $\rho$ is not a normalized probability density (Schmutz et al., 2023; Schwalger et al., 2017). Therefore, it is no longer the strict interpretation of the empirical density of the times since the last infection but $\rho(t, t^*)dt^*$ that 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 \to \infty$, the solution converges to the solution of the macroscopic model. The precise numerical method employed for these simulations is described in (Schmutz et al., 2023) (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 the potentially infectious 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 susceptible again 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.
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Escape noise in the form of (6) with $m = 0.5, 1$ and $2$. The remaining parameter values are as in Fig. 2.</figcaption></figure>

The population waves move in the $t^*-\text{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 $v(t)$, the increase of the $\rho H$ determines the next peak of $v(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 $v(t)$ (Fig. 2A).

For the escape noise in the form of (6), 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 a steep $H$ with $m = 2$

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Epidemiology of Seasonal Coronaviruses. **A**, Data from (Nickbakhsh et al., 2020); specifically, the monthly prevalence of seasonal coronaviruses (sCoVs) detected among patients with respiratory illness viologically tested in NHS Greater Glasgow and Clyde, Scotland, United Kingdom. **B**, The case of oscillating $k$. Escape noise in the form of (6). The remaining parameter values are as in Fig. 2.</figcaption></figure>
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.

**Seasonality.** One of the experimentally observed features of epidemics is its seasonal oscillations, as illustrated in Fig. 4 by data from (Nickbakhsh et al., 2020). 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.

**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 $v(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 $v(t)$ and the fraction of the potentially infectious 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 $v(t)$ and $I(t)$. When $t = t_1 = 900$, the basic reproduction rate $k$ drops due to assumed anti-epidemic measures, and the epidemic wave generation stops.

**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. (11a)-(11d). 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 \to \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.

# 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 (Nickbakhsh et al., 2020). 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 (Jones et al., 2021), 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. Modeling infectiousness depending on infection age

It has been shown that the infectiousness of SARS-CoV-2 depends on the time since the last infection (Jones et al., 2021). This can be easily modeled in our framework by generalizing Eq. (2), Eq. (9d) and Eq. (11b) for the micro-, macro- and mesoscopic fractions of infectious people to the effective quantities

$$I(t) = \frac{1}{N} \int_{0}^{t} \chi(t-s) dN_{c}(s), \quad I(t) = \int_{0}^{t} \chi(t-s)v(s) ds \quad \text{and} \quad I(t) = \int_{0}^{t} \chi(t-s)v_{N}(s) ds, \qquad (12)$$

respectively. The kernel $\chi: \mathbb{R}^{+} \to [0, 1]$ describes the infectiousness depending on the infection age. Our equation (2), (9d) and Eq. (11b) with a fixed infection period $\tau_I$ are recovered as the special case $\chi(t^*) = \eta_{(0,\tau_I)}(t^*)$, allowing only two possible states – infectious and non-infectious. Our infectious period $\tau_I = 10$ days well correspond to the half-duration of cell culture infectivity measured in (Jones et al., 2021) for SARS-CoV-2.

We note that Eq. (12) may also effectively capture heterogeneous infectious periods $\tau_I$ in the mean-field limit if the infectious period $\tau_{I,i}$ of an individual $i$ is much shorter than the time interval between subsequent infections of an individual. Assuming that the infectious periods $\tau_{I,i}$ are distributed with probability density $\bar{p}(\tau_I)$, we have (in analogy to Eq. (9d))
<figure><img src="image_7.png" /><figcaption>**Fig. 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 *σ*(*t*) from 0.5 to 2, at *t* = *t*₀ = 0. A, top to bottom: *k*(*t*), *σ*(*t*), the fraction of potentially infectious population *I*(*t*), the state variable *V*ᵢ(*t*) of a representative individual, and the rate of new cases *v*(*t*). For *I*(*t*) and *v*(*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 *ρ*H, and the density *ρ*, in 4 time moments. C, Response to periodically changing *k* mimicking change of seasons (*V*<sup>*T*</sup> = 4). D, Response to periodically changing *k* mimicking change of seasons (*V*<sup>*T*</sup> = 3.7). Parameters: *τ* = 50, *τ*<sub>I</sub> = 10, *τ*<sub>R</sub> = 1, *τ*<sub>R</sub> = 20, *a* = 50, *V*<sup>*T*</sup> = 4, and *σ* = 2√π.</figcaption></figure>

$$
I(t) = \lim_{N \to \infty} \frac{1}{N} \sum_{i=1}^{N} \mathbb{I}_{\{t_i^*(t) < \tau_{Ii}\}} = \int_{0}^{\infty} \left( \int_{t-\tau_I}^{t} \nu(s|\tau_I) \, ds \right) \hat{p}(\tau_I) \, d\tau_I
$$

where $v(t|\tau_I)$ denotes the conditional rate of infections given the infectious period $\tau_I$ of an individual. If the typical time interval between subsequent infections is much longer than $\tau_I$, the conditional rate is approximately independent of $\tau_I$, i.e., $v(t|\tau_I) \approx v(t)$. Under this plausible assumption, we obtain

$$
I(t) \approx \int_{0}^{t} \left( \int_{t-s}^{\infty} \hat{p}(\tau_I) \, d\tau_I \right) \nu(s) \, ds.
$$

This expression is structurally identical to Eq. (12), provided that we define $\chi(t^*) = \int_{t^*}^{\infty} \hat{p}(\tau_I) \, d\tau_I$. Consequently, heterogeneity in $\tau_I$ can be effectively captured within our model using the generalized Eq. (12) to describe the fraction of infectious individuals $I(t)$.
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Mesoscopic model (6) for different $N = 200, 2000$, and $20000$, and different realizations of noise. The parameters are as in Fig. 2. The solutions with $N \to \infty$ converge to the macroscopic model solution shown in Fig. 2.</figcaption></figure>

### 6.3. Delay or anticipation of infectiousness relative to disease symptoms

In the present version, the onset of disease is the same as the onset of infectiousness. However, the model can be generalized to the case of arbitrary shift between these time moments. In the case of the disease symptoms following the infectiousness, the reaching of the threshold would define the onset of infectiousness, and the course of the disease $D(t*)$ would introduce a delaying zero-phase before the increasing phase. In the case of the disease symptoms preceding the infectiousness, the reaching of the threshold would define the onset of the disease. In both cases, the fraction of infectious people $I(t)$ would be defined with its extended version, Eq. (12).

### 6.4. 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 (Chizhov & Graham, 2021). 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. (Adimy et al., 2023; Ottaviano et al., 2023)). Since the multi-group approach naturally operates with smaller-size subpopulations, the proposed mesoscopic approach could be applied in this case.

### 6.5. 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, the activation and/or inactivation of some processes in the immune system. This concept resembles the consideration of Hodgkin-Huxley-like neurons in neuroscience (Chizhov & Graham, 2007). 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 (Loy & Tosin, 2021) 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 (Schwalger & Chizhov, 2019). 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.6. Age-dependence

The epidemics, to a different extent, affect people of different ages (Jones et al., 2021). While the current version of our model does not distinguish between age groups, it would certainly be interesting to include such a generalization in the future. Age could be included in the model as an additional independent variable; the density $\rho$ and the state variable $U$ would depend on $t$, $t^*$ and age $\hat{a}$. Age could be structured as a discrete (e.g., young, adult, old) parameter (this would be qualitatively equivalent to the multiple population scenario) or a continuous (proper age) variable. Specifically, 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 infectious individuals $I(t)$ would be an integral over $\hat{a}$. This parametrization would allow the consideration of different disease time courses, etc.

The infectiousness $\chi(t^*)$ in Eq. (12) as well as the rate $k(t)$ can also depend on age $\hat{a}$ because people usually interact more strongly within their age groups. In general, the interaction term $k(t)I(t)$ becomes a sum (or integral in the case of continuous $\hat{a}$) over different ages:

$$ \sum_{\hat{a}} k_{\hat{a}}(t) \int_{0}^{t} \chi_{\hat{a}}(t-s) dN_{c, \hat{a}}(s), \quad \sum_{\hat{a}} k_{\hat{a}}(t) \int_{0}^{t} \chi_{\hat{a}}(t-s) \nu_{\hat{a}}(s) ds, \quad \sum_{\hat{a}} k_{\hat{a}}(t) \int_{0}^{t} \chi_{\hat{a}}(t-s) \nu_{N, \hat{a}}(s) ds $$

for micro-, meso- and macroscopic cases, respectively.

## 6.7. 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 (Ansumali et al., 2020; Arino & Portet, 2020; Huo et al., 2021; Ottaviano et al., 2022, 2023). This may be done in multiple ways, depending on which characteristics of symptomatic and asymptomatic spread one wishes to capture in their model.

## 6.8. Heterogeneous severity

It could be of biological interest to assume a heterogeneous setting, e.g. the severity $a$ in Eq. (4) 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.9. 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^I$.

## 6.10. 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^I$ 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.11. 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.

### CRediT authorship contribution statement

**Anton Chizhov:** Writing – original draft, Validation, Methodology, Investigation, Conceptualization. **Laurent Pujo-Menjouet:** Writing – original draft, Validation, Supervision, Conceptualization. **Tilo Schwalger:** Writing – original draft, Validation, Methodology, Investigation, Formal analysis, Conceptualization. **Mattia Sensi:** Writing – original draft, Validation, Methodology, Investigation, Formal analysis, Conceptualization.

### Code availability

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

### Declaration of competing interest

The authors have no competing interest to declare.

### Acknowledgements

Mattia Sensi was partially supported by the Italian Ministry for University and Research (MUR) through the PRIN 2020 project "Integrated Mathematical Approaches to Socio-Epidemiological Dynamics" (No. 2020JLWP23, CUP: E15F21005420006).

## 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 $v(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 $v(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. (8) (Fig.A.7).

<figure><img src="image_9.png" /><figcaption>**Fig. A.7.** Solutions of a test with a stepwise term instead of $k(t)I(t)$ in Eqs. (1) and (9d), obtained with the microscopic model ($N = 2000$) with the white Gaussian noise and the escape noise described with Eq. (8), 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]$.</figcaption></figure>

From the methodical point of view, Fig. A.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.

### References

* Adimy, M., Chekroun, A., Pujo-Menjouet, L., & Sensi, M. (2023). A multigroup approach to delayed prion production. *Discrete and continuous dynamical systems-B*.

* Alfaro-Murillo, J. A., Feng, Z., & Glasser, J. W. (2019). Analysis of an epidemiological model structured by time-since-last-infection. *Journal of Differential Equations*, *267*(10), 5631–5661.

* Anderson, R. M. (1988). The epidemiology of HIV infection: Variable incubation plus infectious periods and heterogeneity in sexual activity. *Journal of the Royal Statistical Society: Series A*, *151*(1), 66–93.

* Ansumali, S., Kaushal, S., Kumar, A., Prakash, M. K., & Vidyasagar, M. (2020). Modelling a pandemic with asymptomatic patients, impact of lockdown and herd immunity, with applications to SARS-CoV-2. *Annual Reviews in Control*, *50*, 432–447.
Arino, J., & Portet, S. (2020). A simple model for COVID-19. *Infectious Disease Modelling, **5***, 309–315.

Brauer, F., & Watmough, J. (2009). Age of infection epidemic models with heterogeneous mixing. *Journal of Biological Dynamics, **3*** (2–3), 324–330.

Chizhov, A. V., & Graham, L. J. (2007). Population model of hippocampal pyramidal neurons, linking a refractory density approach to conductance-based neurons. *Physical Review E, **75*** (1), Article 011924.

Chizhov, A. V., & Graham, L. J. (2008). Efficient evaluation of neuron populations receiving colored-noise current based on a refractory density method. *Physical Review E, **77*** (1), Article 011910.

Chizhov, A. V., & Graham, L. J. (2021). A strategy for mapping biophysical to abstract neuronal network models applied to primary visual cortex. *PLoS Computational Biology, **17*** (8), Article e1009007.

Della Marca, R., Loy, N., & Tosin, A. (2022). An SIR–like kinetic model tracking individuals’ viral load. *Networks and Heterogeneous Media, **17*** (3), 467–494.

d’Onofrio, A., Iannelli, M., Manfredi, P., & Marinoschi, G. (2023). Optimal epidemic control by social distancing and vaccination of an infection structured by time since infection: The COVID-19 case study. *SIAM Journal on Applied Mathematics, **S199***–S224.

Dumont, G., Henry, J., & Tarniceriu, C. O. (2024). Oscillations in a fully connected network of leaky integrate-and-fire neurons with a Poisson spiking mechanism. *Journal of Nonlinear Science, **34*** (1), 18.

Dumont, G., & Tarniceriu, C. O. (2024). Pattern Formation in a spiking neural-field of renewal neurons. *SIAM Journal on Applied Dynamical Systems, **23*** (4), 2670–2694.

Gerstner, W., & Kistler, W. M. (2002). *Spiking neuron models: Single neurons, populations, plasticity*. Cambridge University Press. Google-Books-ID: Rs4oc7HfxIUC.

Huo, X., Chen, J., & Ruan, S. (2021). Estimating asymptomatic, undetected and total cases for the COVID-19 outbreak in wuhan: A mathematical modeling study. *BMC Infectious Diseases, **21*** (1), 476.

Iannelli, M., & Milner, F. (2017). The basic approach to age-structured population dynamics. *Lecture notes on mathematical modelling in the life sciences*. Dordrecht: Springer.

Jones, T. C., Biele, G., Mühlemann, B., Veith, T., Schneider, J., Beheim-Schwarzbach, J., et al. (2021). Estimating infectiousness throughout SARS-CoV-2 infection course. *Science, **373*** (6551), Article eabi5273.

Kermack, W. O., & McKendrick, A. G. (1927). A contribution to the mathematical theory of epidemics. *Proceedings of the Royal Society of London - Series A: Containing Papers of a Mathematical and Physical Character, **115*** (772), 700–721.

Li, L. L., Ferreira, C. P., & Ainseba, B. (2019). Mathematical analysis of an age structured problem modelling phenotypic plasticity in mosquito behaviour. *Nonlinear Analysis: Real World Applications, **48***, 410–423.

Li, M. Y., & Wang, L. (2002). Global stability in some SEIR epidemic models. In *Mathematical approaches for emerging and reemerging infectious diseases: Models, methods, and theory* (pp. 295–311). Springer.

Loy, N., & Tosin, A. (2021). Boltzmann-type equations for multi-agent systems with label switching. *Kinetic and Related Models, **14*** (5), 867–894.

Nickbakhsh, S., Ho, A., Marques, D. F. P., McMenamin, J., Gunson, R. N., & Murcia, P. R. (2020). Epidemiology of seasonal coronaviruses: Establishing the context for the emergence of coronavirus disease 2019. *The Journal of Infectious Diseases, **222*** (1), 17–25.

Ottaviano, S., Sensi, M., & Sottile, S. (2022). Global stability of SAIRS epidemic models. *Nonlinear Analysis: Real World Applications, **65***, Article 103501.

Ottaviano, S., Sensi, M., & Sottile, S. (2023). Global stability of multi-group SAIRS epidemic models. *Mathematical Methods in the Applied Sciences, **46*** (13), 14045–14071.

Pakdaman, K., Perthame, B., & Salort, D. (2009). Dynamics of a structured neuron population. *Nonlinearity, **23*** (1), 55–75.

Schmutz, V., Löcherbach, E., & Schwalger, T. (2023). On a finite-size neuronal population equation. *SIAM Journal on Applied Dynamical Systems, **22*** (2), 996–1029.

Schwalger, T. (2021). Mapping input noise to escape noise in integrate-and-fire neurons: A level-crossing approach. *Biological Cybernetics, **115***, 539–562.

Schwalger, T., & Chizhov, A. V. (2019). Mind the last spike—firing rate models for mesoscopic populations of spiking neurons. *Current Opinion in Neurobiology, **58***, 155–166.

Schwalger, T., Deger, M., & Gerstner, W. (2017). Towards a theory of cortical columns: From spiking neurons to interacting neural populations of finite size. *PLoS Computational Biology, **13*** (4), Article e1005507.

Tarniceriu, C. O. (2020). Age-structure in neuronal models. *Annals of the Alexandru Ioan Cuza University - Mathematics, **66*** (2), 385–396.