# A Kermack-McKendrick model with age of infection starting from a single or
multiple cohorts of infected patients

Jacques Demongeot(a), Quentin Griette(b),111Q.G. acknowledges support from ANR
via the project Indyana under grant agreement ANR-21-CE40-0008. , Yvon
Maday(c), and Pierre Magal(b),
(a)Université Grenoble Alpes, AGEIS EA7407, F-38700 La Tronche, France
(b)Univ. Bordeaux, IMB, UMR 5251, F-33400 Talence, France.
CNRS, IMB, UMR 5251, F-33400 Talence, France.
(c)Sorbonne Université and Université de Paris, CNRS, Laboratoire Jacques-
Louis Lions (LJLL), F-75005 Paris, France
Institut Universitaire de France, 75005 Paris, France Corresponding author.
e-mail<EMAIL_ADDRESS>
###### Abstract

During an epidemic, the infectiousness of infected individuals is known to
depend on the time since the individual was infected, that is called the age
of infection. Here we study the parameter identifiability of the Kermack-
McKendrick model with age of infection which takes into account this
dependency. By considering a single cohort of individuals, we show that the
daily reproduction number can be obtained by solving a Volterra integral
equation that depends on the flow of new infected individuals. We test the
consistency our the method by generating data from deterministic and
stochastic numerical simulations. Finally we apply our method to a dataset
from SARS-CoV-1 with detailed information on a single cluster of patients. We
stress the necessity of taking into account the initial data in the analysis
to ensure the identifiability of the problem.

Keywords: Age of infection; epidemic model; single and multiple cohorts;
cluster of infected patients; daily reproduction number; Volterra equation;
parameter identification.

##  1 Introduction

The existence of an individual infection (or contagiousness) period of
variable length among infected individuals and during the history of an
epidemic is a proven fact in all contagious diseases. The origin of this
variability is multiple. It may be due to

  * (i)

a variation in the symptomatic state of the infected person (such as the
change in frequency and intensity of the cough, source of emission of the
infectious agent) due to variable immune defenses since the beginning of his
infection (e.g., due to the transition between initial innate immune reaction
and secondary adaptive reaction);

  * (ii)

a variation in the environmental conditions of transport and survival of the
infectious agent in the atmosphere (heat, humidity, sunshine, altitude, etc.),
in a more or less favorable socio-sanitary environment presenting different
spreading characteristics (as compliance with distancing or confinement,
existence of preventive vitamin supplementation, educational level, food
quality, demographic isolate, etc.);

  * (iii)

a variation in the state of defense of the final host (immune defense and
self-protection);

  * (iv)

a variation in the virulence of the infectious agent, which can mutate or
possibly change of intermediary host.

Due to the complexity of the contagion mechanisms and of their variations,
accentuated by the fact that the rate of the epidemic may be the result of
several geographically distant clusters starting at different times, it is
preferable, in a first approach, to assume that the initial conditions are
reduced to a single cluster. The challenge is then to estimate, on each day of
the infection period, the transmissibility rate allowing the calculation of
the daily reproduction number of the contagious disease.

Continuous time model: Recall that the age of infection $a$ is the time since
individuals become infected. The major difficulty in comparing the data and
the Kermack-McKendrick model with age of infection is to identify: 1) the
initial distribution of infected with respect to the age of infection; 2) the
daily reproduction number $R_{0}(a)$ which is the reproduction number at the
age of infection $a$ (i.e. the average number of secondary cases produced by a
single infected at the age of infection $a$). We can decompose the daily
reproduction number as follows

$R_{0}(a)=\underset{{\rm(A)}}{\underbrace{\tau_{0}}}\times\underset{{\rm(B)}}{\underbrace{S_{0}}}\times\underset{{\rm(C)}}{\underbrace{\beta(a)}}\times\underset{{\rm(D)}}{\underbrace{e^{-\nu
a}}},$

where (A) is the rate of transmission rate at time $t_{0}$ (we assume the
transmission rate to be constant during the period where $R_{0}(a)$ is
evaluated), (B) is the number of susceptible individuals at time $t_{0}$ (we
assume the number of susceptibles to be constant during the period where
$R_{0}(a)$ is evaluated), (C) is the probability that an infected who have
been for $a$ days is infectious i.e., capable to transmit the pathogen (or the
fraction of infectious among the infected with age of infection $a$), (D) the
probability for an infected for $a$ days to be still infected.

Then the basic reproduction number (i.e. the number of secondary cases
produced by a single infected) is given by

$R_{0}=\int_{0}^{\infty}R_{0}(a)da.$

Here we partly solve the problem of finding the initial distribution of
infected by assuming that we start the epidemic at time $t_{0}$ with a single
cohort of $I_{0}$ new infected patients. That is, the epidemic starts with
$I_{0}$ infected patients all with age of infection $a=0$. The case of an
epidemic starting from a single infected patient (usually called the patient
$0$) corresponds to the case $I_{0}=1$. This is a common assumption in
epidemiology. Note that the time $t_{0}$ at which the first patient becomes
infected is also unknown for most epidemics.

Assume that the epidemic starts at time $t_{0}$ with a cohort of $I_{0}$ new
infected patients (i.e., all with age of infection $a=0$). Then $N(t)$ the
flow of infected at time $t$ satisfies the model starting from a single cohort
of infected

$N(t)=\underset{{\rm(I)}}{\underbrace{R_{0}(t-t_{0})\times
I_{0}}}+\underset{{\rm(II)}}{\underbrace{\int_{0}^{t-t_{0}}R_{0}(s)\times
N(t-s)ds}},\forall t\geq t_{0},$ (1.1)

where (I) is the flow of infected produced directly by the $I_{0}$ infected
individuals already present on day $t_{0}$; and (II) is the flow of new
infected individuals at time $t$ produced by the new infected individuals
since day $t_{0}$.

Here what we call “flow of new infected individuals” is the distribution
$N(t)$ such that

$\int_{t_{1}}^{t_{2}}N(\sigma)d\sigma,$

is the number of new infected individuals during the period of time
$[t_{1},t_{2}]$.

The equation (1.1) will be derived as an extension of the Kermac-McKendrick
model starting from a single cohort of infected (see equation (4.5)). This
model remains valid as long as the transmission rate $\tau(t)$ and the number
of susceptible hosts $S(t)$ remain constant. So this model is valid when the
epidemic starts.

Assume that $I_{0}$ is fixed and the function $a\to R_{0}(a)$ is given. Then
the map $t\to N(t)$ can be obtained by solving (1.1). The goal of the article
is consider the converse problem. That is, assume that $I_{0}$ is fixed and
assume that $t\to N(t)$ is given from the data. Then the map $a\to R_{0}(a)$
can be obtained by solving the Volterra integral equation

$R_{0}(a)=\dfrac{N(a+t_{0})}{I_{0}}-\dfrac{1}{I_{0}}\int_{0}^{a}R_{0}(s)N(a-s+t_{0})ds,\forall
a\geq 0.$ (1.2)

Therefore if the map $t\to N(t)$ is known, we can theoretically derive the
average dynamics of infection at the level of a single patient.

The standard assumption used in the literature is

$R_{0}(a)=0,\forall a\geq a_{+},$ (1.3)

where $a_{+}>0$ the maximal age of infectiousness for an infected patient.

Then for $t\geq t_{0}+a_{+},$ the equation (1.1) becomes

$N(t)=\int_{0}^{a_{+}}R_{0}(s)\times N(t-s)ds.$ (1.4)

Now assume for example that $N(t)=N_{0}e^{\lambda t},$ we obtain after
simplifications a standard characteristic equation

$1=\int_{0}^{a_{+}}R_{0}(s)\times e^{-\lambda s}ds.$ (1.5)

Assume that $\lambda>0$ is given. Consider $a\to\chi(a)$ any non-negative and
non-null continuous function satisfying

$\chi(a)=0,\forall a\geq a_{+}.$

Then

$R_{0}(a)=\dfrac{\chi(a)}{\int_{0}^{a_{+}}\chi(s)\times e^{-\lambda s}ds}$

satisfies (1.5). Therefore, neglecting the initial value in the Volterra
equation leads to a non-identifiable problem (in general). This shows the
crucial role of the initial value in identifying the function $a\to R_{0}(a)$.

Day by day model: The model (1.1) with a single cohort of infected becomes a
discrete Volterra equation

$N(t)=\underset{{\rm(I)}}{\underbrace{R_{0}(t-t_{0})\times
I_{0}}}+\underset{{\rm(II)}}{\underbrace{\sum_{d=1}^{t-t_{0}}R_{0}(d)\times
N(t-d)}},\forall t\geq t_{0},$ (1.6)

where (I) is the number of infected produced directly by the $I_{0}$ infected
individuals already present on day $t_{0}$; and (II) is the number of new
infected individuals at time $t$ produced by the new infected individuals
since day $t_{0}$.

Next, by setting $a=t-t_{0}$, we obtain the day by day equation for the daily
reproduction number

$R_{0}(a)=\dfrac{N(t_{0}+a)}{I_{0}}-\dfrac{1}{I_{0}}\sum_{d=1}^{a}R_{0}(d)\times
N(t_{0}+a-d),\forall a\geq 0.$ (1.7)

In the above formula and throughout the paper, we use the following convention
for the sum

$\displaystyle\sum_{d=k}^{m}=0,\text{ whenever }m<k.$

If we consider the first terms of the discrete time Volterra equation (1.6),
we obtain

$\begin{array}[]{ll}N(t_{0})&=R_{0}(0)\times I_{0},\vspace{0.1cm}\\\
N(t_{0}+1)&=R_{0}(1)\times I_{0}+R_{0}(1)\times N(t_{0}),\vspace{0.1cm}\\\
N(t_{0}+2)&=R_{0}(2)\times I_{0}+R_{0}(2)\times N(t_{0})+R_{0}(1)\times
N(t_{0}+1),\vspace{0.1cm}\\\ N(t_{0}+3)&=R_{0}(3)\times I_{0}+R_{0}(3)\times
N(t_{0})+R_{0}(2)\times N(t_{0}+1)+R_{0}(1)\times N_{0}(t_{0}+2),\\\
\quad\vdots\end{array}$

In practice, we can assume that $R_{0}(0)=0$ since infected individuals are
not infectious immediately after being infected. Under this additional
assumption, we obtain the system

$\begin{array}[]{ll}N(t_{0})&=0,\vspace{0.1cm}\\\ N(t_{0}+1)&=R_{0}(1)\times
I_{0},\vspace{0.1cm}\\\ N(t_{0}+2)&=R_{0}(2)\times I_{0}+R_{0}(1)\times
N(t_{0}+1),\vspace{0.1cm}\\\ N(t_{0}+3)&=R_{0}(3)\times I_{0}+R_{0}(2)\times
N(t_{0}+1)+R_{0}(1)\times N_{0}(t_{0}+2),\\\ \quad\vdots\end{array}$

When reliable information is available on the first cluster(s), the best
formula for calculating daily basic reproduction numbers is equation (1.2) (or
its the discrete time version (1.7)). Based on (1.4), some methods have been
developed in the literature to cope with the lack of a precise information.

For instance in [2, 3], the authors following [27] propose an optimization
algorithm for estimating the daily basic reproduction numbers. In [4], D.
Bernoulli mentions in 1760 the changes in the contagiousness parameters and
places as a crucial challenge for the prediction of the transition between
endemic and epidemic peaks in a prophetic sentence: “Le retour d’une épidémie
longtemps suspendue fait un ravage plus terrible dans une seule année qu’une
endémie uniforme ne pourrait faire pendant un nombre d’années considérable”
(The return of a long-suspended epidemic wreaks more terrible havoc in a
single year than a uniform endemic could do for a considerable number of
years). In [10], the authors use a deconvolution algorithm for calculating the
daily basic reproduction numbers. In each case, the problem of the initial
conditions is evoked at best only through the hypothesis of a unique “patient
zero”. Despite the considerable means of current investigation, in particular
those of the WHO and the members of the government of the WHO, it is rare that
this patient is identified (this was the case for H1N1 in Mexico). The patient
zero, also called index or primary case is the first patient identified in a
given population during an epidemiological investigation. It points out the
source of the spread of a disease in a given reservoir, but this search is in
general very difficult as it was the case for HIV in North America [32].

One of the main difficulties in estimating the $R_{0}(a)$ function is its non-
identifiability in general. Recent studies [25, 8, 16, 17, 9] developed
methods to identify the various parameters for the COVID-19 pandemic by using
cumulative reported cases data and differential equations models. Differential
equations can be written in the form studied here by assuming that $R_{0}(a)$
(or equivalently, $\beta(a)$) is independent of $a$. Suppose that we are
restricted to a period when the data is growing exponentially fast. If we take
a fixed function $\beta(a)$, then by adapting the method developed in [25], we
could identify a transmission rate $\tau$ so that the output of the model
stays very close to the data, for any function $\beta(a)$. The same could be
achieved with a good phenomenological description of the data by using the
method developed in [8, 16, 17] with a time-dependent transmission rate. This
means that the reported cases data is not sufficient to determine accurately
the function $R_{0}(a)$. Without a good description of the initial
distribution, it is hopeless to identify $R_{0}(a)$ by using reported cases
only.

In this article, we first extend the Kermack-McKendrick model to initial
conditions that are a linear combination of Dirac masses. The PDE model (2.1)
can not be extended in the space of measures (due to a lack of time continuity
for the solutions). However, the Volterra integral equation can be extended
and still makes perfect sense whenever we use Dirac masses for the initial
condition. In many real examples, the initial distribution must be a linear
combination of Dirac masses since the data are discrete at the early stage of
an epidemic (in a city, a country). Indeed, at the early beginning of an
epidemic, the epidemic starts from a few cases imported from other places.
Therefore Dirac masses make perfect sense. In practice, the early stage of an
outbreak is often undocumented and generally difficult to determine. But our
study applies to data from a finite collection of clusters that is easier to
determine using contact tracing.

Consequently, for the single cohort model, we can reverse the problem, and by
assuming that the daily number of new infected is known, we can compute the
daily reproduction number by solving a Volterra integral equation. The daily
basic reproduction number informs us about the dynamics of infection at the
level of a single patient. Therefore, knowing $R_{0}(a)$ should help the
medical doctors decide about quarantine measures. Reported case data for
clusters are particularly valuable for reconstructing the dynamics of
infection at the level of a single individual.

In this paper, we also provide an Individual-Based Model (IBM) (see Appendix
B). This IBM converges to the deterministic model whenever the initial number
of infected increases. We use this IBM to generate sample data to test our
method and compute the daily basic reproduction number. This will allow us to
test the effects of the day-to-day discretization (on the data) and the impact
of stochastic perturbations on the daily reproduction number. We conclude the
paper by applying our approach to a cluster of SARS-CoV-1 in Singapore.

The plan of the paper is the following. In Section 2, we recall the Kermack-
McKendrick model with age of infection. We explain how to derive the Volterra
formulation of the model, and we compare it with the Kermack-McKendrick SI
model with age of infection (ODE model). In Section 3, we explain how to
connect the model with the data. In Section 4, we extend the Kermack-
McKendrick model with age of infection in the case where the epidemic starts
from a single or multiple cohorts of infected individuals. In Section 5 we
derive an equation to compute the daily reproduction number from the data. In
Section 6 we consider a day by day discretized Kermack-McKendrick model with
age of infection. In Section 7, we run some numerical simulations, we compare
the deterministic model with a stochastic individual based simulation
presented in appendix. In Section 8, we compare the model with some data from
SARS-CoV-1, and we discuss the data from SARS-CoV-2.

##  2 Kermack-McKendrick model with age of infection

###  2.1 Partial differential equation formulation of the model

The age of infection $a$ is the time since individuals become infected. Let
$a\to i(t,a)$ be the distribution of population of infected individuals at
time $t$ (with respect to $a$ the age of infection). The term distribution of
population means that the integral

$\int_{a_{1}}^{a_{2}}i(t,a)da$

is the number of infected at time $t$ with infection age between $a_{1}$ and
$a_{2}$. Therefore the total number of infected individuals at time $t$ is

$I(t)=\int_{0}^{+\infty}i(t,a)da.$

Let $\beta(a)\in[0,1]$ be the probability to be contagious or infectious (i.e.
capable to transmit the pathogen) at the age of infection $a$. The quantity
$\beta(a)$ can be interpreted as the fraction of infected individuals with age
of infection $a$ that are infectious. Then the total number of contagious
individuals (or also called infectious individuals) (i.e., the individuals
capable of transmitting the pathogen) at time $t$ is

$C(t)=\int_{0}^{+\infty}\beta(a)i(t,a)da.$

The model of Kermack-McKendrick [22] with age of infection is the following,
for each $t\geq t_{0}$

$\left\\{\begin{array}[]{l}S^{\prime}(t)=-\tau(t)\,S(t)\,\displaystyle\int_{0}^{+\infty}\beta(a)\,i(t,a)da,\vspace{0.2cm}\\\
\partial_{t}i+\partial_{a}i=-\nu\,i(t,a),\text{ for }a\geq 0,\vspace{0.2cm}\\\
i(t,0)=\tau(t)\,S(t)\,\displaystyle\int_{0}^{+\infty}\beta(a)\,i(t,a)da,\end{array}\right.$
(2.1)

this system is supplemented by initial data

$S(t_{0})=S_{0}\geq 0,\text{ and }i(t_{0},a)=i_{0}(a)\in L^{1}_{+}(0,\infty).$
(2.2)

In the model, $S(t)$ is the number of susceptible individuals at time $t$, and
$t\to\tau(t)$ is the transmission rate at time $t$, and $\nu\geq 0$ is the
rate at which individuals die or recover. Here, the parameter $\nu$ is assumed
to be independent of the age of infection $a$. This is a simplifying
assumption to improve the readability of the paper. The parameter $\nu$
combines both the specific fatality rate and the recovery rate.

The above equation can be understood first as follows

$I^{\prime}(t)=\underset{{\rm(I)}}{\underbrace{\tau(t)\,S(t)\,\int_{0}^{+\infty}\beta(a)\,i(t,a)da}}-\underset{{\rm(II)}}{\underbrace{\int_{0}^{+\infty}\nu\,i(t,a)da}},$

where (I) is the flow of new infected, and (II) is the flow of individuals who
die or recover.

We make the following assumption.

######  Assumption 2.1.

We assume that

  * (i)

The transmission rate $t\to\tau(t)$ is a bounded continuous map from
$[t_{0},+\infty)$ in $[0,+\infty)$;

  * (ii)

The probability to be infectious at the age of infection $a\to\beta(a)\in
L^{\infty}_{+}(0,+\infty)$ is a non-negative and measurable function of $a$
which is bounded by $1$;

###  2.2 Volterra integral equation formulation of the model

In the model (2.1), the quantity

$N(t):=\tau(t)\,S(t)\,\int_{0}^{+\infty}\beta(a)\,i(t,a)da,$ (2.3)

is the flow of new infected individuals at time $t$.

By using the $S$-equation in system (2.1), we obtain

$S(t)=S_{0}-\int_{t_{0}}^{t}N(\sigma)d\sigma,\forall t\geq t_{0}.$ (2.4)

By integrating the $i$-equation of system (2.1) along the characteristics, we
obtain

$i(t,a)=\left\\{\begin{array}[]{ll}e^{-\nu\,\left(t-t_{0}\right)}\,i_{0}\left(a-\left(t-t_{0}\right)\right),&\text{
if }a\geq t-t_{0},\vspace{0.2cm}\\\ e^{-\nu\,a}N(t-a),&\text{ if }t-t_{0}\geq
a.\\\ \end{array}\right.$ (2.5)

By using (2.5), we deduce that $t\to N(t)$ satisfies the following Volterra
integral equation

$N(t)=\underset{{\rm(I)}}{\underbrace{\tau(t)\,S(t)\,\int_{t-t_{0}}^{+\infty}\beta(a)\,e^{-\nu\,\left(t-t_{0}\right)}i_{0}\left(a-\left(t-t_{0}\right)\right)da}}+\underset{{\rm(II)}}{\underbrace{\tau(t)\,S(t)\,\int_{0}^{t-t_{0}}\beta(a)e^{-\nu\,a}N(t-a)\,da}}$
(2.6)

where ${\rm(I)}$ is the flow of new infected individuals at time $t$ produced
by the infected individuals already present on day $t_{0}$; ${\rm(II)}$ is the
flow of new infected individuals at time $t$ produced by the new infected
individuals since day $t_{0}$.

By using equations (2.4) and (2.6), we can summarize the epidemic model (2.1),
by saying that $t\to N(t)$ is the unique continuous map satisfying

$N(t)=\tau(t)\,S(t)\left[\Lambda(t)+\int_{0}^{t-t_{0}}\beta(a)e^{-\nu\,a}N(t-a)\,da\right],\forall
t\geq t_{0},$ (2.7)

where

$S(t)=S_{0}-\int_{t_{0}}^{t}N(\sigma)d\sigma,\forall t\geq t_{0},$ (2.8)

and

$\Lambda(t):=e^{-\nu\,\left(t-t_{0}\right)}\int_{t-t_{0}}^{+\infty}\beta(a)\,i_{0}\left(a-\left(t-t_{0}\right)\right)da,\forall
t\geq t_{0}.$ (2.9)

The function $\Lambda(t)$ is the number of infectious individuals (capable to
transmit the pathogen) at time $t$ among the infected individuals already
present at time $t_{0}$.

The function $t\to\Lambda(t)$ plays a fundamental role in solving the Volterra
equation. Indeed, the quantity

$\int_{t_{1}}^{t_{2}}\tau(\sigma)\,S(\sigma)\Lambda(\sigma)d\sigma,$

is the number of infected produced between the instants $t_{1}$ and $t_{2}$ by
the infected already present at time $t_{0}$. So, for example, if no new
infected are produced by the infected already present at time $t_{0}$, that is
if

$\Lambda(t)=0,\forall t\geq t_{0},$

then there will be no new infected at all after the time $t_{0}$, that is

$N(t)=0,\forall t\geq t_{0}.$

The function $t\to\Lambda(t)$ can be regarded as the initial distribution of
the Volterra integral equation (2.7).

##  3 Connecting the data and the model

The data are represented by the function $t\to{\rm CR}(t)$ which is the
cumulative number of reported cases at time $t$. We assume that the flow of
reported cases is a fraction $0\leq f\leq 1$ of the flow of recovering
individuals, that is

${\rm CR}^{\prime}(t)=f\,\nu\,\int_{0}^{+\infty}i(t,a)da.$ (3.1)

By using (2.5), we can compute the number of infected at time $t$. That is

$\int_{0}^{+\infty}i(t,a)da=e^{-\nu\,\left(t-t_{0}\right)}I_{0}+\int_{0}^{t-t_{0}}e^{-\nu\,a}N(t-a)\,da,$
(3.2)

where

$I_{0}=\int_{0}^{+\infty}i_{0}(a)da$

is the total number of infected at time $t_{0}$.

By using equations (3.1) and (3.2), we obtain

${\rm CR}^{\prime}(t)=\nu
f\left[e^{-\nu\,\left(t-t_{0}\right)}I_{0}+\int_{0}^{t-t_{0}}e^{-\nu\,a}N(t-a)da\right],$

or equivalently (by using the change of variable $\sigma=t-a$)

${\rm CR}^{\prime}(t)=\nu
f\left[e^{-\nu\,\left(t-t_{0}\right)}I_{0}+\int_{t_{0}}^{t}e^{-\nu\,\left(t-\sigma\right)}N(\sigma)d\sigma\right].$

By choosing $t=t_{0}$ we obtain

$I_{0}=\dfrac{{\rm CR}^{\prime}(t_{0})}{\nu f},$

and

$\int_{t_{0}}^{t}e^{\nu\sigma}N(\sigma)d\sigma=\dfrac{e^{\nu\,t}{\rm
CR}^{\prime}(t)}{\nu f}-e^{\nu\,t_{0}}I_{0},$

and by differentiating both sides of the above equation, we obtain

$e^{\nu t}N(t)=\dfrac{\nu e^{\nu\,t}{\rm CR}^{\prime}(t)+e^{\nu\,t}{\rm
CR}^{\prime\prime}(t)}{\nu f}.$

Therefore we obtain the following connection between the data and the model.

Connection between the data and the model Let $t\to{\rm CR}(t)$ be the
cumulative number of reported cases. Then the initial number of infected is
given by $I_{0}=\dfrac{{\rm CR}^{\prime}(t_{0})}{\nu f},$ (3.3) and the flow
of new infected individuals $N(t)$ at time $t$ is given by
$N(t)=\dfrac{\nu{\rm CR}^{\prime}(t)+{\rm CR}^{\prime\prime}(t)}{\nu
f},\forall t\geq t_{0}.$ (3.4)

##  4 Kermack-McKendrick model starting from a single and multiple cohorts of
infected patients

The major difficulty to compare the model (2.4) with the data is to identify
the functions $a\to i_{0}(a)$ and $a\to\beta(a)$. To simplify the discussion,
let us consider the model at the early stage of the epidemic. When the
epidemic just starts we can assume that the transmission rate $t\to\tau(t)$
remains constant, and the number of susceptible individuals $t\to S(t)$ is
constant and equal to $S_{0}$. Under such a simplifying assumption the
Volterra equation (2.4) becomes

$N(t)=\tau\,S_{0}\left[\Lambda(t)+\int_{0}^{t-t_{0}}\beta(a)e^{-\nu\,a}N(t-a)\,da\right],\forall
t\geq t_{0}.$ (4.1)

###  4.1 Initial distribution for a single cohort of infected with age of
infection $a=0$

In order to understand the mathematical concept of Dirac mass centered at $0$,
we first consider an approximation by an exponential law

$i_{0}(a)=I_{0}\,\kappa\,e^{-\kappa a},$

mean and standard deviation equal to $1/\kappa$. Then a Dirac mass centered at
$0$ can be understood as the limit of such a distribution when $\kappa$ goes
to $+\infty$. The limit of needs some explanations. Recall that

$\int_{a_{1}}^{a_{2}}i_{0}(a)da=I_{0}\,\left[e^{-\kappa a_{1}}-e^{-\kappa
a_{2}}\right],$

is the initial number of infected individuals with infection age $a$ in
between $a_{1}$ and $a_{2}$ at time $t=0$.

We deduce that

$\lim_{\kappa\to\infty}\int_{a_{1}}^{a_{2}}i_{0}(a)da=\left\\{\begin{array}[]{lr}0,&\text{
if }a_{2}>a_{1}>0,\\\ I_{0},&\text{ if }a_{2}>a_{1}=0.\\\ \end{array}\right.$

That is to say that, when $\kappa$ tends to $+\infty$, the initial
distribution of population $i_{0}(a)$ is approaching the case where all the
infected individuals at time $t_{0}$ have the same age of infection $a=0$.

For short, we write

$i_{0}(a)=I_{0}\,\delta_{0}(a),$

where $\delta_{0}(a)$ is called the Dirac mass centered at $0$.

###  4.2 Model starting from a single cohort of infected with age of infection
$a=0$

Recall that

$\Lambda(t)=e^{-\nu\,\left(t-t_{0}\right)}\int_{0}^{+\infty}\beta\left(a+\left(t-t_{0}\right)\right)\,i_{0}\left(a\right)da,$

so when $i_{0}(a)$ is replaced by $I_{0}\times\kappa\times e^{-\kappa a}$ we
obtain

$f_{\kappa}(t):=I_{0}e^{-\nu\,\left(t-t_{0}\right)}\int_{0}^{+\infty}\beta\left(a+\left(t-t_{0}\right)\right)\times\kappa\times
e^{-\kappa a}da.$

In order to derive the Kermack-McKendrick model with Dirac mass initial
distribution as limit, we first need the following result.

######  Lemma 4.1.

Let Assumption 2.1 be satisfied, and assume in addition that $a\to\beta(a)$ is
continuous. Then we have

$\lim_{\kappa\to\infty}f_{\kappa}(t)=I_{0}e^{-\nu\,\left(t-t_{0}\right)}\beta\left(t-t_{0}\right),$

where the limit is uniform in $t\geq t_{0}$. That is

$\lim_{\kappa\to+\infty}\sup_{t\geq
t_{0}}|f_{\kappa}(t)-I_{0}e^{-\nu\,\left(t-t_{0}\right)}\beta\left(t-t_{0}\right)|=0.$

###### Proof.

Let $\varepsilon>0$. We observe that

$\begin{array}[]{ll}f_{\kappa}(t)-I_{0}e^{-\nu\,\left(t-t_{0}\right)}\beta\left(t-t_{0}\right)&=I_{0}\,e^{-\nu\,\left(t-t_{0}\right)}\displaystyle\int_{0}^{+\infty}\left[\beta\left(a+\left(t-t_{0}\right)\right)-\beta\left(t-t_{0}\right)\right]\kappa\,e^{-\kappa
a}da,\vspace{0.2cm}\\\
&=I_{0}\,e^{-\nu\,\left(t-t_{0}\right)}\displaystyle\int_{0}^{\eta}\left[\beta\left(a+\left(t-t_{0}\right)\right)-\beta\left(t-t_{0}\right)\right]\kappa\,e^{-\kappa
a}da,\vspace{0.2cm}\\\ &\hskip
5.69046pt+I_{0}\,e^{-\nu\,\left(t-t_{0}\right)}\displaystyle\int_{\eta}^{+\infty}\left[\beta\left(a+\left(t-t_{0}\right)\right)-\beta\left(t-t_{0}\right)\right]\kappa\,e^{-\kappa
a}da.\end{array}$

Let $t_{1}>t_{0}$ be such that

$I_{0}\,e^{-\nu\,\left(t-t_{0}\right)}\displaystyle\sup_{a\geq
0}\beta(a)\leq\dfrac{\varepsilon}{4},\forall t\geq t_{1}.$

Let $\eta>0$ be such that

$a\leq\eta\Rightarrow|\beta\left(a+t\right)-\beta\left(t\right)|\leq\dfrac{\varepsilon}{2},\forall
t\in[t_{0},t_{1}].$

Then we have

$|I_{0}\,e^{-\nu\,\left(t-t_{0}\right)}\int_{0}^{\eta}\left[\beta\left(a+\left(t-t_{0}\right)\right)-\beta\left(t-t_{0}\right)\right]\kappa\,e^{-\kappa
a}da|\leq\left\\{\begin{array}[]{ll}I_{0}\,e^{-\nu\,\left(t-t_{0}\right)}\,2\,\displaystyle\sup_{a\geq
0}\beta(a),&\text{ if }t\geq t_{1},\vspace{0.2cm}\\\
I_{0}\,e^{-\nu\,\left(t-t_{0}\right)}\int_{0}^{\eta}\dfrac{\varepsilon}{2}\kappa\,e^{-\kappa
a}da&\text{ if }t\in[t_{0},t_{1}],\end{array}\right.$

therefore

$|f_{\kappa}(t)-I_{0}\,e^{-\nu\,\left(t-t_{0}\right)}\beta\left(t-t_{0}\right)|\leq\dfrac{\varepsilon}{2}+|I_{0}\,e^{-\nu\,\left(t-t_{0}\right)}\int_{\eta}^{+\infty}\left[\beta\left(a+\left(t-t_{0}\right)\right)-\beta\left(t-t_{0}\right)\right]\kappa\,e^{-\kappa
a}da|.$

The result from the fact that

$|I_{0}\,e^{-\nu\,\left(t-t_{0}\right)}\int_{\eta}^{+\infty}\left[\beta\left(a+\left(t-t_{0}\right)\right)-\beta\left(t-t_{0}\right)\right]\kappa\,e^{-\kappa
a}da|\leq 2\,I_{0}\,\sup_{a\geq
0}\beta(a)\int_{\eta}^{+\infty}\kappa\,e^{-\kappa a}da$

hence

$|I_{0}e^{-\nu\,\left(t-t_{0}\right)}\int_{\eta}^{+\infty}\left[\beta\left(a+\left(t-t_{0}\right)\right)-\beta\left(t-t_{0}\right)\right]\kappa\,e^{-\kappa
a}da|\leq 2\,I_{0}\sup_{a\geq 0}\beta(a)e^{-\kappa\eta}\to 0,\text{ as
}\kappa\to\infty.\qed$

Define

$\Gamma(a)=e^{-\nu\,a}\beta\left(a\right),\forall a\geq 0.$ (4.2)

Then by using (2.6), the Kermack-McKendrick model can be reformulated for
$t\geq t_{0}$, as the following system

$N_{\kappa}(t)=\tau(t)S_{\kappa}(t)\left[f_{\kappa}(t)+\int_{0}^{t-t_{0}}\Gamma(a)\,N_{\kappa}(t-a)da\right],$

where $f_{\kappa}(t)$ is defined in (4.2), and

$S_{\kappa}(t)=S_{0}-\int_{t_{0}}^{t}N_{\kappa}(\sigma)d\sigma.$

By taking first a formal limit when $\kappa\to+\infty$, we obtain the model
starting from a single cohort of infected.

Kermack-McKendrick model starting from a single cohort of infected Assume that
the initial distribution of infected only contains a single cohort composed of
$I_{0}$ individuals all with age of infection $a=0$ at time $t_{0}$. Then the
flow of new infected $t\to N(t)$ is the unique continuous solution of the
Volterra integral equation
$N(t)=\tau(t)S(t)\left[I_{0}\times\Gamma\left(t-t_{0}\right)+\int_{0}^{t-t_{0}}\Gamma(a)\,N(t-a)da\right],\forall
t\geq t_{0},$ (4.3) where $S(t)$ is obtained from (2.4).

The following theorem says that the model with a single cohort of infected
extends the earlier model of Kermack-McKendrick with initial distribution in
$L^{1}$. This theorem is a consequence of Lemma 4.1 and the continuity of the
semiflow generated by the Volterra integral equation. We refer to Ducrot and
Magal [14] for more results on this topic.

######  Theorem 4.2.

Let Assumption 2.1 be satisfied, and assume in addition that $a\to\beta(a)$ is
continuous. Then

$\lim_{\kappa\to\infty}N_{\kappa}(t)=N(t),$

where the limit is uniform in $t$ on every closed and bounded interval of
$[t_{0},+\infty)$, and the map $t\to N(t)$ is the unique continuous solution
of the Volterra integral equation (4.3)-(4.1).

######  Remark 4.3.

When the initial distribution is a Dirac mass centered at $a=0$, the total
number of infected individuals at time $t$ is

$C(t)=e^{-\nu\,(t-t_{0})}I_{0}+\int_{0}^{t-t_{0}}e^{-\nu\,a}N(t-a)da,\forall
t\geq t_{0},$

and the number of infectious individuals at time $t$ is

$I(t)=\beta\left(t-t_{0}\right)e^{-\nu\,\left(t-t_{0}\right)}I_{0}+\int_{0}^{t-t_{0}}\beta(a)e^{-\nu\,a}N(t-a)da,\forall
t\geq t_{0}.$

Kermack-McKendrick model starting from multiple cohorts of infected Assume
that the initial distribution of infected consists in $n\geq 1$ cohorts of
infected with age of infection $a_{1}<a_{2}<\ldots<a_{n}$ at time $t_{0}$.
That is
$i_{0}(a)=I_{0}^{1}\,\delta_{a_{1}}(a)+\ldots+I_{0}^{n}\,\delta_{a_{n}}(a).$
where $I_{0}^{j}$ is the number of infected in the $j^{th}$-cohort at time
$t_{0}$. Then the flow of infected $t\to N(t)$ satisfies the following
Volterra integral equation
$\begin{array}[]{l}N(t)=\tau(t)S(t)\left[\displaystyle\sum_{j=1}^{n}\Gamma(t-t_{0}+a_{j})\dfrac{I_{0}^{j}}{e^{-\nu
a_{j}}}+\displaystyle\int_{0}^{t-t_{0}}\Gamma(a)\,N(t-a)da\right],\end{array}$
where $S(t)$ is obtained from (2.4).

###  4.3 Basic reproduction number for the extended model (4.3)

In this section, we assume that the transmission $t\to\tau(t)$ is constant
equal to $\tau$, and $t\to S(t)$ is constant equal to $S_{0}$.

Define the daily reproduction numbers

$R_{0}(a)=\tau\times S_{0}\times\Gamma(a)=\tau\times
S_{0}\times\beta\left(a\right)\times e^{-\nu\,a},\forall a\geq 0.$ (4.4)

Assuming that the number of susceptible individuals $t\to S(t)$ is constant
and equal $S_{0}$ in the $N$-equation (4.3), then we obtain

$N(t)=I_{0}\times R_{0}(t-t_{0})+\int_{0}^{t-t_{0}}R_{0}(a)\,N(t-a)da,\hskip
5.69046pt\forall t\geq t_{0}.$ (4.5)

By using the change of variable $s=t-t_{0}$,

$N(s+t_{0})=I_{0}\times R_{0}(s)+\int_{0}^{s}R_{0}(a)\,N(s+t_{0}-a)da,\hskip
5.69046pt\forall s\geq 0.$

Replacing the notation $s$ by $t$, and define

$N_{t_{0}}(t)=N(t+t_{0}),\forall t\geq 0,$

the equation (4.5) becomes

$N_{t_{0}}(t)=\left[I_{0}\times
R_{0}(t)+\int_{0}^{t}R_{0}(a)\,N_{t_{0}}(t-a)da\right],\forall t\geq 0.$ (4.6)

By replacing $N_{t_{0}}(t)$ by the right hand side of (4.6) in the integral
term of (4.6) we obtain

$N_{t_{0}}(t)=I_{0}\times R_{0}(t)+\int_{0}^{t}R_{0}(a)\,I_{0}\times
R_{0}(t-a)da+\int_{0}^{t}R_{0}(a_{0})\int_{0}^{t-a_{0}}R_{0}(a_{1})\,N_{t_{0}}(t-a_{0}-a_{1})da_{1}da_{0}$

and by induction

$N_{t_{0}}(t)=I_{0}R_{0}(t)+I_{0}(R_{0}*R_{0})(t)+I_{0}(R_{0}*R_{0}*R_{0})(t)+\ldots+I_{0}\underbrace{(R_{0}*R_{0}*\ldots*R_{0})}_{\text{n
times}}(t)+\ldots$

where the convolution is defined by

$(U*V)(t)=\int_{0}^{t}U(a)V(t-a)da=\int_{0}^{t}U(t-a)V(a)da.$

We define

$\left(U^{*(2)}\right)(t)=(U*U)(t),$

and for each integer $n\geq 3$,

$\left(U^{*(n)}\right)(t)=\left(U^{*(n-1)}*U\right)(t)=\left(U*U^{*(n-1)}\right)(t)=\underbrace{(U*U*\ldots*U)}_{\text{n
times}}(t).$

We can interpret the $N$-equation (4.5) concretely as follows

$\begin{array}[]{ll}N_{t_{0}}(t)=&\underbrace{I_{0}\,R_{0}(t)}_{\parbox{300.0pt}{\scriptsize\centering
Flow of infected produced at time $t$ by the first $I_{0}$ infected
individuals\<EMAIL_ADDRESS>&+\underbrace{I_{0}\,\left(R_{0}^{*(2)}\right)(t)}_{\parbox{300.0pt}{\scriptsize\centering
Flow of infected produced at time $t$ by the second generation of infected
individuals\<EMAIL_ADDRESS>&+\underbrace{I_{0}\,\left(R_{0}^{*(3)}\right)(t)}_{\parbox{300.0pt}{\scriptsize\centering
Flow of infected produced at time $t$ by the third generation of infected
individuals\<EMAIL_ADDRESS>&+\\\ &\vdots\\\
&+\underbrace{I_{0}\,\left(R_{0}^{*(n)}\right)(t)}_{\parbox{300.0pt}{\scriptsize\centering
Flow of infected produced at time $t$ by the $n^{th}$ generation of infected
individuals\<EMAIL_ADDRESS>&+\\\ &\vdots\end{array}$

Basic reproduction number The total number of the first generation of new
infected produced by a single infected patient with age of infection $a=0$ at
time $t=t_{0}$ is called the basic reproduction number. That is
$R_{0}=\int_{0}^{\infty}R_{0}(a)\,da.$ The flow of the first generation of new
infected produced by a single infected patient have been infected for $a$ days
is called the daily reproduction numbers. When the time scale is one day, the
function $R_{0}(a)$ is also the average daily number of case produced by a
single patient at the age of infection $a$.

######  Proposition 4.4.

The total number of cases produced by the $n^{th}$ generation of infected
resulting from a single infected patient is

$\int_{0}^{\infty}\left(R_{0}^{*(n)}\right)(t)dt=\left(R_{0}\right)^{n}.$

###### Proof.

By using Fubini’s theorem we have

$\int_{0}^{\infty}(R_{0}*V)(t)\,dt=\int_{0}^{\infty}\int_{0}^{t}R_{0}(t-a)V(a)\,da\,dt=\int_{0}^{\infty}\int_{a}^{+\infty}R_{0}(t-a)\,dt\,V(a)\,da$

and by making the change of variable $l=t-a$ we obtain

$\int_{0}^{\infty}\left(R_{0}*V\right)(t)\,dt=\int_{0}^{+\infty}R_{0}(l)\,dl\times\int_{0}^{+\infty}V(a)\,da.$

Replacing $V(t)$ by $(R_{0}^{*(n-1)})(t)$ in the above equation we obtain

$\int_{0}^{\infty}\left(R_{0}^{*(n)}\right)(t)dt=\int_{0}^{+\infty}R_{0}(l)\,dl\int_{0}^{\infty}\left(R_{0}^{*(n-1)}\right)(t)dt,$

and the result follows by induction. ∎

##  5 Computing the age dependent reproduction number $\Gamma(a)$ from the
data

By using (4.3), we obtain the following result.

Computing $\Gamma(a)$ from the data Assume in addition that the parameters
$t_{0},S_{0}>0,I_{0},\nu>0,$ and the function $t\to\tau(t)$ are given. Then
the function $t\to\Gamma(t)$ can be obtained from the flow of new infected
$t\to N(t)$, as the unique solution of the Volterra integral equation
$\Gamma(t-t_{0})=\dfrac{1}{I_{0}}\left(\dfrac{N(t)}{\tau(t)S(t)}-\int_{0}^{t-t_{0}}\Gamma(a)\,N(t-a)da\right),\forall
t\geq t_{0},$ (5.1) where $S(t)$ is obtained by using (2.4).

######  Remark 5.1.

Assume that patients can not transmit the pathogen when the age of infection
is above $a^{+}>0$. That is

$\Gamma(a)=0,\forall a\geq a^{+}.$

Then the equation (5.1) becomes for all $t\geq t_{0}+a^{+}$,

$\dfrac{N(t)}{\tau(t)S(t)}=\int_{0}^{a^{+}}\Gamma(a)\,N(t-a)da\Leftrightarrow
N(t)=\tau(t)S(t)\int_{0}^{a^{+}}\Gamma(a)\,N(t-a)da.$

###  5.1 Some explicit examples of $a\to R_{0}(a)$

We assume that the transmission $t\to\tau(t)$ is constant equal to $\tau$, and
$t\to S(t)$ is constant equal to $S_{0}$. Then

$R_{0}(a)=\tau\times S_{0}\times\Gamma(t-t_{0}),$

and by setting $a=t-t_{0}$, the equation (5.1) becomes

$R_{0}(a)=\dfrac{N(a+t_{0})}{I_{0}}-\dfrac{1}{I_{0}}\int_{0}^{a}R_{0}(\sigma)\,N(a-\sigma+t_{0})d\sigma,\forall
a\geq 0.$ (5.2)

####  5.1.1 Exponential decay

Assume that

$N(t)=\chi_{1}\,I_{0}\,e^{-\chi_{2}\left(t-t_{0}\right)}.$

Then we obtain

$R_{0}(a)=\chi_{1}\,e^{-\chi_{2}a}-\int_{0}^{a}R_{0}(\sigma)\,\chi_{1}\,e^{-\chi_{2}\left(a-\sigma\right)}d\sigma,\forall
a\geq 0.$

Therefore

$R_{0}^{\prime}(a)=-\left(\chi_{1}+\chi_{2}\right)R_{0}(t),\forall a\geq
0,\text{ and }R_{0}(0)=\chi_{1},$

and we obtain

$R_{0}(a)=\chi_{1}e^{-\left(\chi_{1}+\chi_{2}\right)a},\forall a\geq 0.$ (5.3)

####  5.1.2 Exponential decay with latency

Let $t_{1}>t_{0}$. Assume that

$N(t)=\left\\{\begin{array}[]{l}I_{0}\,\chi_{1}\,e^{-\chi_{2}\left(t-t_{0}\right)},\text{
if }t\geq t_{1},\vspace{0.2cm}\\\ 0,\text{ if }t_{0}\leq
t<t_{1}.\end{array}\right.$

Then we obtain

$R_{0}(a)=0,\forall a\in\left[0,t_{1}-t_{0}\right],$

and

$R_{0}(a)=\chi_{1}\,e^{-\chi_{2}a}-\int_{t_{1}-t_{0}}^{a}R_{0}(\sigma)\,\chi_{1}\,e^{-\chi_{2}\left(a-\sigma\right)}\,d\sigma,\forall
t\geq t_{1}.$

By setting $a_{1}=t_{1}-t_{0}$, we obtain

$R_{0}(a)=0,\forall a\in\left[0,a_{1}\right),$
$R_{0}(a)=\chi_{1}\,e^{-\chi_{2}a}-\int_{a_{1}}^{a}R_{0}(\sigma)\,\chi_{1}\,e^{-\chi_{2}\left(a-\sigma\right)}d\sigma,\forall
a\geq a_{1},$

hence by setting $a=\widehat{a}+a_{1}$,

$R_{0}(\widehat{a}+a_{1})=\chi_{1}\,e^{-\chi_{2}\left(\widehat{a}+a_{1}\right)}-\int_{a_{1}}^{\widehat{a}+a_{1}}\chi_{1}\,e^{-\chi_{2}\left(\widehat{a}+a_{1}-\sigma\right)}\,R_{0}(\sigma)d\sigma,\forall\widehat{a}\geq
0,$

and by setting $\widehat{a}=a-a_{1}$, and
$\widehat{R}_{0}(\widehat{a})=R_{0}(\widehat{a}+a_{1})$, we obtain

$\widehat{R}_{0}(\widehat{a})=C_{1}e^{-\chi_{2}\widehat{a}}-\int_{0}^{\widehat{a}}\chi_{1}\,e^{-\chi_{2}\left(\widehat{a}-\sigma\right)}\,\widehat{R}_{0}(\sigma)d\sigma,\forall\widehat{a}\geq
0,$

where

$C_{1}=\chi_{1}\,e^{-\chi_{2}\,a_{1}}.$

Thus

$\widehat{R}_{0}(\widehat{a})=C_{1}e^{-\left(\chi_{1}+\chi_{2}\right)\widehat{a}},\forall\widehat{a}\geq
0.$

We conclude that

$R_{0}(a)=\left\\{\begin{array}[]{ll}0,&\text{ if }0\leq a\leq
a_{1},\vspace{0.2cm}\\\ \chi_{3}e^{-\left(\chi_{1}+\chi_{2}\right)\,a},&\text{
if }a>a_{1},\end{array}\right.$ (5.4)

where

$\chi_{3}=\chi_{1}\,e^{\chi_{1}\,a_{1}},\text{ and }a_{1}=t_{1}-t_{0}.$

##  6 Day by day Kermack-McKendrick model with age of infection

The variation of the number of susceptible individuals $S(t)$ is given each
day $t=t_{0},t_{0}+1,\ldots$, by

$S(t)=S_{0}-\sum_{d=t_{0}}^{t-1}N(d),$ (6.1)

where $S_{0}$ is the number of susceptible on day $0$, $S(t)$ is the number of
susceptible on day $t$ and $N(d)$ is the daily number of new infected
individuals on day $d$. By analogy with the equation (2.6), the daily number
of new infected individuals satisfies the following discrete time Volterra
integral equation for all $\forall t=t_{0},t_{0}+1,t_{0}+2,\ldots,$

$N(t)=\tau(t)S(t)\sum_{d=t-t_{0}}^{+\infty}\Gamma(d)\dfrac{I_{0}(d-(t-t_{0}))}{e^{-\nu\,(d-(t-t_{0}))}}+\tau(t)S(t)\sum_{d=1}^{t-t_{0}}\Gamma(d)\times
N(t-d),$ (6.2)

where

$\Gamma(d):=\beta(d)\,e^{-\nu\,d},\forall d=0,1,2,\ldots,$

and $\tau(t)$ is the transmission rate, $\beta(d)$ is the probability to be
infectious (i.e. capable to transmit the pathogen) after $d$ days of infection
and $e^{-\nu\,d}$ is the probability to stay infected after $d$ days of
infection (i.e. the probability not to recover or to die after $d$ days of
infection). The quantities $I_{0}(d)$ is the number of infected on day $0$
which have being infected $d$ days ago.

The model (4.3) with a single cohort of infected becomes

$N(t)=\tau(t)S(t)\left[\Gamma(t-t_{0})I_{0}(0)+\sum_{d=1}^{t-t_{0}}\Gamma(d)\times
N(t-d)\right].$ (6.3)

Day by day single cohort model and daily basic reproduction number Assume that
$t\to\tau(t)$ equal $\tau_{0}$, and $t\to S(t)$ is constant equal to $S_{0}$.
Assume that the epidemic starts at time $t_{0}$ with a cohort of $I_{0}$ new
infected patient (i.e. with age of infection $a=0$). The model (4.5) with a
single cohort of infected becomes a discrete Volterra equation
$N(t)=\left[R_{0}(t-t_{0})\times I_{0}+\sum_{d=1}^{t-t_{0}}R_{0}(d)\times
N(t-d)\right],\forall t\geq t_{0}.$ (6.4) We obtain the day by day equation
for the daily reproduction number
$R_{0}(a)=\dfrac{N(t_{0}+a)}{I_{0}}-\dfrac{1}{I_{0}}\sum_{d=1}^{a}R_{0}(d)\times
N(t_{0}+a-d),\forall a\geq 0.$ (6.5)

##  7 Numerical simulations

###  7.1 Comparison of deterministic and stochastic simulations

In the simulations, the unit of time is one day, and we fix

$S_{0}=10^{7}=10\,000\,000,\;1/\nu=9\text{ days},\text{ and }R_{0}=1.1.$

For each function $\beta(a)$ described below, the parameter $\tau$ is obtained
numerically by using the following formula

$\tau=\dfrac{R_{0}}{S_{0}\int_{0}^{\infty}\beta(a)e^{-\nu a}da},$

where the integral is computed by using the Simpson integration method.

In the following, we use the numerical scheme described in Appendix A to run
the simulation of the Volterra integral equation (4.3)-(4.1). We use the
Individual Based Model (IBM) described in Appendix B to run the stochastic
simulations of the model. In the following, we illustrate the convergence of
the IBM to the deterministic model whenever $I_{0}$ increases.

####  7.1.1 Example 1

We assume that the probability to be infectious is a shifted gamma like
distribution. That is

$\beta(a)=\beta_{0}(a-a_{0})^{+}e^{-\beta_{1}(a-a_{0})^{+}},$ (7.1)

with

$a_{0}=3\text{ days},\beta_{0}=e/2=1.3591,\text{ and }\beta_{1}=1/2=0.5.$

Figure 1:  On the left-hand side, we plot the function $a\to\beta(a)$. On the
right-hand side, we plot the function $a\to R_{0}(a)=\tau_{0}\times
S_{0}\times\beta(a)\times e^{-\nu a}$.

In the Figures 2-3 we use the IBM to investigate some properties of the
clusters obtained from the stochastic simulations. We compare such a
stochastic sample with the original $a\to R_{0}(a)$.

By comparing Figures 4-5 (for $I_{0}=10$) and Figures 9-10 (for
$I_{0}=1\,000$), we observe the convergence of the IBM to the deterministic
model.

Then in Figures 6-8 (for $I_{0}=10$) and Figures 11-13 (for $I_{0}=1\,000$),
we apply the discrete-time equation (6.5) to reconstruct $a\to R_{0}(a)$ from
the trajectories of the deterministic or stochastic models. In the
deterministic model, we observe the effect of the day-by-day discretization
(which corresponds to the daily reported data). In the stochastic case, we
observe the effect of the stochasticity of the IBM.

First generation of secondary cases produced by a single infected:

(a)

(b)

Figure 2:  In these figures, we present sets of 500 samples of secondary cases
produced by a single infected individual in a population of $S=10^{7}$
susceptible hosts. Theses samples are produced by using the IBM. (a)
Statistical summary: the blue curve represents the average number of cases at
age of infection $a$; the dark blue area is the 95% confidence interval of
this average obtained by fitting a Gaussian distribution to the data; the
light blue area corresponds to the standard deviation; the orange curve is the
deterministic daily basic reproductive number at age $a$. (b) Bar graph of the
average number of secondary cases as a function of the age since infection.
Figure 3:  In these figures, we present sets of 500 samples of secondary cases
produced by a single infected individual in a population of $S=10^{7}$
susceptible hosts. Theses samples are produced by using the IBM. We plot a
histogram of the total number of secondary cases produced during the whole
infection. This estimates the probability of a single infected to generate $n$
secondary cases (with $n$ in the abscissa).

Simulations for $I_{0}=10$:

Figure 4:  On the left-hand side, we plot the function $t\to N(t)$ solution of
(4.3) with (2.4). On the right-hand side, we plot the function
$t\to\int_{t-1}^{t}N(s)ds$ (for $t=1,2,\ldots$) which corresponds to the daily
number of cases obtained from by solving (4.3) with (2.4), and we compare it
with the daily number of cases obtained from $500$ runs of the IBM.

Figure 5:  On the left-hand side, we plot the function
$t\to\int_{0}^{t}N(s)ds$ (for $t=0,1,2,\ldots$) which corresponds to the
cumulative number of cases obtained from by solving (4.3) with (2.4), and we
compare it with the cumulative number of cases obtained from $500$ runs of the
IBM. On the right-hand side, we plot the average values of the $500$ runs
obtained from the IBM as well as the quantiles ($10\%-90\%$ (light blue) and
$25\%-75\%$ (blue)).

In Figures 6-8 we focus on the reconstruction of the daily reproduction number
$R_{0}(a)=\tau\,S_{0}\,e^{-\nu\,a}\,\beta\left(a\right)$. In Figure 6, we
focus on the reconstruction of the daily reproduction number from
deterministic simulations, while in Figures 7-8 we focus on the reconstruction
of the daily reproduction number from stochastic simulations.

Figure 6:  On the left hand side, we plot the daily basic reproduction number
by using the original formula $R_{0}(a)$ with (7.1). On the right-hand side,
we apply formula (6.5) to the flow of new infected obtained from the
deterministic model. We vary $I_{0}=6,10,14$. The value $I_{0}=10$ corresponds
to the value used for the simulation of the deterministic model. The yellow
curve gives the best visual fit, and the $R_{0}(a)$ becomes negative whenever
$I_{0}$ becomes too small.

Figure 7:  On the left-hand side, we plot the daily number of cases $t\to
N(t)$ (for $t=0,1,2,\ldots$) obtained from a single run of the IBM. On the
right-hand side, we apply formula (6.5) (with $I_{0}=10$) to the daily number
of cases of case obtained from the IBM.

Figure 8:  On the left-hand side, we plot the daily number of cases $t\to
N(t)$ (for $t=0,1,2,\ldots$) obtained by summing the daily number of cases for
$500$ runs of the IBM. On the right-hand side, we apply formula (6.5) (with
$I_{0}=500\times 10$) to the daily number of cases obtained from the IBM.

Simulations for $I_{0}=1000$:

Figure 9:  On the left-hand side, we plot the function $t\to N(t)$ solution of
(4.3) with (2.4). On the right-hand side, we plot the function
$t\to\int_{t-1}^{t}N(s)ds$ (for $t=1,2,\ldots$) which corresponds to the daily
number of cases obtained from by solving (4.3) with (2.4), and we compare it
with the daily number of cases obtained from $500$ runs of the IBM.

Figure 10:  On the left-hand side, we plot the function
$t\to\int_{0}^{t}N(s)ds$ (for $t=0,1,2,\ldots$) which corresponds to the
cumulative number of cases obtained from by solving (4.3) with (2.4), and we
compare it with the cumulative number of cases obtained from $500$ runs of the
IBM. On the right-hand side, we plot the average values of the $500$ runs
obtained from the IBM as well as the quantiles ($10\%-90\%$ (light blue) and
$25\%-75\%$ (blue)).

In Figures 11-13 we focus on the reconstruction of the daily reproduction
number $R_{0}(a)=\tau\,S_{0}\,e^{-\nu\,a}\,\beta\left(a\right)$. In Figure 11,
we focus on the reconstruction of the daily reproduction number from
deterministic simulations, while in Figures 12-13 we focus on the
reconstruction of the daily reproduction number from stochastic simulations.

Figure 11:  On the left hand side, we plot the daily basic reproduction number
by using the original formula $R_{0}(a)$ with (7.1). On the right-hand side,
we apply formula (6.5) to the flow of new infected obtained from the
deterministic model. We vary $I_{0}=600,1000,1400$. The value $I_{0}=1000$
corresponds to the value used for the simulation of the deterministic model.
The yellow curve gives the best visual fit, and the $R_{0}(a)$ becomes
negative whenever $I_{0}$ becomes too small.

Figure 12:  On the left-hand side, we plot the daily number of cases $t\to
N(t)$ (for $t=0,1,2,\ldots$) obtained from a single run of the IBM. On the
right-hand side, we apply formula (6.5) (with $I_{0}=1000$) to the daily
number of cases of case obtained from the IBM.

Figure 13:  On the left-hand side, we plot the daily number of cases $t\to
N(t)$ (for $t=0,1,2,\ldots$) obtained by summing the daily number of cases for
$500$ runs of the IBM. On the right-hand side, we apply formula (6.5) (with
$I_{0}=500\times 1000$) to the daily number of cases obtained from the IBM.

####  7.1.2 Example 2

It is common to see biphasic flu clinically: after incubation of one day,
there is a high fever, then a drop in temperature before rising again, hence
the term "V" fever [7]. Such a biphasic contagiousness is also observed in
Covid-19. The viral load in throat swab and sputum has been measured for
Covid-19 patients, which leads to biphasic contagiousness [10, 28]. To cover
these type of infectious diseases, we introduce the following form for the
probability to be infectious

$\beta(a)=0.5\times
4q\left\\{(a-a_{0})\left(1-q(a-a_{0})\right)\right\\}^{+}+4q\left\\{(a-pa_{0})\left(1-q(a-pa_{0})\right)\right\\}^{+},$
(7.2)

with

$a_{0}=3\text{ days},p=2.5,\text{ and }q=0.3.$

Figure 14:  On the left-hand side, we plot the function $a\to\beta(a)$. On the
right-hand side, we plot the function $a\to R_{0}(a)=\tau_{0}\times
S_{0}\times\beta(a)\times e^{-\nu a}$.

In the Figures 15-16 we use the IBM to investigate some properties of the
clusters obtained from the stochastic simulations. We compare such a
stochastic sample with the original $a\to R_{0}(a)$.

By comparing Figures 17-18 (for $I_{0}=10$) and Figures 22-23 (for
$I_{0}=1\,000$), we observe the convergence of the IBM to the deterministic
model.

Then in Figures 19-21 (for $I_{0}=10$) and Figures 24-26 (for $I_{0}=1\,000$),
we apply the discrete-time equation (6.5) to reconstruct $a\to R_{0}(a)$ from
the trajectories of the deterministic or stochastic models. In the
deterministic model, we observe the effect of the day-by-day discretization
(which corresponds to the daily reported data). In the stochastic case, we
observe the effect of the stochasticity of the IBM.

First generation of secondary cases produced by a single infected:

(a)

(b)

Figure 15:  In these figures, we present sets of 500 samples of secondary
cases produced by a single infected individual in a population of $S=10^{7}$
susceptible hosts. Theses samples are produced by using the IBM. (a)
Statistical summary: the blue curve represents the average number of cases at
age of infection $a$; the dark blue area is the 95% confidence interval of
this average obtained by fitting a Gaussian distribution to the data; the
light blue area corresponds to the standard deviation; the orange curve is the
deterministic daily basic reproductive number at age $a$. (b) Bar graph of the
average number of secondary cases as a function of the age since infection.
Figure 16:  In these figures, we present sets of 500 samples of secondary
cases produced by a single infected individual in a population of $S=10^{7}$
susceptible hosts. Theses samples are produced by using the IBM. We plot a
histogram of the total number of secondary cases produced during the whole
infection. This estimates the probability of a single infected to generate $n$
secondary cases (with $n$ in the abscissa).

Simulations for $I_{0}=10$:

Figure 17:  On the left-hand side, we plot the function $t\to N(t)$ solution
of (4.3) with (2.4). On the right-hand side, we plot the function
$t\to\int_{t-1}^{t}N(s)ds$ (for $t=1,2,\ldots$) which corresponds to the daily
number of cases obtained from by solving (4.3) with (2.4), and we compare it
with the daily number of cases obtained from $500$ runs of the IBM.

Figure 18:  On the left-hand side, we plot the function
$t\to\int_{0}^{t}N(s)ds$ (for $t=0,1,2,\ldots$) which corresponds to the
cumulative number of cases obtained from by solving (4.3) with (2.4), and we
compare it with the cumulative number of cases obtained from $500$ runs of the
IBM. On the right-hand side, we plot the average values of the $500$ runs
obtained from the IBM as well as the quantiles ($10\%-90\%$ (light blue) and
$25\%-75\%$ (blue)).

In Figures 19-21 we focus on the reconstruction of the daily reproduction
number $R_{0}(a)=\tau\,S_{0}\,e^{-\nu\,a}\,\beta\left(a\right)$. In Figure 19,
we focus on the reconstruction of the daily reproduction number from
deterministic simulations, while in Figures 20-21 we focus on the
reconstruction of the daily reproduction number from stochastic simulations.

Figure 19:  On the left hand side, we plot the daily basic reproduction number
by using the original formula $R_{0}(a)$ with (7.2). On the right-hand side,
we apply formula (6.5) to the flow of new infected obtained from the
deterministic model. We vary $I_{0}=6,10,14$. The value $I_{0}=10$ corresponds
to the value used for the simulation of the deterministic model. The yellow
curve gives the best visual fit, and the $R_{0}(a)$ becomes negative whenever
$I_{0}$ becomes too small.

Figure 20:  On the left-hand side, we plot the daily number of cases $t\to
N(t)$ (for $t=0,1,2,\ldots$) obtained from a single run of the IBM. On the
right-hand side, we apply formula (6.5) (with $I_{0}=10$) to the daily number
of cases of case obtained from the IBM.

Figure 21:  On the left-hand side, we plot the daily number of cases $t\to
N(t)$ (for $t=0,1,2,\ldots$) obtained by summing the daily number of cases for
$500$ runs of the IBM. On the right-hand side, we apply formula (6.5) (with
$I_{0}=500\times 10$) to the daily number of cases obtained from the IBM.

Simulations for $I_{0}=1000$:

Figure 22:  On the left-hand side, we plot the function $t\to N(t)$ solution
of (4.3) with (2.4). On the right-hand side, we plot the function
$t\to\int_{t-1}^{t}N(s)ds$ (for $t=1,2,\ldots$) which corresponds to the daily
number of cases obtained from by solving (4.3) with (2.4), and we compare it
with the daily number of cases obtained from $500$ runs of the IBM.

Figure 23:  On the left-hand side, we plot the function
$t\to\int_{0}^{t}N(s)ds$ (for $t=0,1,2,\ldots$) which corresponds to the
cumulative number of cases obtained from by solving (4.3) with (2.4), and we
compare it with the cumulative number of cases obtained from $500$ runs of the
IBM. On the right-hand side, we plot the average values of the $500$ runs
obtained from the IBM as well as the quantiles ($10\%-90\%$ (light blue) and
$25\%-75\%$ (blue)).

In Figures 24-26 we focus on the reconstruction of the daily reproduction
number $R_{0}(a)=\tau\,S_{0}\,e^{-\nu\,a}\,\beta\left(a\right)$. In Figure 24,
we focus on the reconstruction of the daily reproduction number from
deterministic simulations, while in Figures 25-26 we focus on the
reconstruction of the daily reproduction number from stochastic simulations.

Figure 24:  On the left hand side, we plot the daily basic reproduction number
by using the original formula $R_{0}(a)$ with (7.2). On the right-hand side,
we apply formula (6.5) to the flow of new infected obtained from the
deterministic model. We vary $I_{0}=600,1000,1400$. The value $I_{0}=1000$
corresponds to the value used for the simulation of the deterministic model.
The yellow curve gives the best visual fit, and the $R_{0}(a)$ becomes
negative whenever $I_{0}$ becomes too small.

Figure 25:  On the left-hand side, we plot the daily number of cases $t\to
N(t)$ (for $t=0,1,2,\ldots$) obtained from a single run of the IBM. On the
right-hand side, we apply formula (6.5) (with $I_{0}=1000$) to the daily
number of cases of case obtained from the IBM.

Figure 26:  On the left-hand side, we plot the daily number of cases $t\to
N(t)$ (for $t=0,1,2,\ldots$) obtained by summing the daily number of cases for
$500$ runs of the IBM. On the right-hand side, we apply formula (6.5) (with
$I_{0}=500\times 1000$) to the daily number of cases obtained from the IBM.

##  8 Application to SARS-CoV-1

In practice, the Kermack-McKendrick model starting from a Dirac mass means
that the epidemic starts from a single patient at time $t_{0}$ (whenever
$I_{0}=1$) or from a group of $I_{0}$ infected patients all with the same age
of infection $a=0$ at time $t_{0}$. This assumption corresponds to the
standard conception of a cluster in epidemiology. An example of such a cluster
is obtained [34] for the SARS-CoV-1 epidemic in Singapore in 2003. The cluster
is represented by a network of contact between individuals in Figure 27.

Figure 27:  On the left-hand side, we plot the contact network of the five
super spreader cases in the SARS epidemic in Singapore in 2003 [34]. The super
spreaders are patient 1, patient 6, patient 35, patient 130 and patient 127.
On the right-hand side, we plot the daily reported cases from Singapore for
the epidemic of SARS in 2003. Case 1 generated 21 cases and 3 suspected cases,
case 2 generated 23 cases and 5 suspected cases, case 3 generated 23 cases and
18 suspected cases, case 4 generated 40 cases and 22 suspected cases, case 5
generated 15 cases and 0 suspected cases [34]. The cases 1,2,3,4,5 correspond
respectively to the patients 1, 6, 35, 130 and 127.

Figure 27 presents the time series of reported cases by source of infection
and date of fever onset. In Figure 28 we present three representations of
these data in continuous time: as a step function, regularized by Gaussian
average and rolling weekly average. In Figure 29 we apply the continuous-time
model to the rolling weekly regularization of the data. Similar to the
reconstruction of $R_{0}(a)$ presented in Figures 6-8, 11-13, 19-21, 24-26,
the basic reproduction number $R_{0}(a)$ becomes negative after a given age.
Our interpretation is that the data are far from perfect and involves sampling
errors and probably a large number of undetected cases. The fact that the
transmission rate is subject to variations in time could also explain this
negativity. In Figure 30 we apply the discrete model (1.6) to the original
data for different values of $I_{0}$. Finally, in Figure 31, we transform the
data by taking advantage of the information on the source of infection given
in [34]. We fix an incubation period of $5$ days, which corresponds to the
average incubation period reported in [34]. Then we shift all secondary cases
produced by the six sources identified in the article [34] to the same origin,
as if all cases had been produced by the same cluster of six individuals. We
present the data on the left-hand side of Figure 31 and apply the method to
obtain $R_{0}(a)$ for parameters $I_{0}=30$, $I_{0}=50$, $I_{0}=100$.

Figure 28:  Regularizations of the daily cases data from the SARS-CoV-1
outbreak in Singapore [34]. The applications in Figure 29 are done with the
“Rolling Weekly” regularization.

Figure 29:  Left: Regularized data of the SARS-CoV-1 outbreak in Singapore in
2003 [34] (black line) and the numerical solution of the model (1.1) with
$I_{0}=30$ and $R_{0}(a)$ computed by (1.6) (blue line). The solutions $N(t)$
of the model (1.1) with $I_{0}=50$ and $I_{0}=100$ are exactly the same when
we use the corresponding $R_{0}(a)$, therefore they are not represented here.
Right: numerical solution of the $R_{0}(a)$ function computed by using the
continuous model (1.4) with $I_{0}=30$, $I_{0}=50$ and $I_{0}=100$.

Figure 30:  Left: Original daily reported cases data of the SARS-CoV-1
outbreak in Singapore in 2003 [34] (blue line). Right: Numerical solution of
the $R_{0}(a)$ function computed by using the discrete model (1.6) with
$I_{0}=30$, $I_{0}=50$ and $I_{0}=100$.

Figure 31:  Left: Daily reported cases of a theoretical cluster based on the
data from the SARS-CoV-1 outbreak in Singapore in 2003 [34] (blue line). The
secondary cases produced by the six patients identified in [34] are shifted to
the same initial infection date of the infector, as if all were produced by
the same cluster starting at $t=0$. We also use an incubation period of 5
days, as indicated in [34]. Right: Numerical solution of the $R_{0}(a)$
function computed by using the discrete model (1.6) with $I_{0}=30$,
$I_{0}=50$ and $I_{0}=100$.

##  9 Discussion

We see from the numerical simulations in Section 7 that the initial number of
infected $I_{0}$ has a very significant influence on the value of the basic
daily reproduction numbers $R_{0}(d)$: these decrease sharply with $I_{0}$,
until they become negative and their fluctuations increase in stochastic
simulations. This tendency to negativity for small $I_{0}$ and these
fluctuations are only corrected when the results are averaged for a large
number of stochastic simulations (500). It can also be noted that the
stochastic simulations lead to a behavior of the hyper exponential type in the
coefficient of variation of the secondary cases produced by an infectious
individual, that is to say that it is relatively constant and much greater
than 1. This phenomenon is to be related to the exponential character of the
gamma distribution used in the simulations.

In deterministic simulations, one can observe a great dependence of the
results with regard to the value of $I_{0}$ concerning the reconstruction of
the daily basic reproduction numbers $R_{0}(a)$ along the contagiousness
period of an individual. In particular, the set of days at which $R_{0}(a)$
fluctuates becoming in some cases negative is important for the great values
of $I_{0}$. In stochastic simulations, we observe the same behavior for the
different curves related to the $R_{0}(a)$ curves sample, but the expectation
of this curves sample considerably attenuates these fluctuations and the
coefficient of variation of the curves remains approximately constant, while
being greater than 1, as in the case of hyperexponential distributions, in
agreement with the exponential character of the part D of the equation (1.7)
defining $R_{0}(a)$.

Concerning the clusters, from observations made during investigations of the
start of the outbreak in some countries [33, 1, 11, 29, 15, 21, 5, 18, 30, 6,
20, 19, 24, 31], it is possible to get spatial and temporal information on the
start of the epidemic, but these studies rarely allow the estimation of the
parameters $S_{0}$ and $\tau$ in the concerned population and worse, they give
no indication of how long they remain constant. Here we assumed that they
remained constant only during the period of exponential growth of new cases
observed.

Our work provides a method to reconstruct the daily basic reproduction number
$R_{0}(a)$ from the daily reported cases data, as long as we consider a
cluster starting from a single infected. This is a strong assumption which is
usually neglected. It is extremely hard to find in the literature a datasets
which satisfies this assumption. For COVID-19, we did not find any publication
including suitable data. While not published yet, we believe that this kind of
data could be gathered by a detailed contact-tracing and – duly anonymized –
could be made available by request. That would allow the development of more
realistic and accurate methods for the analysis and forecast of epidemics.

Appendix

##  Appendix A Euler approximation for the Volterra integral equation

We use a numerical scheme for the equation (4.3) and (4.1).

$N(t)=\tau\,S\left(\max\left(t-\Delta
t,t_{0}\right)\right)\left[I_{0}\times\Gamma\left(t-t_{0}\right)+\int_{0}^{\max\left(t-\Delta
t-t_{0},0\right)}\Gamma(a)\,N(t-a)da\right],$ (A.1)

with

$S(t)=S_{0}-\int_{t_{0}}^{t}N(\sigma)d\sigma,$ (A.2)

and we use the Simpson’s rule to compute the integrals.

##  Appendix B Stochastic simulations: Individual Based Model

In order to estimate the uncertainty expected in real datasets, we use
stochastic simulations that reproduce the first stages of the epidemic in
finite populations. We consider a population composed of a finite number
$N=S_{0}+I_{0}$ of individuals. We start the simulation a time $t=0$ with
$S_{0}\in\mathbb{N}$ susceptible individuals and $I_{0}\in\mathbb{N}$ infected
individuals all with age of infection $a=0$. For each infected individuals we
also compute the time spent in the $I$-compartment which follows an
exponential law with parameters $1/\nu$. The principles of the simulations are
as follows:

  1. 1.

Individuals meet at random at rate $\tau>0$. In other words, each pair of
individual in the population has a contact which occurs at a time following an
exponential law of average $1/\tau$.

  2. 2.

When a contact occurs between an infected individual of age $a$ and a
susceptible individual, the contact results in a newly infected individual of
age $0$ with probability $\beta(a)$. When the infection occurs, the newly
infected individual is assigned a duration of infection which follows an
exponential law of rate $\nu$. Therefore individuals stay infected on average
for a duration of $1/\nu$.

  3. 3.

The age of all individuals is updated at fixed intervals of time of size
$\Delta t$. Simultaneously the life-span of each infected inviduals is
decreased by $\Delta t$ and individual whose life-span has become negative are
removed from the system.

The MATLAB code of the IBM is available online at:

https://github.com/romainvieme/2022-kermack-mckendrick-single-cohort.

## References

  * [1] D. C. Adam, et al., Clustering and superspreading potential of SARS-CoV-2 infections in Hong Kong, Nature Medicine, 26.11 (2020), 1714-1719.
  * [2] L. Alvarez, M. Colom, J. D. Morel, & J. M. Morel, Computing the daily reproduction number of COVID-19 by inverting the renewal equation using a variational technique. Proceedings of the National Academy of Sciences, 118(50) (2021).
  * [3] L. Alvarez, J-D. Morel, J-M. Morel, Modeling Covid-19 incidence by the renewal equation after removal of administrative bias and noise, Biology 11(4), 540 (2022).
  * [4] D. Bernoulli, Essai d”une nouvelle analyse de la mortalité causé par la petite vérole et des avantages de l’inoculation pour la prévenir, Mém. Math. Phys. Acad. Roy. Sci. Paris, (1760), 1-45.
  * [5] M. M. Böhmer, U. Buchholz, V. M. Corman, M. Hoch, K. Katz, D. V. Marosevic, et al. Investigation of a COVID-19 outbreak in Germany resulting from a single travel-associated primary case: a case series. The Lancet Infectious Diseases, 20.8 (2020), 920-928.
  * [6] T. C. Chan, and C. C. King. Surveillance and epidemiology of infectious diseases using spatial and temporal lustering methods. In: Infectious disease informatics and biosurveillance (pp. 207-234). Springer, Boston, MA (2011).
  * [7] D. L. Chao, M. E. Halloran, V. J. Obenchain, & Jr, I. M Longini, FluTE, a publicly available stochastic influenza epidemic simulation model. PLoS computational biology, 6(1) (2010), e1000656.
  * [8] J. Demongeot, Q. Griette, and P. Magal.  SI epidemic model applied to COVID-19 data in mainland China.  R. Soc. Open Sci. 7.12 (2020), 201878.
  * [9] J. Demongeot, Q. Griette, P. Magal, and G. Webb. Modeling vaccine efficacy for COVID-19 outbreak in New York city. Biology, 11.3 (2022), 345.
  * [10] J. Demongeot, K. Oshinubi, M. Rachdi, H. Seligmann, F. Thuderoz & J. Waku, Estimation of Daily Reproduction Rates in COVID-19 Outbreak. Computation, 9, 109 (2021).
  * [11] M. R. Desjardins, A. Hohl, and E. M. Delmelle. Rapid surveillance of COVID-19 in the United States using a prospective space-time scan statistic: Detecting and evaluating emerging clusters. Applied geography 118 (2020), 102202.
  * [12] K. Dietz and J. A. P. Heesterbeek, Daniel Bernoulli’s epidemiological model revisited, Math. Biosci., 180, (2002), 1-21.
  * [13] K. Dietz and J. A. P. Heesterbeek, Daniel Bernoulli’s epidemiological model revisited, Math. Biosci., 180, 1-21, (2002).
  * [14] A. Ducrot, and P. Magal, Volterra integral equations in the framework of integrated semigroups (in preparation).
  * [15] T. Ganyani, C. Kremer, D. Chen, A. Torneri, C. Faes, J. Wallinga, and N. Hens Estimating the generation interval for coronavirus disease (COVID-19) based on symptom onset data, March 2020. Eurosurveillance, 25.17 (2020) 2000257.
  * [16] Q. Griette, J. Demongeot, and P. Magal.  A robust phenomenological approach to investigate COVID-19 data for France.  Math. Appl. Sci. Eng., 2021.
  * [17] Q. Griette, J. Demongeot, and P. Magal. What can we learn from COVID-19 data by using epidemic models with unidentified infectious cases?. Mathematical Biosciences and Engineering, 19.1 (2022), 537-594.
  * [18] A. Guttmann, L. Ouchchane, X. Li, I. Perthus, J. Gaudart, J. Demongeot, and J. Y. Boire. Performance map of a cluster detection test using extended power. International Journal of Health Geographics, 12.1 (2013), 1-10.
  * [19] L. Han, P. Shen, J. Yan, Y. Huang, X. Ba, W. Lin, et al. Exploring the Clinical Characteristics of COVID-19 Clusters Identified Using Factor Analysis of Mixed Data-Based Cluster Analysis. Frontiers in medicine, 8 (2021).
  * [20] S. Hisada, T. Murayama, K. Tsubouchi, S. Fujita, S. Yada, S. Wakamiya, and E. Aramaki. Surveillance of early stage COVID-19 clusters using search query logs and mobile device-based location information. Scientific Reports, 10.1 (2020), 1-8.
  * [21] Q. L. Jing, M. J. Liu, Z. B. Zhang, L. Q. Fang, J. Yuan, A. R. Zhang, et al. Household secondary attack rate of COVID-19 and associated determinants in Guangzhou, China: a retrospective cohort study. The Lancet Infectious Diseases, 20.10 (2020), 1141-1150.
  * [22] W. O. Kermack and A. G. McKendrick, Contributions to the mathematical theory of epidemics: II, Proc. R. Soc. Lond. Ser. B, 138, (1932), 55-83.
  * [23] C. Kremer, T. Ganyani, D. Chen, A. Torneri, C. Faes, J. Wallinga, and N. Hens. Authors’ response: estimating the generation interval for COVID-19 based on symptom onset data. Eurosurveillance, 25.29) (2020), 2001269.
  * [24] A. Ladoy, O. Opota, P. N. Carron, I. Guessous, S. Vuilleumier, S. Joost, and G. Greub. Size and duration of COVID-19 clusters go along with a high SARS-CoV-2 viral load: A spatio-temporal investigation in Vaud state, Switzerland. Science of The Total Environment, 787 (2021), 147483.
  * [25] Z. Liu, P. Magal, O. Seydi, & G. Webb, Understanding unreported cases in the COVID-19 epidemic outbreak in Wuhan, China, and the importance of major public health interventions. Biology, 9(3), 50 (2020).
  * [26] H. Nishiura, Time variations in the transmissibility of pandemic influenza in Prussia, Germany, from 1918-19. Theor. Biol. Med. Model. 4, 20 (2007).
  * [27] H. Nishiura, G. Chowell, “The effective reproduction number as a prelude to statistical estimation of time-dependent epidemic trends” in Mathematical and Statistical Estimation Approaches in Epidemiology, G. Chowell, J. M. Hyman, L. M. A. Bettencourt, C. Castillo-Chavez, Eds. (Springer, Dordrecht, Netherlands, 2009), pp. 103-121.
  * [28] Y. Pan, D. Zhang, P. Yang, L. L. M. Poon, and Q. Wang, Viral load of sars-cov-2 in clinical samples. The Lancet Infectious Diseases, 20(4) (2020), 411–412.
  * [29] R. Pung, C. J. Chiew, B. E. Young, S. Chin, M. I. Chen, H. E. Clapham, et al. Investigation of three clusters of COVID-19 in Singapore: implications for surveillance and response measures. The Lancet, 395.10229 (2020), 1039-1046.
  * [30] F. Shams, A. Abbas, W. Khan, U. S. Khan, and R. Nawaz. A death, infection, and recovery (DIR) model to forecast the COVID-19 spread. Computer Methods and Programs in Biomedicine Update, 2 (2022), 100047.
  * [31] A. Tariq, Y. Lee, K. Roosa, S. Blumberg, P. Yan, S. Ma, and G. Chowell. Real-time monitoring the transmission potential of COVID-19 in Singapore, March 2020. BMC medicine, 18.1 (2020), 1-14.
  * [32] M. Worobey, T. Watts, R. McKay, et al. 1970s and ‘Patient 0’ HIV-1 genomes illuminate early HIV/AIDS history in North America. Nature 539 (2016), 98–101. https://doi.org/10.1038/nature19827
  * [33] S. E. F. Yong, et. al., Connecting clusters of COVID-19: an epidemiological and serological investigation, The Lancet Infectious Diseases 20.7 (2020), 809-815.
  * [34] Centers for Disease Control and Prevention (CDC), Severe acute respiratory syndrome Singapore, MMWR. Morbidity and mortality weekly report 52.18, 2005. https://www.cdc.gov/mmwr/preview/mmwrhtml/mm5218a1.htm