# Final size and partial distance estimate for a two-group SEIRD model

Alison M.V.D.L. Melo  A.M.V.D.L.M.: Universidade Federal do Vale do São
Francisco - UNIVASF, 56304-917, Petrolina, Brazil<EMAIL_ADDRESS>and  Matheus C. Santos  M.C.S.: Departamento de Matemática Pura e Aplicada–
IME, Universidade Federal do Rio Grande do Sul - UFRGS, 91509-900, Porto
Alegre, Brazil<EMAIL_ADDRESS>
(Date: June 23, 2021)

###### Abstract.

In this paper we consider a SEIRD epidemic model for a population composed by
two groups of individuals with asymmetric interaction. Given an approximate
solution for the two-group model, we estimate the error of this approximation
to the unknown solution to the second group based on the known error that the
approximation has with respect to the solution to the first group. We also
study the final size of the epidemic for each group. We illustrate our results
with the spread of the coronavirus disease 2019 (COVID-19) pandemic in the New
York County (USA) for the initial stage of the contamination, and in the
cities of Petrolina and Juazeiro (Brazil).

###### Key words and phrases:

Epidemic mathematical model, latency period, final size, distance of solutions

###### 2010 Mathematics Subject Classification:

92D25, 92D30, 34C60

Submitted on

##  1\. Introduction

The models for outbreak and spread of diseases have a long history of studies
and analysis since the pioneering work of John Graunt, who made a statistical
approach to understand public health and causes of deaths in England ([16]),
and Daniel Bernoulli, who proposed the first mathematical model describing an
infectious disease ([8]) and the impact of inoculation for the smallpox
control ([9]). The models we use nowadays (SIR models and extensions)
originated in the works of Sir Ronald Ross ([32]), who formulated a system of
differential equations after his studies on the spread of malaria, and later,
refined by William Kermack and Anderson McKendrick, who generalized the
approach and established the condition for an epidemic to occur ([19]), which
is now known as the basic reproduction number. Many related models and
approaches have been studied since then ([11, 18, 26, 29]), and were specially
used recently to analyze the spread and to forecast the number of cases of
COVID-19 epidemic ([10, 15, 17, 23, 24, 30, 31, 36, 37] and many others).

In this paper we study the case of a two-group epidemic model, where the
infectious individuals of each group transmit the disease to the susceptible
individuals of both groups according to an asymmetric interaction. The
division in groups may happen for biological, geographic or even socioeconomic
reasons. For instance, in [6] an age structured SIRS model was applied to
study seasonal evolution of Respiratory Syncytial Virus (RSV) in Valencia,
Spain. Infection by RSV tends to be more severe in babies under one year old.
On the other hand, COVID-19 affects elderly people more aggressively ([33,
4]). A disease may also spread differently between locations. In [12] the
authors study the spread of sleeping sickness by using a differential equation
model where both human and vector populations are divided into two patches
according to location: plantation and village. Finally, socioeconomic
inequality might play an important role in the spread of diseases such as
COVID-19 as shown in [13] and [27].

After a susceptible individual becomes infected by interacting with an
infectious one, there is a period of time before he can transmit the disease
himself. This period of latency or exposure may or may not be relevant for the
dynamics of a disease, depending on how long it may be compared to the total
duration of the infectious period. In the case of COVID-19, the latency period
is on average 3-5 days ([21, 38]) up to 14 days, which is 1-2 days less than
the incubation period, i.e., the period needed for the symptoms onset, which
is on average 5-6 days, also up to 14 days ([4]). Therefore, since the
infectious period may vary between 5 days (mild cases) up to 6 weeks (severe
or critical cases) ([4, 22]), this exposure period is relevant to the
progression of the epidemic. That is why in this work we have chosen to study
the two group SEIRD model, which includes the exposure period.

In this article, we study the system of equations for the two-group SEIRD
model. The system considered here is the following:

(1)
$\left\\{\begin{array}[]{l}S^{\prime}(t)=-\mathrm{diag}(S(t))\mathcal{B}I(t)\\\
E^{\prime}(t)=\mathrm{diag}(S(t))\mathcal{B}I(t)-\mathcal{A}E(t)\\\
I^{\prime}(t)=\mathcal{A}E(t)-(\Gamma+\mathcal{M})I(t)\\\ R^{\prime}(t)=\Gamma
I(t)\\\ D^{\prime}(t)=\mathcal{M}I(t)\end{array}\right.$

with initial conditions

$\displaystyle
S(0)=S_{0},\;E(0)=E_{0},\;I(0)=I_{0},\;R(0)=R_{0},\;D(0)=D_{0}\in\mathbb{R}^{2}_{+},$

where the vectorial functions $\displaystyle S(t),E(t),I(t),R(t)$ and
$\displaystyle D(t)$ represent the classes of susceptible, exposed,
infectious, recovered and deceased individuals, respectively, for groups 1 and
2 as

$\displaystyle S(t)\\!=\\!\begin{pmatrix}S_{1}(t)\\\
S_{2}(t)\end{pmatrix},E(t)\\!=\\!\begin{pmatrix}E_{1}(t)\\\
E_{2}(t)\end{pmatrix},I(t)\\!=\\!\begin{pmatrix}I_{1}(t)\\\
I_{2}(t)\end{pmatrix},R(t)\\!=\\!\begin{pmatrix}R_{1}(t)\\\
R_{2}(t)\end{pmatrix},D(t)\\!=\\!\begin{pmatrix}D_{1}(t)\\\
D_{2}(t)\end{pmatrix}.$

Note that, if $\displaystyle N=(N_{1},N_{2})$ is the number of individuals at
each group at time $\displaystyle t=0$, by the system of equations (1) we have

(2) $S(t)+E(t)+I(t)+R(t)+D(t)=N\;,\;\;\;\mbox{ for all }\;t\geqslant 0.$

The transmission of the disease is given by the interaction of the infectious
individuals with the susceptible ones and expressed by the infection matrix,
composed by the infection rates $\displaystyle\beta_{ij}$ of the susceptible
individuals $\displaystyle S_{i}$ by the infectious individuals $\displaystyle
I_{j}$

$\displaystyle\mathcal{B}=\begin{pmatrix}\beta_{11}&\beta_{12}\\\
\beta_{21}&\beta_{22}\end{pmatrix}.$

Once infected, each individual becomes an asymptomatic noninfectious
individual in class $\displaystyle E$. The exit of the exposed class is given
by the latency matrix

$\displaystyle\mathcal{A}=\left(\begin{array}[]{cc}\alpha_{1}&0\\\
0&\alpha_{2}\end{array}\right),$

where $\displaystyle\alpha_{1},\alpha_{2}>0$. This means that each newly
infected individual takes the average time of $\displaystyle 1/\alpha_{i}$
(according to the group) to become infectious, i.e, to be in the class
$\displaystyle I$ and contribute with the communicability of the disease. The
dynamics of individuals leaving the class $\displaystyle I$ is governed by the
recovering rate matrix $\displaystyle\Gamma$ and the death rate matrix
$\displaystyle\mathcal{M}$ given by

$\displaystyle\Gamma=\left(\begin{array}[]{cc}\gamma_{1}&0\\\
0&\gamma_{2}\end{array}\right)\;,\;\;\mathcal{M}=\left(\begin{array}[]{cc}\mu_{1}&0\\\
0&\mu_{2}\end{array}\right).$

Thus the exit flux of infectious individuals in each group is composed by
$\displaystyle\gamma_{i}I_{i}$, which represents those who recovered from the
disease and become immune (per time unit), and $\displaystyle\mu_{i}I_{i}$
which gives the rate of deaths caused by the disease per unit of time. Hence
the average time of infectiousness in each group is $\displaystyle
1/(\gamma_{i}+\mu_{i})$.

The flux of individuals for system (1) is shown in Figure 1.

$\displaystyle S_{1}(t)$$\displaystyle S_{2}(t)$$\displaystyle
E_{1}(t)$$\displaystyle E_{2}(t)$$\displaystyle I_{1}(t)$$\displaystyle
I_{2}(t)$$\displaystyle R_{1}(t)$$\displaystyle R_{2}(t)$$\displaystyle
D_{1}(t)$$\displaystyle
D_{2}(t)$$\displaystyle\beta_{11}I_{1}+\beta_{12}I_{2}$$\displaystyle\alpha_{1}$$\displaystyle\gamma_{1}$$\displaystyle\mu_{1}$$\displaystyle\beta_{21}I_{1}+\beta_{22}I_{2}$$\displaystyle\alpha_{2}$$\displaystyle\gamma_{2}$$\displaystyle\mu_{2}$$\displaystyle\beta_{21}$$\displaystyle\beta_{12}$$\displaystyle\beta_{11}$$\displaystyle\beta_{22}$
Figure 1. Flow diagram of individuals in system (3). The solid arrows
represents the flux of individuals and the dashed arrows represent the
interaction of infectious classes $\displaystyle I_{1}$ and $\displaystyle
I_{2}$.

Therefore the system (1) can be explicitly rewritten as

(3)
$\left\\{\begin{array}[]{lll}S_{1}^{\prime}=-(\beta_{11}I_{1}+\beta_{12}I_{2})S_{1}&,&S_{2}^{\prime}=-(\beta_{21}I_{1}+\beta_{22}I_{2})S_{2}\\\
E_{1}^{\prime}=(\beta_{11}I_{1}+\beta_{12}I_{2})S_{1}-\alpha_{1}E_{1}&,&E_{2}^{\prime}=(\beta_{21}I_{1}+\beta_{22}I_{2})S_{2}-\alpha_{2}E_{2}\\\
I_{1}^{\prime}=\alpha_{1}E_{1}-(\gamma_{1}+\mu_{1})I_{1}&,&I_{2}^{\prime}=\alpha_{2}E_{2}-(\gamma_{2}+\mu_{2})I_{2}\\\
R_{1}^{\prime}=\gamma_{1}I_{1}&,&R_{2}^{\prime}=\gamma_{2}I_{2}\\\
D_{1}^{\prime}=\mu_{1}I_{1}&,&D_{2}^{\prime}=\mu_{2}I_{2}\end{array}\right.$

The present work deals with two different problems concerning the two-group
SEIRD model. The first problem is estimating the final size of the epidemic.
The study of the final size SIR models with multi-group population was made in
[7] and [25], for example. We discuss it in Section 2 and prove that for both
groups the number of individuals who escape the epidemic is always positive.
We also present a way of calculating these limit values in Theorem 2.8.

For the second problem, consider two sets of parameters,
$\displaystyle\widetilde{\mathcal{A}}$,
$\displaystyle\widetilde{\mathcal{B}}$, $\displaystyle\widetilde{\Gamma}$,
$\displaystyle\widetilde{\mathcal{M}}$ and
$\displaystyle\mathcal{A},\mathcal{B}$, $\displaystyle\Gamma$,
$\displaystyle\mathcal{M}$, for the system (1) with the same initial
conditions, and let
$\displaystyle\widetilde{S},\widetilde{E},\widetilde{I},\widetilde{R},\widetilde{D}$
and $\displaystyle S,E,I,R,D$ be their respective sets of solutions. We would
like to obtain estimates for the distance between the solutions to second
group, i.e., for the quantities
$\displaystyle\|\widetilde{S}_{2}-S_{2}\|_{[0,T]}$,
$\displaystyle\|\widetilde{E}_{2}-E_{2}\|_{[0,T]}$,
$\displaystyle\|\widetilde{I}_{2}-I_{2}\|_{[0,T]}$,
$\displaystyle\|\widetilde{R}_{2}-R_{2}\|_{[0,T]}$ and
$\displaystyle\|\widetilde{D}_{2}-D_{2}\|_{[0,T]}$, based on the respective
distances between the solutions to first group. In other words, if the
respective solutions to the first group are close to each other, than we would
like to know how close the respective solutions to the second group are. We
can also think of this problem as follows: suppose that the parameters
$\displaystyle\widetilde{\mathcal{A}}$,
$\displaystyle\widetilde{\mathcal{B}}$, $\displaystyle\widetilde{\Gamma}$,
$\displaystyle\widetilde{\mathcal{M}}$ and the solutions to the second group
$\displaystyle\widetilde{S}_{2},\widetilde{E}_{2},\widetilde{I}_{2},\widetilde{R}_{2},\widetilde{D}_{2}$
are unknown, but the solutions to the first group are known. Also, suppose
that the parameters $\displaystyle\mathcal{A},\mathcal{B}$,
$\displaystyle\Gamma$, $\displaystyle\mathcal{M}$ and the respective solutions
for both groups are known. Then, from the error
$\displaystyle\|(\widetilde{S}_{1},\widetilde{E}_{1},\widetilde{I}_{1},\widetilde{R}_{1},\widetilde{D}_{1})-(S_{1},E_{1},I_{1},R_{1},D_{1})\|_{[0,T]}$,
we would like to estimate the respective error for the solutions to the second
group, i.e., we would like to specify a range around the solutions
$\displaystyle S_{2},E_{2},I_{2},R_{2},D_{2}$ where the unknown
$\displaystyle\widetilde{S}_{2},\widetilde{E}_{2},\widetilde{I}_{2},\widetilde{R}_{2},\widetilde{D}_{2}$
must be. Unfortunately, we were not able to provide such estimates in this
full generality. However, in Section 3, we present a way to solve this problem
under the additional hypothesis that the distance
$\displaystyle\|\widetilde{D}_{2}-D_{2}\|_{[0,T]}$ of number of deaths in
group 2 is known. Thus, in Theorem 3.4 we present these estimates under this
extra information on group 2 and with some hypotheses on the parameters of the
system.

Finally, in Section 4 we illustrate these results with the evolution of
COVID-19 epidemic in two two-group populations: New York County and its
neighboring counties in USA in Subsection 4.1, and the neighboring cities of
Petrolina and Juazeiro in Brazil, in Subsection 4.2. In the first example, we
are able to estimate the evolution of infectious cases in New York County only
using its data of deaths and the evolution of cases and deaths in the
neighboring counties. On the other hand, the same method applied to the second
example of the cities in Brazil, suggests that the cases in Juazeiro may be
under-reported, based on the evolution of cases and deaths in Petrolina and on
the deaths in Juazeiro.

Estimating the distance for the solutions in the second group as mentioned
above might be useful in situations where reporting on part of the cases of
infection is compromised. This could happen, for instance, to an unprivileged
minority or to poor people who might have no proper access to the health
system, or even a structural feature of the health system, where a relative
under-development of the health system in one group could lead to a lower
number of reported cases.

It is worth mentioning that we are assuming that infections by infected
individuals from group $\displaystyle j$ are proportional to $\displaystyle
I_{j}/N_{j}$, since the rate in which individuals in group $\displaystyle i$
become infected, as presented in (3), can be rewritten as
$\displaystyle\beta_{ij}I_{j}=\beta_{ij}N_{j}(I_{j}/N_{j})$. For many
infections that can cause deaths, this is not a realistic assumption since it
does not take into account the change in number of alive individuals. A more
accurate model should exclude the total number $\displaystyle D_{j}$ of
deceased individuals from the infecting process and thereby suppose that the
force of infection is proportional $\displaystyle I_{j}/(N_{j}-D_{j})$. For
many applications, this is not a significant change if the number of deceased
individuals is small relatively to the population size, as it is at the
beginning of an epidemic for example. But in a long-time analysis, these two
models may be considerably different, specially in cases where the mortality
rates $\displaystyle\mu_{i}$ are large. In the following section we present
results on the final size analysis not considering the effects that the change
in alive population size has on the infectious rates, and the results
presented does not seem adaptable to that case. Although this is a very
interesting problem, the authors are not aware of any final size analysis that
includes deaths apart from [14].

We are also assuming with model (3) that individuals from different groups
interact freely in a homogeneous space, as studied in [7] and [25]. The
probability of an infected individual from group $\displaystyle i$ entering
into contact with individuals from group $\displaystyle j$ is proportional to
the group size. For a much more general approach that takes into account
different population nodes and the flux of individuals between them, see for
example [10] and [15].

Finally, another limitation of model (1) that we would like to mention is that
it considers $\displaystyle 1/(\gamma_{i}+\mu_{i})$ as the average time of
infection for each group, whether it ends up with death or recovery. A more
precise model should take into account different periods of time for infection
depending on the severity and the outcome of the cases, since, as we mentioned
before, it may vary from 5 days to 6 weeks in the case of the COVID-19
disease. In [35], the authors present a more accurate model that considers
extra compartments that covers different age of infection depending on the
severity and also on the outcome of the disease.

##  2\. Final Size

In this section we analyze the final size of the epidemic for the two-group
model. We prove that the number of exposed and infectious individuals always
goes to zero as time goes to infinity, and that there is always a positive
number of individuals who escape the epidemic. We remark that our results for
the final size of the two-group SEIRD model and the respective proof are very
similar to the results given in [25] for the SIR model, since the models are
close. Therefore, our proof is adapted from the one given by the authors in
that reference. The main difference here is the inclusion of the class of
exposed individuals, which affects the time in which an individual starts to
contribute to the spread of the disease.

For the results in this section, we will need the following definition:

######  Definition 2.1.

The $\displaystyle n\times n$ matrix $\displaystyle A$ is irreducible when
$\displaystyle A$ cannot be transformed into block upper-triangular form by
simultaneous row-column permutations, i.e, when there is no permutation matrix
$\displaystyle P$ such that $\displaystyle PAP^{-1}$ is block upper-
triangular.

Note that, if $\displaystyle n=2$, $\displaystyle A$ is irreducible if and
only if $\displaystyle A_{12},A_{21}\neq 0$.

We will assume the following hypotheses on the matrices
$\displaystyle\mathcal{B}$ and $\displaystyle\Gamma+\mathcal{M}$:

######  Hypothesis 2.2.

$\displaystyle\mathcal{B}$ is a nonnegative irreducible matrix. By Definition
2.1, this is equivalent to assuming that
$\displaystyle\beta_{12},\beta_{21}>0$.

######  Hypothesis 2.3.

$\displaystyle\alpha_{1},\alpha_{2}>0$ ,
$\displaystyle\gamma_{1},\gamma_{2},\mu_{1},\mu_{2}\geq 0$ and
$\displaystyle\gamma_{1}+\mu_{1}>0$ and $\displaystyle\gamma_{2}+\mu_{2}>0$

We start with the following representation formulas for the solutions
$\displaystyle E_{i}$ and $\displaystyle I_{i}$ of the exposed and infectious
individuals, respectively, of each group:

######  Lemma 2.4.

For each $\displaystyle i=1,2$, let $\displaystyle S_{i}$, $\displaystyle
E_{i}$ and $\displaystyle I_{i}$ be the components of the solution to the
system (3). Then

(4) $\displaystyle\displaystyle E_{i}(t)$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle-
S_{i}(t)+\alpha_{i}e^{-\alpha_{i}t}\int_{0}^{t}S_{i}(s)e^{\alpha_{i}s}\,ds+(S_{i}(0)+E_{i}(0))e^{-\alpha_{i}t},$
(5) $\displaystyle\displaystyle I_{i}(t)$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle
I_{i}(0)e^{-(\gamma_{i}+\mu_{i})t}+\alpha_{i}e^{-(\gamma_{i}+\mu_{i})t}\int_{0}^{t}e^{(\gamma_{i}+\mu_{i})s}E_{i}(s)\,ds.$

###### Proof.

Firstly, adding the equations for $\displaystyle S$ and $\displaystyle E$, we
obtain

$\displaystyle S^{\prime}(t)+E^{\prime}(t)=-\mathcal{A}E(t).$

Componentwise, it holds for each $\displaystyle i=1,2$ that

(6) $E_{i}^{\prime}(t)+\alpha_{i}E_{i}(t)=-S^{\prime}_{i}(t).$

Thus, multiplying by $\displaystyle e^{\alpha_{i}t}$ we can rewrite it as

$\displaystyle\frac{d}{dt}\left(e^{\alpha_{i}t}E(t)\right)=-e^{\alpha_{i}t}S^{\prime}_{i}(t).$

Finally, an integration by parts leads to the result.

The expression for $\displaystyle I_{i}$ follows straight from the
$\displaystyle I$-equation. ∎

With Lemma 2.4 above, we can obtain the asymptotic values for the solutions
$\displaystyle E$ and $\displaystyle I$:

######  Lemma 2.5.

For any initial conditions, there exist the limits

$\displaystyle\displaystyle S^{\infty}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\lim_{t\to+\infty}S(t)\;,\;\;E^{\infty}=\lim_{t\to+\infty}E(t)\;,\;\;I^{\infty}=\lim_{t\to+\infty}I(t)\;,$
$\displaystyle\displaystyle R^{\infty}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\lim_{t\to+\infty}R(t)\;\;\mbox{ and
}\;\;D^{\infty}=\lim_{t\to+\infty}D(t).$

Furthermore, by Lemma 2.4, it follows that $\displaystyle E^{\infty}=0$ and
$\displaystyle I^{\infty}=0$.

###### Proof.

Since $\displaystyle\beta_{i,j}\geqslant 0$, $\displaystyle\gamma_{i}>0$ and
$\displaystyle\mu_{i}>0$ for all $\displaystyle i,j\in\\{1,2\\}$, it follows
from the system (1) that $\displaystyle S_{i}(t),R_{i}(t)$ and $\displaystyle
D_{i}(t)$ are monotone. Furthermore, by the restriction (2), it follows that
these functions are also bounded in $\displaystyle[0,\infty)$. Therefore,
there exists the limits $\displaystyle S_{i}^{\infty}$, $\displaystyle
R_{i}^{\infty}$ and $\displaystyle D_{i}^{\infty}$ for all $\displaystyle
i=1,2$.

For the functions $\displaystyle E_{i}$ and $\displaystyle I_{i}$,
$\displaystyle i=1,2$, we use the representation given by Lemma 2.4. Then
applying L’Hôpital’s Rule, we obtain that $\displaystyle
E_{i}^{\infty},I_{i}^{\infty}=0$ for $\displaystyle i=1,2$. ∎

######  Remark 2.6.

Another way to see that $\displaystyle E^{\infty}=0$ is estimating

$\displaystyle\displaystyle\alpha_{i}e^{-\alpha_{i}t}\int_{0}^{t}S_{i}(s)e^{\alpha_{i}s}\,ds$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\alpha_{i}e^{-\alpha_{i}t}\int_{0}^{t/2}S_{i}(s)e^{\alpha_{i}s}\,ds+\alpha_{i}e^{-\alpha_{i}t}\int_{t/2}^{t}S_{i}(s)e^{\alpha_{i}s}\,ds$
$\displaystyle\displaystyle\leq$ $\displaystyle\displaystyle
S_{i}(0)e^{-\alpha_{i}t}(e^{\alpha_{i}t/2}-1)+S_{i}(t/2)e^{-\alpha_{i}t}(e^{\alpha_{i}t}-e^{\alpha_{i}t/2})$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\left(S_{i}(0)-S_{i}(t/2)\right)e^{-\alpha_{i}t/2}+S_{i}(0)e^{-\alpha_{i}t}+S_{i}\left(t/2\right).$

Therefore

$\displaystyle 0\leq E_{i}(t)\leq
S_{i}(t/2)-S_{i}(t)+\left(S_{i}(0)-S_{i}(t/2)\right)e^{-\alpha_{i}t/2}+E_{i}(0)e^{-\alpha_{i}t}$

and the result follows. This also shows that $\displaystyle S(t)\to
S^{\infty}$ no faster than $\displaystyle E(t)\to 0$.

In the rest of this section we show that the asymptotic values $\displaystyle
S_{1}^{\infty}$ and $\displaystyle S_{2}^{\infty}$ are always positive in the
two-group model and we will obtain a formula to estimate these values.

Let us begin with the equation for $\displaystyle S(t)$ in (1). Since
$\displaystyle S_{1}(0),S_{2}(0)>0$, by continuity we can take
$\displaystyle\ln S(t)$ at least for $\displaystyle t$ sufficiently small.
Thus, we have that

$\displaystyle\frac{d}{dt}(\ln S(t))=-\mathcal{B}I(t)$

where the logarithm is taken componentwise. By adding the first three
equations we obtain

$\displaystyle\frac{d}{dt}(S(t)+E(t)+I(t))=-(\Gamma+\mathcal{M})I(t).$

Therefore,

$\displaystyle\frac{d}{dt}\Big{(}\mathcal{B}(\Gamma+\mathcal{M})^{-1}(S(t)+E(t)+I(t))-\ln
S(t)\Big{)}=0.$

Thus, we obtain that the function $\displaystyle\mathcal{F}$ defined by

(7) $\mathcal{F}(t):=\mathcal{B}(\Gamma+\mathcal{M})^{-1}(S(t)+E(t)+I(t))-\ln
S(t)$

must be constant equal to $\displaystyle\mathcal{F}(0)$. By continuity, the
logarithm above can be taken for all values of $\displaystyle t$, and the fact
that $\displaystyle S^{\infty}$ is finite implies that we can take the limit
$\displaystyle t\to+\infty$ and conclude that $\displaystyle
S_{1}^{\infty},S_{2}^{\infty}>0$. Taking the limit $\displaystyle t\to\infty$
we also have the relation

(8) $\mathcal{B}(\Gamma+\mathcal{M})^{-1}S^{\infty}-\ln
S^{\infty}=\mathcal{B}(\Gamma+\mathcal{M})^{-1}(S_{0}+E_{0}+I_{0})-\ln S_{0},$

since $\displaystyle I^{\infty}=0=E^{\infty}$. Denoting $\displaystyle
X_{0}:=S_{0}+E_{0}+I_{0}$, we can write

$\displaystyle\displaystyle\mathcal{B}(\Gamma+\mathcal{M})^{-1}S^{\infty}-\ln
S^{\infty}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathcal{B}(\Gamma+\mathcal{M})^{-1}X_{0}-\ln
S_{0}$ $\displaystyle\displaystyle\ln S^{\infty}$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle\ln
S_{0}+\mathcal{B}(\Gamma+\mathcal{M})^{-1}(S^{\infty}-X_{0})$ (9)
$\displaystyle\displaystyle S^{\infty}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle
S_{0}\exp\Big{(}\mathcal{B}(\Gamma+\mathcal{M})^{-1}(S^{\infty}-X_{0})\Big{)},$

where the exponential is taken componentwise. Equation (9) shows that
$\displaystyle S^{\infty}$ is a fixed point of the map $\displaystyle
T:\mathbb{R}^{2}\to\mathbb{R}^{2}$ defined by

(10)
$T(X)=S_{0}\exp\Big{(}\mathcal{B}(\Gamma+\mathcal{M})^{-1}(X-X_{0})\Big{)}$

or, denoting $\displaystyle X=(x,y)$ and $\displaystyle
T(X)=(T_{1}(x,y),T_{2}(x,y))$,

(11) $\displaystyle\displaystyle T_{1}(x,y)$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle
S_{1}(0)\exp\left(\frac{\beta_{11}}{\gamma_{1}+\mu_{1}}(x-x_{0})+\frac{\beta_{12}}{\gamma_{2}+\mu_{2}}(y-y_{0})\right),$
(12) $\displaystyle\displaystyle T_{2}(x,y)$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle
S_{2}(0)\exp\left(\frac{\beta_{21}}{\gamma_{1}+\mu_{1}}(x-x_{0})+\frac{\beta_{22}}{\gamma_{2}+\mu_{2}}(y-y_{0})\right),$

with $\displaystyle x_{0}=S_{1}(0)+E_{1}(0)+I_{1}(0)$ and $\displaystyle
y_{0}=S_{2}(0)+E_{2}(0)+I_{2}(0)$. Before we continue with the analysis of the
epidemic final size, let us define the following notation for partial ordering
of vectors in $\displaystyle\mathbb{R}^{2}$: given $\displaystyle
X=(X_{1},X_{2}),Y=(Y_{1},Y_{2})\in\mathbb{R}^{2}$, se say that

$\displaystyle\begin{matrix}[l]X\leq Y\Longleftrightarrow X_{i}\leq
Y_{i}&\mbox{ for all }&i=1,2\,,\\\ X<Y\Longleftrightarrow X\leq Y\mbox{ and
}X_{i}<Y_{i}&\mbox{ for some }&i=1,2\,,\\\ X\ll Y\Longleftrightarrow
X_{i}<Y_{i}&\mbox{ for all }&i=1,2\,.\end{matrix}$

######  Theorem 2.7.

Let $\displaystyle T:\mathbb{R}^{2}\to\mathbb{R}^{2}$ be the map defined in
(10) with $\displaystyle X_{0}=S_{0}+I_{0}+E_{0}$. Then, under the hypotheses
2.2 and 2.3, we have

  1. a)

$\displaystyle T$ is componentwise increasing

  2. b)

$\displaystyle T(S_{0})=S_{0}\Leftrightarrow E_{0}=0=I_{0}$

  3. c)

$\displaystyle E_{0}+I_{0}>0\Leftrightarrow 0\ll T(0)\ll T(S_{0})\ll S_{0}$

  4. d)

The derivative of $\displaystyle T$ is componentwise increasing and given by

(13)
$DT(x,y)=\begin{pmatrix}\frac{\beta_{11}}{\gamma_{1}+\mu_{1}}T_{1}(x,y)&\frac{\beta_{12}}{\gamma_{2}+\mu_{2}}T_{1}(x,y)\\\
\frac{\beta_{21}}{\gamma_{1}+\mu_{1}}T_{2}(x,y)&\frac{\beta_{22}}{\gamma_{2}+\mu_{2}}T_{2}(x,y)\end{pmatrix}.$

###### Proof.

The point a) follows from hypotheses 2.2 and 2.3.

For point b), note that $\displaystyle T(S_{0})=S_{0}$ is equivalent to

$\displaystyle\frac{\beta_{i1}}{\gamma_{1}+\mu_{1}}(E_{1}(0)+I_{1}(0))+\frac{\beta_{i2}}{\gamma_{2}+\mu_{2}}(E_{2}(0)+I_{2}(0))=0\,\qquad
i=1,2.$

Again by hypotheses 2.2 and 2.3, this is true only when $\displaystyle
E(0)=0=I(0)$, since these functions are nonnegative.

Point c) follows from points a) and b), and point d) follows from the
definition of $\displaystyle T$.

∎

It follows from the points (a) and (c) in the theorem above that, for every
$\displaystyle n\in\mathbb{N}$

$\displaystyle 0\ll T(0)\ll\cdots\ll T^{n}(0)\ll\cdots\ll
T^{n}(S_{0})\ll\cdots\ll T(S_{0})\ll S_{0}.$

Therefore, there exist the limits

$\displaystyle
S^{-}=\lim_{n\to\infty}T^{n}(0),\;\;\;S^{+}=\lim_{n\to\infty}T^{n}(S_{0}),$

they satisfy $\displaystyle S^{-}\leq S^{+}$ and are fixed points of
$\displaystyle T$ on $\displaystyle[0,S_{1}(0)]\times[0,S_{2}(0)]$, since
$\displaystyle T$ is continuous. Note that, due to the monotonicity of
$\displaystyle T$, there is no other fixed points of $\displaystyle T$ in
$\displaystyle[0,S_{1}(0)]\times[0,S_{2}(0)]\backslash[S_{1}^{-},S_{1}^{+}]\times[S_{2}^{-},S_{2}^{+}]$.
Thus, the limit value $\displaystyle S^{\infty}$ must satisfy $\displaystyle
S^{-}\leq S^{\infty}\leq S^{+}$.

######  Theorem 2.8 (Final size).

Let us assume that (2.2) and (2.3) hold. If $\displaystyle S_{0}\gg 0$ and
$\displaystyle E_{0}+I_{0}>0$, then the final state of the epidemic model is

  1. i)

$\displaystyle 0\neq S^{\infty}=\lim_{n\to\infty}T^{n}(0)$

  2. ii)

$\displaystyle I^{\infty}=0=E^{\infty}$

  3. iii)

$\displaystyle R^{\infty}+D^{\infty}=N-S^{\infty}$

  4. iv)

$\displaystyle
R^{\infty}=\Gamma(\Gamma+\mathcal{M})^{-1}(N-S^{\infty}-R_{0}-D_{0})+R_{0}$

  5. v)

$\displaystyle
D^{\infty}=\mathcal{M}(\Gamma+\mathcal{M})^{-1}(N-S^{\infty}-R_{0}-D_{0})+D_{0}$

###### Proof.

For item i), let us show that under Hypotheses 2.2 and 2.3, we have
$\displaystyle S^{-}=S^{+}$, i.e., the map $\displaystyle T$ has only one
fixed point.

Let us suppose that $\displaystyle S^{-}<S^{+}$. Then, by the Taylor’s Theorem
we can write

$\displaystyle\displaystyle
S^{+}-S^{-}=T(S^{+})-T(S^{-})=\int_{0}^{1}DT(S^{-}+\tau(S^{+}-S^{-}))(S^{+}-S^{-})\;d\tau.$

From (13), we see that $\displaystyle DT$ is componentwise increasing by the
monotonicity of $\displaystyle T$. It follows that $\displaystyle
DT(S^{-}+\tau(S^{+}-S^{-}))\leq DT(S^{+})$ componentwise for all
$\displaystyle\tau\in[0,1]$. Therefore, we obtain

(14) $\displaystyle\displaystyle S^{+}-S^{-}\leq DT(S^{+})(S^{+}-S^{-}).$

Since $\displaystyle
DT(S^{+})=\mathrm{diag}(S^{+})\mathcal{B}(\Gamma+\mathcal{M})^{-1}$, we have
from Hypothesis 2.2 and 2.3 that $\displaystyle DT(S^{+})$ is nonnegative
irreducible. Hence, by the Perron-Frobenius Theorem, let
$\displaystyle\lambda>0$ be the dominant eigenvalue and $\displaystyle W\gg 0$
an associated left eigenvector. From (14) we have

(15) $\displaystyle\displaystyle W^{T}(S^{+}-S^{-})\leq
W^{T}DT(S^{+})(S^{+}-S^{-})=\lambda W^{T}(S^{+}-S^{-}),$

and thus, $\displaystyle\lambda\geqslant 1$. By the definition of
$\displaystyle S^{+}$, we have $\displaystyle S_{0}\geqslant S^{+}$ and then
$\displaystyle DT(S^{+}+\tau(S_{0}-S^{+}))\geqslant DT(S^{+})$ componentwise
for all $\displaystyle\tau\geqslant 0$. Therefore,

$\displaystyle\displaystyle T(S_{0})-S^{+}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle
T(S_{0})-T(S^{+})=\int_{0}^{1}DT(S^{+}+\tau(S_{0}-S^{+}))(S_{0}-S^{+})\;d\tau$
$\displaystyle\displaystyle\geqslant$ $\displaystyle\displaystyle
DT(S^{+})(S_{0}-S^{+}).$

Multiplying the inequality above by $\displaystyle W^{T}$,

$\displaystyle\displaystyle W^{T}(T(S_{0})-S^{+})$
$\displaystyle\displaystyle\geqslant$ $\displaystyle\displaystyle
W^{T}DT(S^{+}))(S_{0}-S^{+})=\lambda W^{T}(S_{0}-S^{+})$
$\displaystyle\displaystyle\geqslant$ $\displaystyle\displaystyle
W^{T}(S_{0}-S^{+})$

and, since $\displaystyle W\gg 0$ this implies that $\displaystyle
T(S_{0})\geqslant S_{0}$, which cannot happen if $\displaystyle E_{0}+I_{0}>0$
by item (c) of Theorem 2.7.

The item ii) was already proved.

Item iii) follows from item ii) and (2).

For the item iv), we can use the equations for $\displaystyle R$ and
$\displaystyle D$ in (1) and obtain

$\displaystyle
R^{\prime}+D^{\prime}=(\Gamma+\mathcal{M})I=(\Gamma+\mathcal{M})\Gamma^{-1}\Gamma
I=(\Gamma+\mathcal{M})\Gamma^{-1}R^{\prime},$

which implies that

$\displaystyle
R(t)+D(t)-R_{0}-D_{0}=(\Gamma+\mathcal{M})\Gamma^{-1}(R(t)-R_{0}),\qquad\forall\;t>0.$

Therefore, as $\displaystyle t\to+\infty$ we conclude that

$\displaystyle
R^{\infty}+D^{\infty}-R_{0}-D_{0}=(\Gamma+\mathcal{M})\Gamma^{-1}(R^{\infty}-R_{0}).$

Using item iii), we can write

$\displaystyle
N-S^{\infty}-R_{0}-D_{0}=(\Gamma+\mathcal{M})\Gamma^{-1}(R^{\infty}-R_{0})$

and thus

$\displaystyle
R^{\infty}=\Gamma(\Gamma+\mathcal{M})^{-1}(N-S^{\infty}-R_{0}-D_{0})+R_{0}.$

Finally, item v) follows from items iii) and iv). ∎

######  Remark 2.9.

It follows from items c) and d) that the final size of recovered and deceased
classes are $\displaystyle
R^{\infty}_{i}=\frac{\gamma_{i}}{\gamma_{i}+\mu_{i}}(N_{i}-S^{\infty}_{i}-R_{i}(0)-D_{i}(0))+R_{i}(0)$
and $\displaystyle
D^{\infty}_{i}=\frac{\mu_{i}}{\gamma_{i}+\mu_{i}}(N_{i}-S^{\infty}_{i}-R_{i}(0)-D_{i}(0))+D_{i}(0)$
for each group $\displaystyle i$.

######  Remark 2.10.

It follows from item i) that the final size $\displaystyle S^{\infty}$ can
also be estimated by using only the map $\displaystyle T$, instead of solving
numerically the system (1).

##  3\. Distance estimates for the second group

In this section we study the following problem: given two sets of parameters,
$\displaystyle\widetilde{\mathcal{A}}$,
$\displaystyle\widetilde{\mathcal{B}}$, $\displaystyle\widetilde{\Gamma}$,
$\displaystyle\widetilde{\mathcal{M}}$ and
$\displaystyle\mathcal{A},\mathcal{B}$, $\displaystyle\Gamma$,
$\displaystyle\mathcal{M}$, for the system (1), let
$\displaystyle\widetilde{S},\widetilde{E},\widetilde{I},\widetilde{R},\widetilde{D}$
and $\displaystyle S,E,I,R,D$ be their respective set of solutions, we want to
estimate the distance related to the solutions to second group
$\displaystyle\|(\widetilde{S}_{2},\widetilde{E}_{2},\widetilde{I}_{2},\widetilde{R}_{2})-(S_{2},E_{2},I_{2},R_{2})\|_{[0,T]}$
based on the the distance
$\displaystyle\|(\widetilde{S}_{1},\widetilde{E}_{1},\widetilde{I}_{1},\widetilde{R}_{1},\widetilde{D}_{1},\widetilde{D}_{2})-(S_{1},E_{1},I_{1},R_{1},D_{1},D_{2})\|_{[0,T]}$,
i.e., based on the distances between solutions of the first group plus on the
distance between the evolution of deaths in group 2 too.

We start with the following lemma on the ratio of susceptible class:

######  Lemma 3.1.

The solutions $\displaystyle S_{1}(t)$ and $\displaystyle S_{2}(t)$ of
susceptible individuals satisfy

$\displaystyle\displaystyle\frac{S_{1}(t)}{S_{2}(t)}=\frac{S_{1}(0)}{S_{2}(0)}\exp\left\\{\int_{0}^{t}\left(\left(\beta_{21}-\beta_{11}\right)I_{1}(s)+\left(\beta_{22}-\beta_{12}\right)I_{2}(s)\right)\,ds\right\\}.$

###### Proof.

From the equations for $\displaystyle S_{1}$ and $\displaystyle S_{2}$ in (3),
we have

$\displaystyle\displaystyle\displaystyle\left[\frac{S_{1}}{S_{2}}\right]^{\prime}$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\frac{(\beta_{21}I_{1}+\beta_{22}I_{2})S_{1}S_{2}-(\beta_{11}I_{1}+\beta_{12}I_{2})S_{1}S_{2}}{S_{2}^{2}}$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\big{(}\left(\beta_{21}-\beta_{11}\right)I_{1}+\left(\beta_{22}-\beta_{12}\right)I_{2}\big{)}\frac{S_{1}}{S_{2}}.$

and thus, the result follows. ∎

The following two elementary results will be useful for the estimates in
Theorem 3.4, therefore we include them with proofs.

######  Remark 3.2.

Let $\displaystyle a,b>0$ and $\displaystyle f:[0,+\infty)\to\mathbb{R}$
defined by

$\displaystyle f(\varepsilon)=(1-b\varepsilon)e^{-a\varepsilon}.$

Then, the Mean Value Theorem implies that, for every
$\displaystyle\varepsilon>0$, there holds

$\displaystyle\left|1-(1-b\varepsilon)e^{-a\varepsilon}\right|\leq\max_{\xi\in[0,\varepsilon]}|f^{\prime}(\xi)|\varepsilon=|f^{\prime}(0)|\varepsilon=(a+b)\varepsilon.$

######  Lemma 3.3.

Let $\displaystyle a,b>0$ and $\displaystyle g:[0,T]\to\mathbb{R}$ defined by

$\displaystyle g(s)=ae^{a(s-T)}-be^{b(s-T)}.$

Then there exists at most one $\displaystyle\tau\in[0,T)$ such that
$\displaystyle g(\tau)=0$. Also, if $\displaystyle h:[0,T]\to\mathbb{R}_{+}$
is nonincreasing, then

(16)
$\int_{0}^{T}g(s)h(s)\,ds\leq\left\\{\begin{matrix}h(0)(e^{-bT}-e^{-aT}),\hfill&\mbox{if
}0<b<a\\\
h(0)(e^{-bT}-e^{-aT})+(h(0)-h(T))\ln\left(\frac{b}{a}\right),&\mbox{if
}0<a<b\end{matrix}\right..$

###### Proof.

To see this, note that $\displaystyle g(\tau)=0$ if and only if
$\displaystyle\tau=T-\frac{1}{a-b}\ln\left(\frac{a}{b}\right)\in[0,T)$, so it
changes sign only once. If $\displaystyle 0<b<a$, we have $\displaystyle
g(T)=a-b>0$ and then $\displaystyle g$ nonnegative on $\displaystyle[0,T]$ or
it changes sign at $\displaystyle s=\tau$. In this case, for any
$\displaystyle h$ nonincreasing nonnegative function we can write

$\displaystyle\int_{0}^{T}\\!g(s)h(s)ds=\int_{0}^{\tau}\\!g(s)h(s)\,ds+\int_{\tau}^{T}\\!g(s)h(s)\,ds\leq
h(\tau)\int_{0}^{\tau}\\!g(s)\,ds+h(\tau)\int_{\tau}^{T}\\!g(s)\,ds$

and the estimate follows. In the case where $\displaystyle g$ is nonnegative,
the estimate is immediate.

If $\displaystyle 0<a<b$, we have $\displaystyle g(T)=a-b<0$ and then

$\displaystyle\displaystyle\int_{0}^{T}g(s)h(s)\,ds$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\int_{0}^{\tau}g(s)h(s)\,ds+\int_{\tau}^{T}g(s)h(s)\,ds$
$\displaystyle\displaystyle\leq$ $\displaystyle\displaystyle
h(0)\left(e^{a(\tau-T)}-e^{-aT}-e^{b(\tau-T)}+e^{-bT}\right)\\!+\\!h(T)\left(e^{b(\tau-T)}\\!-\\!e^{a(\tau-T)}\right)$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle
h(0)\left(e^{-bT}-e^{-aT}\right)+(h(0)-h(T))\left(e^{a(\tau-T)}-e^{b(\tau-T)}\right).$

By the Mean Value Theorem, there exists a $\displaystyle a<\xi<b$ satisfying

$\displaystyle\displaystyle e^{a(\tau-T)}-e^{b(\tau-T)}$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle(T-\tau)(b-a)e^{\xi(\tau-T)}\leq(T-\tau)(b-a)e^{a(\tau-T)}$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\ln\left(\frac{b}{a}\right)\left(\frac{a}{b}\right)^{\frac{a}{b-a}}\leq\ln\left(\frac{b}{a}\right),$

where for the last equality we used that $\displaystyle
T-\tau=\frac{1}{a-b}\ln\left(\frac{a}{b}\right)$. ∎

For the main result of this section, we will denote the norm $\displaystyle
L^{\infty}(0,T)$ by $\displaystyle\|\,.\,\|_{[0,T]}$.

######  Theorem 3.4.

Let
$\displaystyle(\widetilde{S},\widetilde{E},\widetilde{I},\widetilde{R},\widetilde{D})$
and $\displaystyle(S,E,I,R,D)$ be two solutions for (1) corresponding to two
sets of parameters $\displaystyle\widetilde{\mathcal{A}}$,
$\displaystyle\widetilde{\mathcal{B}}$, $\displaystyle\widetilde{\Gamma}$,
$\displaystyle\widetilde{\mathcal{M}}$ and
$\displaystyle\mathcal{A},\mathcal{B}$, $\displaystyle\Gamma$,
$\displaystyle\mathcal{M}$ respectively, with the same initial conditions. If

(17)
$\|(\widetilde{S}_{1},\widetilde{E}_{1},\widetilde{I}_{1},\widetilde{R}_{1},\widetilde{D}_{1})-(S_{1},E_{1},I_{1},R_{1},D_{1})\|_{[0,T]}\leq\varepsilon$

and

  1. i)

$\displaystyle\frac{\widetilde{\beta}_{2i}-\widetilde{\beta}_{1i}}{\widetilde{\mu}_{i}}=\frac{\beta_{2i}-\beta_{1i}}{\mu_{i}},\mbox{
for }i=1,2;$

  2. ii)

$\displaystyle\frac{\widetilde{\alpha}_{1}}{\widetilde{\alpha}_{2}}=\frac{\alpha_{1}}{\alpha_{2}};$

  3. iii)

$\displaystyle\widetilde{\gamma}_{1}-\widetilde{\gamma}_{2}=\gamma_{1}-\gamma_{2};$

  4. iv)

$\displaystyle\widetilde{\mu}_{1}-\widetilde{\mu}_{2}=\mu_{1}-\mu_{2};$

  5. v)

$\displaystyle\|D_{2}-\widetilde{D}_{2}\|_{[0,T]}<\varepsilon$

then

$\displaystyle\|(\widetilde{S}_{2},\widetilde{E}_{2},\widetilde{I}_{2},\widetilde{R}_{2})-(S_{2},E_{2},I_{2},R_{2})\|_{[0,T]}<K(\varepsilon)\varepsilon,$

where $\displaystyle K(\varepsilon)$ is bounded as
$\displaystyle\varepsilon\to 0$. Moreover, if
$\displaystyle\widetilde{\beta}_{22}-\widetilde{\beta}_{12}=\beta_{22}-\beta_{12}$,
then the hypothesis (v) is not necessary.

###### Proof.

For simplicity let us define $\displaystyle
V_{i}=\frac{\widetilde{\beta}_{2i}-\widetilde{\beta}_{1i}}{\widetilde{\mu}_{i}}=\frac{\beta_{2i}-\beta_{1i}}{\mu_{i}},\mbox{
for }i=1,2$.

Step 1: For the susceptible part, by Lemma 3.1 and the equation for
$\displaystyle D_{1}$ and $\displaystyle D_{2}$ we have

$\displaystyle\displaystyle\frac{\widetilde{S}_{2}}{S_{2}}=\displaystyle\frac{\widetilde{S}_{1}}{S_{1}}\exp\\{V_{1}(D_{1}(t)-\widetilde{D}_{1}(t))+V_{2}(D_{2}(t)-\widetilde{D}_{2}(t))\\}.$

Firstly, let us assume that for a given $\displaystyle t>0$ we have
$\displaystyle\widetilde{S}_{2}(t)\leq S_{2}(t)$. Since
$\displaystyle|S_{1}(t)-\widetilde{S}_{1}(t)|<\varepsilon$ by hypothesis, we
obtain

(18) $\displaystyle\displaystyle 0$ $\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle
S_{2}(t)-\widetilde{S}_{2}(t)\leq{S}_{2}(t)\left(1-\frac{\widetilde{S}_{2}(t)}{S_{2}(t)}\right)$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle{S}_{2}(t)\left(1-\frac{\widetilde{S}_{1}(t)}{S_{1}(t)}e^{V_{1}(D_{1}(t)-\widetilde{D}_{1}(t))+V_{2}(D_{2}(t)-\widetilde{D}_{2}(t))}\right)$
$\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle{S}_{2}(t)\left(1-\left(1-\frac{\varepsilon}{S_{1}(t)}\right)e^{V_{1}(D_{1}(t)-\widetilde{D}_{1}(t))+V_{2}(D_{2}(t)-\widetilde{D}_{2}(t))}\right).$

Since $\displaystyle|D_{i}(t)-\widetilde{D}_{i}(t)|\leq\varepsilon$, we can
write

(19) $\displaystyle\displaystyle 0$ $\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle
S_{2}(t)-\widetilde{S}_{2}(t)\leq{S}_{2}(t)\left(1-\left(1-\frac{\varepsilon}{S_{1}(t)}\right)e^{-(|V_{1}|+|V_{2}|)\varepsilon}\right)$
$\displaystyle\displaystyle\leq$ $\displaystyle\displaystyle
S_{2}(t)\left(|V_{1}|+|V_{2}|+\frac{1}{S_{1}(t)}\right)\varepsilon,$

where, for the last inequality, we used the Remark 3.2 for $\displaystyle
a=|V_{1}|+|V_{2}|$ and $\displaystyle b=1/S_{1}(t)$.

In the case where $\displaystyle S_{2}(t)\leq\widetilde{S}_{2}(t)$ for a given
$\displaystyle t>0$, we proceed analogously and obtain

(20) $\displaystyle\displaystyle 0\leq\widetilde{S}_{2}(t)-S_{2}(t)$
$\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle\widetilde{S}_{2}(t)\left(1-\frac{{S}_{1}(t)}{\widetilde{S}_{1}(t)}e^{-V_{1}(D_{1}(t)-\widetilde{D}_{1}(t))-V_{2}(D_{2}(t)-\widetilde{D}_{2}(t))}\right)$
$\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle\widetilde{S}_{2}(t)\left(1-\frac{{S}_{1}(t)}{{S}_{1}(t)+\varepsilon}e^{-(|V_{1}|+|V_{2}|)\varepsilon}\right)$
$\displaystyle\displaystyle\leq$ $\displaystyle\displaystyle
S_{2}(0)\left(1-\left(1-\frac{\varepsilon}{{S}_{1}(t)+\varepsilon}\right)e^{-(|V_{1}|+|V_{2}|)\varepsilon}\right),$

since the initial condition is the same for $\displaystyle S_{2}$ and
$\displaystyle\widetilde{S}_{2}$. Thus, for every fixed $\displaystyle T>0$,
we can apply (20) and (19) to each $\displaystyle t\in[0,T]$ and obtain that

(21) $\|\widetilde{S}_{2}-S_{2}\|_{[0,T]}\leq K_{S}\varepsilon,$

where $\displaystyle
K_{S}=S_{2}(0)\left(|V_{1}|+|V_{2}|+\frac{1}{S_{1}(T)}\right).$

Step 2: Suppose that $\displaystyle E_{2}(t)\geq\widetilde{E}_{2}(t)$. Note
that

$\displaystyle
S_{2}(t)=S_{2}(t)e^{-\alpha_{2}t}+\alpha_{2}\int_{0}^{t}S(t)e^{\alpha_{2}(s-t)}$

and then, by Lemma 3.1 and Lemma 2.4 we have

$\displaystyle\displaystyle E_{2}(t)-\widetilde{E}_{2}(t)=$
$\displaystyle\displaystyle-
S_{2}(t)+\alpha_{2}\int_{0}^{t}e^{\alpha_{2}(s-t)}S_{2}(s)ds+(S_{2}(0)+E_{2}(0))e^{-\alpha_{2}t}$
$\displaystyle\displaystyle+\widetilde{S}_{2}(t)-\widetilde{\alpha}_{2}\int_{0}^{t}e^{\widetilde{\alpha}_{2}(s-t)}\widetilde{S}_{2}(s)ds-(S_{2}(0)+E_{2}(0))e^{-\widetilde{\alpha}_{2}t}$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\alpha_{2}\int_{0}^{t}e^{\alpha_{2}(s-t)}(S_{2}(s)-S_{2}(t))ds+(S_{2}(0)+E_{2}(0)-S_{2}(t))e^{-\alpha_{2}t}$
$\displaystyle\displaystyle-\widetilde{\alpha}_{2}\int_{0}^{t}e^{\widetilde{\alpha}_{2}(s-t)}(\widetilde{S}_{2}(s)-\widetilde{S}_{2}(t))ds-(S_{2}(0)+E_{2}(0)-\widetilde{S}_{2}(t))e^{-\widetilde{\alpha}_{2}t}$
$\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle\alpha_{2}\int_{0}^{t}e^{\alpha_{2}(s-t)}(S_{2}(s)-S_{2}(t))ds+(S_{2}(0)+E_{2}(0)-S_{2}(t))e^{-\alpha_{2}t}$
$\displaystyle\displaystyle-\widetilde{\alpha}_{2}\int_{0}^{t}e^{\widetilde{\alpha}_{2}(s-t)}(S_{2}(s)-S_{2}(t)-2K_{S}\varepsilon)ds$
$\displaystyle\displaystyle-(S_{2}(0)+E_{2}(0)-S_{2}(t)-K_{S}\varepsilon)e^{-\widetilde{\alpha}_{2}t}$
$\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle\int_{0}^{t}\left(\alpha_{2}e^{\alpha_{2}(s-t)}-\widetilde{\alpha}_{2}e^{\widetilde{\alpha}_{2}(s-t)}\right)(S_{2}(s)-S_{2}(t))ds+K_{S}\varepsilon
e^{-\widetilde{\alpha}_{2}t}$
$\displaystyle\displaystyle+(S_{2}(0)+E_{2}(0)-S_{2}(t))\left(e^{-\alpha_{2}t}-e^{-\widetilde{\alpha}_{2}t}\right)+2\widetilde{\alpha}_{2}K_{S}\varepsilon(1-e^{-\widetilde{\alpha}_{2}t}).$

Now, supposing that $\displaystyle\widetilde{\alpha}_{2}<\alpha_{2}$, by the
Lemma 3.3 for $\displaystyle a=\alpha_{2}$, $\displaystyle
b=\widetilde{\alpha}_{2}$ and $\displaystyle h(s)=S_{2}(s)-S_{2}(t)$, we can
estimate

$\displaystyle\displaystyle E_{2}(t)\\!-\\!\widetilde{E}_{2}(t)$
$\displaystyle\displaystyle\\!\leq\\!$
$\displaystyle\displaystyle(S_{2}(0)\\!-\\!S_{2}(t))\\!\left(e^{-\widetilde{\alpha}_{2}t}\\!-\\!e^{-\alpha_{2}t}\right)\\!+\\!(S_{2}(0)\\!+\\!E_{2}(0)\\!-\\!S_{2}(t))\\!\left(e^{-\alpha_{2}t}\\!-\\!e^{-\widetilde{\alpha}_{2}t}\right)$
$\displaystyle\displaystyle\;\;\;\;+2\widetilde{\alpha}_{2}K_{S}\varepsilon+(1-2\widetilde{\alpha}_{2})K_{S}\varepsilon
e^{-\widetilde{\alpha}_{2}t}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle-
E_{2}(0)\left(e^{-\widetilde{\alpha}_{2}t}-e^{-\alpha_{2}t}\right)+\left(2\widetilde{\alpha}_{2}+(1-2\widetilde{\alpha}_{2})e^{-\widetilde{\alpha}_{2}t}\right)K_{S}\varepsilon$
$\displaystyle\displaystyle\leq$ $\displaystyle\displaystyle-
E_{2}(0)(\alpha_{2}-\widetilde{\alpha}_{2})te^{-\alpha_{2}t}+K_{S}\varepsilon$

and, if $\displaystyle\widetilde{\alpha}_{2}>\alpha_{2}$ the Lemma 3.3 implies
that

$\displaystyle\displaystyle E_{2}(t)-\widetilde{E}_{2}(t)$
$\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle(S_{2}(0)-S_{2}(t))\left(e^{-\widetilde{\alpha}_{2}t}-e^{-\alpha_{2}t}\right)+(S_{2}(0)-S_{2}(t))\ln\left(\frac{\widetilde{\alpha}_{2}}{{\alpha}_{2}}\right)$
$\displaystyle\displaystyle+(S_{2}(0)+E_{2}(0)-S_{2}(t))\left(e^{-\alpha_{2}t}-e^{-\widetilde{\alpha}_{2}t}\right)+K_{S}\varepsilon$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle-
E_{2}(0)\left(e^{-\widetilde{\alpha}_{2}t}-e^{-\alpha_{2}t}\right)+(S_{2}(0)-S_{2}(t))\ln\left(\frac{\widetilde{\alpha}_{2}}{{\alpha}_{2}}\right)+K_{S}\varepsilon$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle
E_{2}(0)(\widetilde{\alpha}_{2}-\alpha_{2})te^{-\alpha_{2}t}+(S_{2}(0)-S_{2}(t))\ln\left(\frac{\widetilde{\alpha}_{2}}{{\alpha}_{2}}\right)+K_{S}\varepsilon.$

It remains to show that the parameters $\displaystyle\widetilde{\alpha}_{2}$
and $\displaystyle{\alpha}_{2}$ are near each other under the hypothesis of
the Theorem. To see this, we can integrate (6) for $\displaystyle i=1$ and
obtain

(22) $\alpha_{1}\int_{0}^{t}E_{1}(s)\,ds=S_{1}(0)-S_{1}(t)+E_{1}(0)-E_{1}(t).$

Thus, using the same argument for the parameter
$\displaystyle\widetilde{\alpha}_{1}$ and the hypothesis (ii) and supposing
without loss of generality that
$\displaystyle\widetilde{\alpha}_{2}>\alpha_{2}$, we have for every
$\displaystyle t>0$ that

$\displaystyle\displaystyle 1\leq\frac{\widetilde{\alpha}_{2}}{\alpha_{2}}$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\frac{\widetilde{\alpha}_{1}}{\alpha_{1}}=\left(\frac{S_{1}(0)-\widetilde{S}_{1}(t)+E_{1}(0)-\widetilde{E}_{1}(t)}{S_{1}(0)-S_{1}(t)+E_{1}(0)-E_{1}(t)}\right)\left(\frac{\int_{0}^{t}E_{1}ds}{\int_{0}^{t}\widetilde{E}_{1}ds}\right)$
$\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle\left(1+\frac{2\varepsilon}{S_{1}(0)-S_{1}(t)+E_{1}(0)-E_{1}(t)}\right)\left(1+\frac{\varepsilon
t}{\int_{0}^{t}\widetilde{E}_{1}ds}\right).$

Therefore, using for example $\displaystyle t=1$, there exists constants
$\displaystyle C_{1},C_{2}>0$ such that

$\displaystyle\displaystyle 1$
$\displaystyle\displaystyle\leq\frac{\widetilde{\alpha}_{2}}{\alpha_{2}}$
$\displaystyle\displaystyle\leq 1+C_{1}\varepsilon+C_{2}\varepsilon^{2}$
$\displaystyle\displaystyle\mbox{and thus}\hskip 102.0307pt0$
$\displaystyle\displaystyle<\widetilde{\alpha}_{2}-\alpha_{2}$
$\displaystyle\displaystyle\leq\alpha_{2}C_{1}\varepsilon+\alpha_{2}C_{2}\varepsilon^{2}\hskip
155.2697pt.$

Therefore, we conclude that, for
$\displaystyle\widetilde{\alpha}_{2}<\alpha_{2}$, there holds

$\displaystyle E_{2}(t)-\widetilde{E}_{2}(t)\leq-
E_{2}(0)(\alpha_{2}-\widetilde{\alpha}_{2})te^{-\alpha_{2}t}+K_{S}\varepsilon$

and, for $\displaystyle\widetilde{\alpha}_{2}>\alpha_{2}$ we have

$\displaystyle
E_{2}(t)-\widetilde{E}_{2}(t)\leq\left(\alpha_{2}E_{2}(0)te^{-\alpha_{2}t}+S_{2}(0)-S_{2}(t)\right)(C_{1}+C_{2}\varepsilon)\varepsilon+K_{S}\varepsilon.$

Therefore,

(23) $\|E_{2}-\widetilde{E}_{2}\|_{[0,T]}\leq K_{E}(\varepsilon)\varepsilon,$

where

$\displaystyle
K_{E}(\varepsilon)=\left(\frac{E_{2}(0)}{e}+S_{2}(0)-S_{2}(T)\right)(C_{1}+C_{2}\varepsilon)+K_{S},$

since for every $\displaystyle\lambda>0$, the function $\displaystyle
g(t):=te^{-\lambda t}$ satisfies $\displaystyle g(t)\leq\frac{1}{\lambda e}$.

Step 3: For the functions $\displaystyle I_{2}$ and
$\displaystyle\widetilde{I}_{2}(t)$, we can use (5) and write

(24) $\displaystyle\displaystyle I_{2}(t)-\widetilde{I}_{2}(t)$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle
I_{2}(0)\left(e^{-(\gamma_{2}+\mu_{2})t}-e^{-(\widetilde{\gamma}_{2}+\widetilde{\mu}_{2})t}\right)$
$\displaystyle\displaystyle+\int_{0}^{t}\left(\alpha_{2}e^{(\gamma_{2}+\mu_{2})(s-t)}E_{2}(s)-\alpha_{2}e^{(\widetilde{\gamma}_{2}+\widetilde{\mu}_{2})(s-t)}\widetilde{E}_{2}(s)\right)\,ds$
$\displaystyle\displaystyle=:$ $\displaystyle\displaystyle J_{1}(t)+J_{2}(t),$

where $\displaystyle J_{1}(t)$ and $\displaystyle J_{2}(t)$ are defined by the
first and second term on the right-hand side of the first equality. Firstly we
can use the Mean Value Theorem to obtain

$\displaystyle\displaystyle|J_{1}(t)|=|I_{2}(0)|\left|e^{-(\gamma_{2}+\mu_{2})t}-e^{-(\widetilde{\gamma}_{2}+\widetilde{\mu}_{2})t}\right|\leq|I_{2}(0)||\gamma_{2}+\mu_{2}-\widetilde{\gamma}_{2}-\widetilde{\mu}_{2}|te^{-\min\\{\gamma_{2}+\mu_{2},\widetilde{\gamma}_{2}+\widetilde{\mu}_{2}\\}t}.$

In order to prove that the first term on the right-hand side of the inequality
is bounded by $\displaystyle\varepsilon$, note that the equations for
$\displaystyle R_{1}(t)$ and $\displaystyle D_{1}(t)$ imply

$\displaystyle
R^{\prime}_{1}(t)+D^{\prime}_{1}(t)=(\gamma_{1}+\mu_{1})I_{1}(t)$

and therefore, supposing w.l.g.
$\displaystyle\gamma_{1}+\mu_{1}<\widetilde{\gamma}_{1}+\widetilde{\mu}_{1}$,

$\displaystyle\displaystyle
1<\frac{\widetilde{\gamma}_{1}+\widetilde{\mu}_{1}}{\gamma_{1}+\mu_{1}}$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\frac{\widetilde{R}_{1}(t)+\widetilde{D}_{1}(t)-R_{1}(0)-D_{1}(0)}{R_{1}(t)+D_{1}(t)-R_{1}(0)-D_{1}(0)}\;\frac{\int_{0}^{t}I_{1}(s)\,ds}{\int_{0}^{t}\widetilde{I}_{1}(s)\,ds}$
$\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle\left(1+\frac{2\varepsilon}{R_{1}(t)+D_{1}(t)-R_{1}(0)-D_{1}(0)}\right)\left(1+\frac{\varepsilon
t}{\int_{0}^{t}\widetilde{I}_{1}(s)\,ds}\right)$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\left(1+\frac{2\varepsilon}{R_{1}(t)+D_{1}(t)-R_{1}(0)-D_{1}(0)}\right)\left(1+\frac{\varepsilon\widetilde{\mu}_{1}t}{\widetilde{D}_{1}(t)-\widetilde{D}_{1}(0)}\right),$

where the last equality comes from the equation for
$\displaystyle\widetilde{D}_{1}(t)$.

Since the inequality above holds for every $\displaystyle t>0$, there exist
constants $\displaystyle C_{3}$ and $\displaystyle C_{4}$ such that

$\displaystyle\displaystyle\frac{\widetilde{\gamma}_{1}+\widetilde{\mu}_{1}}{\gamma_{1}+\mu_{1}}$
$\displaystyle\displaystyle\leq$ $\displaystyle\displaystyle
1+C_{3}\varepsilon+C_{4}\varepsilon^{2}$

and then

$\displaystyle\displaystyle
0\leq\widetilde{\gamma}_{1}+\widetilde{\mu}_{1}-\gamma_{1}-\mu_{1}$
$\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle(\gamma_{1}+\mu_{1})\left(C_{3}\varepsilon+C_{4}\varepsilon^{2}\right).$

Since by hypotheses (iii) and (iv) we have
$\displaystyle\widetilde{\gamma}_{2}+\widetilde{\mu}_{2}-\gamma_{2}-\mu_{2}=\widetilde{\gamma}_{1}+\widetilde{\mu}_{1}-\gamma_{1}-\mu_{1}$,
we conclude that

(25)
$\displaystyle\displaystyle|J_{1}(t)|\leq|I_{2}(0)|te^{-\min\\{\gamma_{2}+\mu_{2},\widetilde{\gamma}_{2}+\widetilde{\mu}_{2}\\}t}(\gamma_{1}+\mu_{1})\left(C_{3}\varepsilon+C_{4}\varepsilon^{2}\right).$

Now for the term $\displaystyle J_{2}(t)$, we can decompose it as

(26) $\displaystyle\displaystyle J_{2}(t)=H_{1}(t)+H_{2}(t)+H_{3}(t)$

where

$\displaystyle\displaystyle H_{1}(t)$ $\displaystyle\displaystyle:=$
$\displaystyle\displaystyle(\alpha_{2}-\widetilde{\alpha}_{2})\int_{0}^{t}e^{(\gamma_{2}+\mu_{2})(s-t)}E_{2}(s)\,ds,$
$\displaystyle\displaystyle H_{2}(t)$ $\displaystyle\displaystyle:=$
$\displaystyle\displaystyle\widetilde{\alpha}_{2}\int_{0}^{t}\left(e^{(\gamma_{2}+\mu_{2})(s-t)}-e^{(\widetilde{\gamma}_{2}+\widetilde{\mu}_{2})(s-t)}\right)E_{2}(s)\,ds,$
$\displaystyle\displaystyle H_{3}(t)$ $\displaystyle\displaystyle:=$
$\displaystyle\displaystyle\widetilde{\alpha}_{2}\int_{0}^{t}e^{(\widetilde{\gamma}_{2}+\widetilde{\mu}_{2})(s-t)}\left(E_{2}(s)-\widetilde{E}_{2}(s)\right)\,ds.$

For the $\displaystyle H_{1}(t)$, we can use (22) and obtain

(27) $\displaystyle\displaystyle\left|H_{1}(t)\right|$
$\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle|\alpha_{2}-\widetilde{\alpha}_{2}|\frac{S_{2}(0)-S_{2}(t)+E_{2}(0)-E_{2}(t)}{\alpha_{2}}$
$\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle\left(S_{2}(0)-S_{2}(t)+E_{2}(0)-E_{2}(t)\right)(C_{1}+C_{2}\varepsilon)\varepsilon.$

For the $\displaystyle H_{2}(t)$, we can use a similar argument as the one
used for $\displaystyle J_{1}(t)$ and $\displaystyle H_{1}(t)$:

(28) $\displaystyle\displaystyle\left|H_{2}(t)\right|$
$\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle\widetilde{\alpha}_{2}|\gamma_{2}+\mu_{2}-\widetilde{\gamma}_{2}-\widetilde{\mu}_{2}|\int_{0}^{t}(t-s)e^{\min\\{\gamma_{2}+\mu_{2},\widetilde{\gamma}_{2}+\widetilde{\mu}_{2}\\}(s-t)}E_{2}(s)\,ds$
$\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle\widetilde{\alpha}_{2}\frac{|\gamma_{2}+\mu_{2}-\widetilde{\gamma}_{2}-\widetilde{\mu}_{2}|}{\min\\{\gamma_{2}+\mu_{2},\widetilde{\gamma}_{2}+\widetilde{\mu}_{2}\\}e}\int_{0}^{t}E_{2}(s)\,ds$
$\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle\frac{\widetilde{\alpha}_{2}}{\alpha_{2}}\frac{|\gamma_{2}+\mu_{2}-\widetilde{\gamma}_{2}-\widetilde{\mu}_{2}|}{\min\\{\gamma_{2}+\mu_{2},\widetilde{\gamma}_{2}+\widetilde{\mu}_{2}\\}e}\left(S_{2}(0)-S_{2}(t)+E_{2}(0)-E_{2}(t)\right)$
$\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle\frac{\widetilde{\alpha}_{2}}{\alpha_{2}}\,\frac{S_{2}(0)-S_{2}(t)+E_{2}(0)-E_{2}(t)}{\min\\{\gamma_{2}+\mu_{2},\widetilde{\gamma}_{2}+\widetilde{\mu}_{2}\\}e}(\gamma_{1}+\mu_{1})\left(C_{3}\varepsilon+C_{4}\varepsilon^{2}\right).$

Finally, we use (23) to estimate $\displaystyle H_{3}(t)$ as

(29)
$\displaystyle\displaystyle\left|H_{3}(t)\right|\leq\frac{\widetilde{\alpha}_{2}}{\widetilde{\gamma}_{2}-\widetilde{\mu}_{2}}(1-e^{-(\widetilde{\gamma}_{2}+\widetilde{\mu}_{2})t})K_{E}\varepsilon$

and, from (24), (25), (26), (27), (28) and (29), we obtain

$\displaystyle\displaystyle\left|I_{2}(t)-\widetilde{I}_{2}(t)\right|$
$\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle|I_{2}(0)|te^{-\min\\{\gamma_{2}+\mu_{2},\widetilde{\gamma}_{2}+\widetilde{\mu}_{2}\\}t}(\gamma_{1}+\mu_{1})\left(C_{3}\varepsilon+C_{4}\varepsilon^{2}\right)$
$\displaystyle\displaystyle+\left(S_{2}(0)-S_{2}(t)+E_{2}(0)-E_{2}(t)\right)(C_{1}+C_{2}\varepsilon)\varepsilon$
$\displaystyle\displaystyle+\frac{\widetilde{\alpha}_{2}}{\alpha_{2}}\,\frac{S_{2}(0)-S_{2}(t)+E_{2}(0)-E_{2}(t)}{\min\\{\gamma_{2}+\mu_{2},\widetilde{\gamma}_{2}+\widetilde{\mu}_{2}\\}e}(\gamma_{1}+\mu_{1})\left(C_{3}\varepsilon+C_{4}\varepsilon^{2}\right)$
$\displaystyle\displaystyle+\frac{\widetilde{\alpha}_{2}}{\widetilde{\gamma}_{2}-\widetilde{\mu}_{2}}(1-e^{-(\widetilde{\gamma}_{2}+\widetilde{\mu}_{2})t})K_{E}\varepsilon.$

Therefore, we have $\displaystyle\|I_{2}-\widetilde{I}\|_{[0,T]}\leq
K_{I}(\varepsilon)\varepsilon$ where

$\displaystyle\displaystyle K_{I}(\varepsilon)$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\left(|I_{2}(0)|+\frac{\widetilde{\alpha}_{2}}{\alpha_{2}}\left((S_{2}(0)-S_{2}(T)+E_{2}(0)\right)\right)\frac{(\gamma_{1}+\mu_{1})\left(C_{3}+C_{4}\varepsilon\right)}{\min\\{\gamma_{2}+\mu_{2},\widetilde{\gamma}_{2}+\widetilde{\mu}_{2}\\}e}$
$\displaystyle\displaystyle+\left(S_{2}(0)-S_{2}(T)+E_{2}(0)\right)(C_{1}+C_{2}\varepsilon)+\frac{\widetilde{\alpha}_{2}}{\widetilde{\gamma}_{2}-\widetilde{\mu}_{2}}K_{E}.$

Step 4: For the solutions $\displaystyle R_{2}$ and
$\displaystyle\widetilde{R}_{2}$, we use (2)

$\displaystyle
S_{2}+E_{2}+I_{2}+R_{2}+D_{2}=N_{2}=\widetilde{S}_{2}+\widetilde{E}_{2}+\widetilde{I}_{2}+\widetilde{R}_{2}+\widetilde{D}_{2}$

and thus

$\displaystyle\|R_{2}-\widetilde{R}_{2}\|_{[0,T]}\leq\left(K_{S}+K_{E}+K_{I}+1\right)\varepsilon.$

∎

######  Remark 3.5.

In many cases the hypothesis ii-iv can be retrieved from clinical data under
the assumption that some of the parameters are equal for both populations.
Hence, in general, condition i seems to be the most difficult to verify for
practical applications.

######  Remark 3.6.

Note that, given a fixed model with parameters
$\displaystyle\widetilde{\mathcal{A}}$,
$\displaystyle\widetilde{\mathcal{B}}$, $\displaystyle\widetilde{\Gamma}$,
$\displaystyle\widetilde{\mathcal{M}}$, the conclusion of the theorem holds
only for all the solutions $\displaystyle S_{2},E_{2},I_{2},R_{2},D_{2}$ to
models whose parameters $\displaystyle\mathcal{A},\mathcal{B}$,
$\displaystyle\Gamma$, $\displaystyle\mathcal{M}$ satisfy (i)-(iv) with
respect to the fixed ones.

######  Remark 3.7.

Instead of using the equation for $\displaystyle D_{2}$ in Step 1, we could
use the equation for $\displaystyle R_{2}$. Then, the estimates would depend
on $\displaystyle\frac{\beta_{2i}-\beta_{1i}}{\gamma_{i}}$ instead of
$\displaystyle\frac{\beta_{2i}-\beta_{1i}}{\mu_{i}}$. In this case, hypotheses
(i) and (v) would change. We choose to base our estimates on death rates
because data on number of deaths are more reliable for applications.

######  Remark 3.8.

Alternatively, if
$\displaystyle\mathcal{B}(\Gamma+\mathcal{M})^{-1}=\widetilde{\mathcal{B}}(\widetilde{\Gamma}+\widetilde{\mathcal{M}})^{-1}$,
i.e., if

$\displaystyle\frac{\beta_{ij}}{\gamma_{j}+\mu_{j}}=\frac{\widetilde{\beta}_{ij}}{\widetilde{\gamma}_{j}+\widetilde{\mu}_{j}}\mbox{
for }i,j=1,2\;,$

we can prove the estimates in Step 3 by using (7) and then

$\displaystyle\displaystyle\frac{\beta_{11}}{\gamma_{1}+\mu_{1}}(S_{1}(t)+E_{1}(t)+I_{1}(t))+\frac{\beta_{12}}{\gamma_{2}+\mu_{2}}(S_{2}(t)+E_{2}(t)+I_{2}(t))$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle P_{1},$
$\displaystyle\displaystyle\frac{\widetilde{\beta}_{11}}{\widetilde{\gamma}_{1}+\widetilde{\mu}_{1}}(\widetilde{S}_{1}(t)+\widetilde{E}_{1}(t)+\widetilde{I}_{1}(t))+\frac{\widetilde{\beta}_{12}}{\widetilde{\gamma}_{2}+\widetilde{\mu}_{2}}(\widetilde{S}_{2}(t)+\widetilde{E}_{2}(t)+\widetilde{I}_{2}(t))$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle\widetilde{P}_{1},$

where $\displaystyle P_{1}$ and $\displaystyle\widetilde{P}_{1}$ are constants
depending on the matrices $\displaystyle\mathcal{B}(\Gamma+\mathcal{M})^{-1}$
and
$\displaystyle\widetilde{\mathcal{B}}(\widetilde{\Gamma}+\widetilde{\mathcal{M}})^{-1}$,
and on the initial conditions. By hypothesis, $\displaystyle
P_{1}=\widetilde{P}_{1}$. Then, from these two identities, we obtain

$\displaystyle\displaystyle I_{2}(t)-\widetilde{I}_{2}(t)$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\frac{\gamma_{2}+\mu_{2}}{\beta_{12}}\ln\left(\frac{S_{1}(t)}{\widetilde{S}_{1}(t)}\right)-S_{2}(t)+\widetilde{S}_{2}(t)-E_{2}(t)+\widetilde{E}_{2}(t)$
$\displaystyle\displaystyle-\frac{\gamma_{2}+\mu_{2}}{\gamma_{1}+\mu_{1}}\frac{\beta_{11}}{\beta_{12}}\left(S_{1}(t)-\widetilde{S}_{1}(t)+E_{1}(t)-\widetilde{E}_{1}(t)+I_{1}(t)-\widetilde{I}_{1}(t)\right),$

and we can estimate

$\displaystyle\displaystyle\left|I_{2}(t)-\widetilde{I}_{2}(t)\right|$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\left(\frac{\gamma_{2}+\mu_{2}}{\beta_{12}}\frac{1}{\widetilde{S}_{1}(t)}+3\frac{\gamma_{2}+\mu_{2}}{\gamma_{1}+\mu_{1}}\frac{\beta_{11}}{\beta_{12}}+K_{S}+K_{E}\right)\varepsilon.$

######  Remark 3.9.

Although the argument in Step 4 is enough for the proof of the theorem, the
estimate is not very good since it considers each of the previous classes
separately. One can improve the the estimate by noticing from the
$\displaystyle R$ and $\displaystyle D$ equations in (3) that

$\displaystyle
R_{2}^{\prime}(t)=\gamma_{2}I_{2}(t)=\frac{\gamma_{2}}{\mu_{2}}D_{2}^{\prime}(t)$

and therefore

$\displaystyle
R_{2}(t)-R_{2}(0)=\frac{\gamma_{2}}{\mu_{2}}(D_{2}(t)-D_{2}(0))\;,\;\;\forall
t\geq 0.$

Thus, using the same argument for $\displaystyle\widetilde{R}_{2}(t)$, we
obtain

(30)
$\left|R_{2}(t)-\widetilde{R}_{2}(t)\right|\leq\left|\frac{\gamma_{2}}{\mu_{2}}-\frac{\widetilde{\gamma}_{2}}{\widetilde{\mu}_{2}}\right|\left|D_{2}(t)-D_{2}(0)\right|+\frac{\widetilde{\gamma}_{2}}{\widetilde{\mu}_{2}}\left|D_{2}(t)-\widetilde{D}_{2}(t)\right|,$

and thus using the hypotheses (iii), (iv) and (v), we obtain the estimate.

######  Remark 3.10.

We could obtain better estimate for the distances between solutions for the
second group by considering separate distances on the first group. For
example, defining

(31)
$\displaystyle\varepsilon_{S}:=\|S_{1}-\widetilde{S}_{1}\|_{[0,T]}\,,\;\varepsilon_{E}:=\|E_{1}-\widetilde{E}_{1}\|_{[0,T]}\,,\;\varepsilon_{I}:=\|I_{1}-\widetilde{I}_{1}\|_{[0,T]},$
$\displaystyle\varepsilon_{R}:=\|R_{1}-\widetilde{R}_{1}\|_{[0,T]}\,,\;\varepsilon_{D_{1}}:=\|D_{1}-\widetilde{D}_{1}\|_{[0,T]},\;\varepsilon_{D_{2}}:=\|D_{2}-\widetilde{D}_{2}\|_{[0,T]},$

we can rewrite the estimate for the susceptible class in the second group as
$\displaystyle\|S_{2}-\widetilde{S}_{2}\|_{[0,T]}\leq\mathcal{K}_{S}(T)$ where

(32)
$\mathcal{K}_{S}(T):={S}_{2}(T)\left(1-\left(1-\frac{\varepsilon_{S}}{S_{1}(T)}\right)e^{-\left(|V_{1}|\varepsilon_{D_{1}}+|V_{2}|\varepsilon_{D_{2}}\right)}\right).$

Analogously, we have
$\displaystyle\|E_{2}-\widetilde{E}_{2}\|_{[0,T]}\leq\mathcal{K}_{E}(T)$ and
$\displaystyle\|I_{2}-\widetilde{I}_{2}\|_{[0,T]}\leq\mathcal{K}_{I}(T)$ where

$\displaystyle\displaystyle\mathcal{K}_{E}(T)$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle
E_{2}(0)|\widetilde{\alpha}_{2}-\alpha_{2}|te^{-\alpha_{2}t}$
$\displaystyle\displaystyle+(S_{2}(0)-S_{2}(T))\max\left\\{0,\ln\left(\frac{\widetilde{\alpha}_{2}}{\alpha_{2}}\right)\right\\}+\mathcal{K}_{S}(T)$

and

(34) $\displaystyle\mathcal{K}_{I}(T)$
$\displaystyle=\frac{|\gamma_{2}+\mu_{2}-\widetilde{\gamma}_{2}-\widetilde{\mu}_{2}|}{\min\\{\gamma_{2}+\mu_{2},\widetilde{\gamma}_{2}+\widetilde{\mu}_{2}\\}e}\left(I_{2}(0)+\frac{\widetilde{\alpha}_{2}}{\alpha_{2}}\left(S_{2}(0)-S_{2}(T)+E_{2}(0)\right)\right)$
$\displaystyle+\left|1-\frac{\widetilde{\alpha}_{2}}{\alpha_{2}}\right|\left(S_{2}(0)-S_{2}(T)+E_{2}(0)\right)\frac{\widetilde{\alpha}_{2}}{\widetilde{\gamma}_{2}-\widetilde{\mu}_{2}}\|E_{2}-\widetilde{E}_{2}\|_{[0,T]}.$

Note that, by the proof of the Theorem 3.4, the difference
$\displaystyle|\widetilde{\alpha_{2}}-\alpha_{2}|$ is controlled by
$\displaystyle\varepsilon_{S}$ and $\displaystyle\varepsilon_{E}$. And the
difference
$\displaystyle|\gamma_{2}+\mu_{2}-\widetilde{\gamma}_{2}-\widetilde{\mu}_{2}|$
is controlled by $\displaystyle\varepsilon_{R}$ and
$\displaystyle\varepsilon_{D}$.

Finally, for the class of recovered individuals,
$\displaystyle\|R_{2}-\widetilde{R}_{2}\|_{[0,T]}$, we can use Remark 3.9 to
obtain

(35)
$\|R_{2}-\widetilde{R}_{2}\|_{[0,T]}\leq\mathcal{K}_{R}(T):=\left|\frac{\gamma_{2}}{\mu_{2}}D_{2}(T)-\frac{\widetilde{\gamma}_{2}}{\widetilde{\mu}_{2}}\widetilde{D}_{2}(T)+\left(\frac{\widetilde{\gamma}_{2}}{\widetilde{\mu}_{2}}-\frac{\gamma_{2}}{\mu_{2}}\right)D_{2}(0)\right|.$

######  Remark 3.11.

Consider the cumulative number of infectious individuals in both groups for
the first set of parameters

(36) $C_{1}(t):=I_{1}(t)+R_{1}(t)+D_{1}(t)\quad\mbox{ and }\quad
C_{2}(t):=I_{2}(t)+R_{2}(t)+D_{2}(t).$

and the respective functions for the second set of parameters

(37)
$\widetilde{C}_{1}(t)=\widetilde{I}_{1}(t)+\widetilde{R}_{1}(t)+\widetilde{D}_{1}(t)\quad\mbox{
and
}\quad\widetilde{C}_{2}(t):=\widetilde{I}_{2}(t)+\widetilde{R}_{2}(t)+\widetilde{D}_{2}(t).$

Estimating the difference
$\displaystyle(S_{2}+E_{2})(t)-(\widetilde{S}_{2}+\widetilde{E}_{2})(t)$ as we
did in Step 2 of the proof of Theorem 3.4, and using (2), we obtain

(38)
$\|C_{2}-\widetilde{C}_{2}\|_{[0,T]}=\|S_{2}+E_{2}-\widetilde{S}_{2}-\widetilde{E}_{2}\|_{[0,T]}\leq
K_{E}(T).$

Notice that hypothesis (i) is related to social interaction between the two
groups, since it involves $\displaystyle\beta_{ij}$. On the other hand the
hypotheses (ii), (iii) and (iv) are on the parameters
$\displaystyle\alpha_{i}$, $\displaystyle\gamma_{i}$ and
$\displaystyle\mu_{i}$, which depend on biological features of the disease,
and also on how each group deals with the infected individuals, such as access
to the health system and access to medicines.Thus, these conditions may play a
role whether the definition of groups 1 and 2 refer to social or biological
features.

##  4\. Examples with the COVID-19

In this section we propose a method for applying the theoretical results
presented in Section 3. The method is developed in the Subsection 4.1 and
applied to New York State early epidemics of COVID-19 for validation of its
predictions, since in this case data is fully known. In Subsection 4.2 we
apply the same method to the COVID-19 epidemics in the cities of Juazeiro and
Petrolina, which are neighboring cities in the northeast of Brazil that face
different health policies. Our goal with this example is to detect a potential
problem in the number of reported cases, since the city of Juazeiro is under a
different health system and applied proportionally much less tests for
detecting infected individuals compared to the city of Petrolina.

###  4.1. New York County

In this subsection, we show a numerical simulation of the results obtained in
Section 3 based on data for the early stage of the COVID-19 epidemic in the
New York County and neighboring counties. We will illustrate the estimates in
Theorem 3.4 using data for the number of cases and deaths in the neighboring
counties and deaths in the New York County to estimate the curve of infection
in the New York County population and an error range for this estimation. The
choice of these populations was made only by geographic proximity.

We consider as Group 1 the inhabitants of the following counties: Richmond,
Kings, Queens, Bronx, Monmouth, Hudson and Bergen. And as Group 2 the
inhabitants of New York County. Due to the restrictions imposed by the “New
York State on Pause” Executive Order, which affected social distancing and
economic activities from late March 2020 onwards, and thus changed the
interaction rates between groups, we restrict the analysis to the first 30
days of the outbreak of the disease: from March 5 (first day with registered
cases in both groups) till April 3. Therefore, we will use the data presented
in Table 1.

Day (Mar) | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 | 10 | 11 | 12 | 13
---|---|---|---|---|---|---|---|---|---|---|---|---
Cases Group 1 | 0 | 0 | 0 | 2 | 2 | 5 | 7 | 14 | 17 | 17 | 46 | 66
Cases Group 2 | 1 | 1 | 1 | 3 | 4 | 8 | 8 | 11 | 17 | 18 | 39 | 48
Deaths Group 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1
Deaths Group 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
Day (Mar) | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | 23 | 24 | 25
Cases Group 1 | 109 | 151 | 246 | 475 | 1017 | 1833 | 2449 | 3887 | 4877 | 6497 | 7925 | 11364
Cases Group 2 | 71 | 72 | 111 | 277 | 590 | 1038 | 1314 | 1863 | 2072 | 2572 | 2887 | 3616
Deaths Group 1 | 1 | 2 | 4 | 4 | 4 | 6 | 7 | 10 | 46 | 81 | 125 | 173
Deaths Group 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 10 | 19 | 35 | 43
Day (Mar-Apr) | 26 | 27 | 28 | 29 | 30 | 31 | 01 | 02 | 03 | 04 | 05 | 06
Cases Group 1 | 13764 | 15613 | 18986 | 21014 | 24526 | 27382 | 30575 | 34036 | 39407 | 43801 | 47325 | 51323
Cases Group 2 | 4046 | 4478 | 5237 | 5582 | 6060 | 6539 | 7022 | 7398 | 8452 | 9300 | 9898 | 10440
Deaths Group 1 | 229 | 287 | 418 | 491 | 570 | 712 | 935 | 1123 | 1310 | 1573 | 1775 | 2012
Deaths Group 2 | 55 | 65 | 93 | 103 | 119 | 129 | 165 | 178 | 215 | 264 | 390 | 436

Table 1. Cumulative daily reported cases and deaths confirmed by testing from
March 3, 2020 to April 5, 2020, for New York County (Group 2), and Richmond,
Kings, Queens, Bronx, Monmouth, Hudson and Bergen counties (Group 1) according
to [5].

In order to illustrate Theorem 3.4, we are assuming that data in Table 1 comes
from a model (1) (possibly with noises) for some unknown set of parameters.
This set of unknown parameters will be used as the first set mentioned in
Theorem 3.4 and will be represented as $\displaystyle\widetilde{\mathcal{A}}$,
$\displaystyle\widetilde{\mathcal{B}}$, $\displaystyle\widetilde{\Gamma}$,
$\displaystyle\widetilde{\mathcal{M}}$. Also, the respective solutions to the
model will be indicated by tilde. Therefore, from the data we will obtain the
(assumed) model solutions $\displaystyle\widetilde{D}_{1}$ and
$\displaystyle\widetilde{D}_{2}$, and the cumulative function
$\displaystyle\widetilde{C}_{1}:=\widetilde{S}_{1}+\widetilde{E}_{1}+\widetilde{I}_{1}$
(we give more details later). The compartments
$\displaystyle\widetilde{S}_{1},\widetilde{E}_{1}$ and
$\displaystyle\widetilde{R}_{1}$ will be extrapolated from these two. The
solutions $\displaystyle\widetilde{S}_{2},\widetilde{E}_{2},\widetilde{I}_{2}$
and $\displaystyle\widetilde{R}_{2}$ are assumed to be unknown. We can think
of it as if the data for group 2 is full of measurement errors and thus
unreliable for the number of cases, but reliable for number of deaths. Thus,
we will use the data on the evolution of cases of infection and deaths in
Group 1 and also the evolution of deaths in Group 2 to estimate the parameters
for the model. These estimated parameters will be used as the second set of
parameters mentioned in Theorem 3.4 and will be represented as
$\displaystyle\mathcal{A},\mathcal{B}$, $\displaystyle\Gamma$,
$\displaystyle\mathcal{M}$, without tilde, with the respective solutions
represented also without the tilde. For the group 1, we will have the
solutions obtained from the data (with the tilde), since the parameters are
unknown, and the solutions obtained from the model with estimated parameters
(without the tilde), and thus we can measure the distances between them, as
presented in (17). For the group 2, the same can be done, only for the
compartment $\displaystyle D$, since the only reliable data for group 2 are
the number of deaths, and thus will be measured as in (v) of Theorem 3.4.
Therefore, we will use Theorem 3.4 to estimate the distance between the
solutions for the other compartments of group 2. Thus, we will be able to
measure how far the solutions for the estimated model (without the tilde) are
from the unknown solutions (with tilde). The theorem gives a measure on how
the curves obtained by plotting the solutions of the estimated model
accurately depict the evolution of the cases in Group 2, since it gives us an
error estimate where the actual curve must be.

The hypotheses of Theorem 3.4 are imposed in the fitting process as follows.
Firstly, for the hypothesis (ii) we will suppose that the latency rate is the
same for both groups, i.e. $\displaystyle\alpha_{1}=\alpha_{2}=\alpha$, since
this is a parameter related exclusively to the disease itself. Furthermore, by
the same reason, after fitting this parameter to the dataset, we will consider
only solutions with this same latency. This means that we will consider also
$\displaystyle\widetilde{\alpha}_{1}=\widetilde{\alpha}_{2}=\alpha$. We will
make the same assumption with respect to the hypothesis iii. Therefore we will
consider that
$\displaystyle\widetilde{\gamma}_{1}=\widetilde{\gamma}_{2}=\gamma=\gamma_{1}=\gamma_{2}$.
We will not impose $\displaystyle\mu_{1}=\mu_{2}$ because the mortality rates
may not depend only on the disease itself, but be influenced by the quality of
medical treatment and access to it in each group.

Let us clarify how the data from Table 1 will be used. Let
$\displaystyle\widetilde{C}_{1}$ and $\displaystyle\widetilde{D}_{1}$ be the 7
days centered moving average of the cumulative number of cases and cumulative
number of deaths in Group 1, respectively, as given in Table 1. And
$\displaystyle\widetilde{C}_{2}$ and $\displaystyle\widetilde{D}_{2}$ be the 7
days centered moving average of the cumulative number of cases and cumulative
number of deaths in Group 2, respectively as given in Table 1. Thus, in order
to fit the parameters for the SEIRD model, we will use the values of
$\displaystyle\widetilde{C}_{1}$, $\displaystyle\widetilde{D}_{1}$ and
$\displaystyle\widetilde{D}_{2}$ only. After fitting them, we will use the
solution of the fitted SEIRD model for the number of cases in Group 2 and
Theorem 3.4 to estimate $\displaystyle\widetilde{C}_{2}$.

From the functions defined above, we can define the approximated number of
susceptible individuals in each group by

(39)
$\widetilde{S}_{i}(t)=N_{i}-\frac{\widetilde{C}_{i}(t+1)-\widetilde{C}_{i}(t)}{\alpha}-\widetilde{C}_{i}(t),\qquad
i=1,2$

since the second term on the right-hand side approximates the number of
exposed individuals.

For the hypothesis (i), we will estimate the constants $\displaystyle
V_{1}=(\widetilde{\beta}_{21}-\widetilde{\beta}_{11})/\widetilde{\mu}_{1}$ and
$\displaystyle
V_{2}=(\widetilde{\beta}_{22}-\widetilde{\beta}_{12})/\widetilde{\mu}_{2}$ as
follows: from Lemma 3.1, we have

$\displaystyle\displaystyle\frac{\widetilde{S}_{1}(t)}{\widetilde{S}_{1}(0)}=\frac{\widetilde{S}_{2}(t)}{\widetilde{S}_{2}(0)}\exp\left\\{\int_{0}^{t}\left(\left(\widetilde{\beta}_{21}-\widetilde{\beta}_{11}\right)\widetilde{I}_{1}(s)+\left(\widetilde{\beta}_{22}-\widetilde{\beta}_{12}\right)\widetilde{I}_{2}(s)\right)\,ds\right\\}.$

Since we are not using the data on the cases of infection in Group 2, we need
to estimate $\displaystyle\widetilde{S}_{2}(t)/\widetilde{S}_{2}(0)$. Note
that for small $\displaystyle t$ we have
$\displaystyle\widetilde{S}_{2}(t)/\widetilde{S}_{2}(0)\approx 1$, by
continuity. Therefore, using the equations for
$\displaystyle\widetilde{S}_{1}$ and $\displaystyle\widetilde{S}_{2}$, for
small $\displaystyle t$ we have

(40)
$\frac{\widetilde{S}_{1}(t)}{\widetilde{S}_{1}(0)}\approx\exp\left\\{V_{1}(\widetilde{D}_{1}(t)-\widetilde{D}_{1}(0))+V_{2}(\widetilde{D}_{2}(t)-\widetilde{D}_{2}(0))\,ds\right\\}.$

For a pair of time values $\displaystyle t=t_{1}$ and $\displaystyle t=t_{2}$,
the equation (40) generates the linear system, that has solutions:

(41) $\displaystyle\displaystyle
V_{1}\approx\displaystyle\frac{\widetilde{D}_{2}(t_{2})\ln\left(\frac{\widetilde{S}_{1}(t_{1})}{\widetilde{S}_{1}(0)}\right)-\widetilde{D}_{2}(t_{1})\ln\left(\frac{\widetilde{S}_{1}(t_{2})}{\widetilde{S}_{1}(0)}\right)}{\widetilde{D}_{1}(t_{1})\widetilde{D}_{2}(t_{2})-\widetilde{D}_{1}(t_{2})\widetilde{D}_{2}(t_{1})};$
(42) $\displaystyle\displaystyle
V_{2}\approx\displaystyle\frac{\widetilde{D}_{1}(t_{1})\ln\left(\frac{\widetilde{S}_{1}(t_{2})}{\widetilde{S}_{1}(0)}\right)-\widetilde{D}_{1}(t_{2})\ln\left(\frac{\widetilde{S}_{1}(t_{1})}{\widetilde{S}_{1}(0)}\right)}{\widetilde{D}_{1}(t_{1})\widetilde{D}_{2}(t_{2})-\widetilde{D}_{1}(t_{2})\widetilde{D}_{2}(t_{1})},$

as long as
$\displaystyle\widetilde{D}_{1}(t_{1})\widetilde{D}_{2}(t_{2})-\widetilde{D}_{1}(t_{2})\widetilde{D}_{2}(t_{1})\not=0$.
The quotients $\displaystyle\widetilde{S}_{1}(t_{j})/\widetilde{S}_{1}(0)$,
$\displaystyle j=1,2$, can be approximated from the data using (39).

Therefore, one can calculate $\displaystyle V_{1}$ and $\displaystyle V_{2}$
for any pair of time values $\displaystyle t=t_{1}$ and $\displaystyle
t=t_{2}$. We took the average of the respective absolute values of
$\displaystyle V_{1}$ and $\displaystyle V_{2}$ for all pairs of integer times
$\displaystyle t$ between $\displaystyle t=1$ and $\displaystyle t=30$ days
and then set the signs such that $\displaystyle V_{1}<0$ and $\displaystyle
V_{2}>0$. Thus we obtain $\displaystyle V_{1}=-3.023433\times 10^{-7}$ and
$\displaystyle V_{2}=1.991887\times 10^{-5}$. These values are used to fix the
relation between the transmission rates $\displaystyle\beta_{ij}$ and the
mortality rates, due to hypothesis (i). Hence, given $\displaystyle\beta_{i1}$
and $\displaystyle\mu_{i}$, the corresponding value for
$\displaystyle\beta_{i2}$ is given by
$\displaystyle\beta_{i2}=\beta_{i1}-\mu_{i}V_{i}$. We observed that taking the
average of (41)-(42) directly led to unrealistic values of $\displaystyle
V_{i}$. Furthermore, using the absolute values to estimate $\displaystyle
V_{1}$ and $\displaystyle V_{2}$ and setting them to have opposite signs gave
best adjustments to the data. That is compatible with the fact that the local
transmissions are stronger than the infection between individuals of different
groups.

We consider the following set of parameters: $\displaystyle\beta_{11},$
$\displaystyle\beta_{12}$, $\displaystyle\beta_{21},$
$\displaystyle\beta_{22},$ $\displaystyle\alpha$, $\displaystyle\gamma_{1}$,
$\displaystyle\gamma_{2}$, $\displaystyle\mu_{1}$, $\displaystyle\mu_{2}$,
$\displaystyle E_{1}(0)$, $\displaystyle E_{2}(0)$. We are considering the
initial conditions as a parameter too because it is not clear from the dataset
how to obtain a good approximation for the evolution of the exposed class.
Thus, we will use the initial conditions for the number of susceptible
individuals in each group as $\displaystyle S_{1}(0)=N_{1}-I_{1}(0)-E_{1}(0)$
and $\displaystyle S_{2}(0)=N_{2}-I_{2}(0)-E_{2}(0)$, where $\displaystyle
N_{1}=0.84N$ and $\displaystyle N_{2}=0.16N$, with $\displaystyle N=10.5\times
10^{6}$ being the total population of Group 1 and 2 together (see [1]).

Under the above conditions, we fitted all the parameters mentioned above
together using the Weighted Root Mean Square Error (WRMSE). The weights were
used to compensate for the fact that the first group is much bigger.
Therefore, we used weight 1 for the cases in Group 1, and weights 40 and 240
for deaths in Group 1 and Group 2, respectively. The numeric calculations were
made by the software R-CRAN, using the packages deSolve ([34]) and DEoptim
([28]) for solving the ODEs and optimizing all the parameters together,
respectively. In order to reduce overfitting, we have already set
$\displaystyle\alpha_{1}=\alpha_{2}=\alpha$,
$\displaystyle\gamma_{1}=\gamma_{2}$, and fixed the relation between
$\displaystyle\beta_{1i},\beta_{2i}$ and $\displaystyle\mu_{i}$ through
$\displaystyle V_{i}$, $\displaystyle i=1,2$, but we also impose
$\displaystyle\alpha=0.3$ (value chosen from the range observed in [4]) and a
restriction to the optimizing process by considering only parameters
satisfying $\displaystyle\beta_{i,j}\in(0,10/N)$, $\displaystyle i,j=1,2$, and
$\displaystyle\gamma\in(1/14,1/4)$. These bounds are based on the observed
values in [4, 21, 22, 38]. The Table 2 presents the fitted parameters for the
SEIRD system (3).

Parameters
---
Symbol | Value | Symbol | Value | Symbol | Value
$\displaystyle\beta_{11}$ | $\displaystyle 5.817\times 10^{-8}$ | $\displaystyle\mu_{1}$ | 0.010349 | $\displaystyle E_{1}(0)$ | 94.64
$\displaystyle\beta_{12}$ | $\displaystyle 5.609\times 10^{-8}$ | $\displaystyle\mu_{2}$ | 0.006640 | $\displaystyle E_{2}(0)$ | 81.69
$\displaystyle\beta_{21}$ | $\displaystyle 5.504\times 10^{-8}$ | $\displaystyle\gamma$ | 0.192012 |  |
$\displaystyle\beta_{22}$ | $\displaystyle 1.883\times 10^{-7}$ | $\displaystyle\alpha$ | 0.3 |  |
Table 2. Fitted parameters for the system (3) for the data of reported cases
and deaths in Group 1 and reported deaths in Group 2 as presented in Table 1.
The measurement units of the parameters are: $\displaystyle\beta_{ij}$ is in
$\displaystyle(\mbox{individuals}\cdot\mbox{days})^{-1}$,
$\displaystyle\mu_{i},\gamma$ and $\displaystyle\alpha$ are in
$\displaystyle(\mbox{days})^{-1}$, and $\displaystyle E_{i}(0)$ are in
individuals, for all $\displaystyle i,j\in\\{1,2\\}$.

Let $\displaystyle S_{1},S_{2},E_{1},E_{2},I_{1},I_{2},R_{1},R_{2},D_{1}$ and
$\displaystyle D_{2}$ the solutions to (1) with the parameters showed in Table
2. Then the function of cumulative number of infections individuals for both
groups is

(43) $C_{1}(t):=I_{1}(t)+R_{1}(t)+D_{1}(t)\quad\mbox{ and }\quad
C_{2}(t):=I_{2}(t)+R_{2}(t)+D_{2}(t).$

The Figure 2 shows the graphs of $\displaystyle C_{1}(t)$, $\displaystyle
C_{2}(t)$, $\displaystyle D_{1}(t)$ and $\displaystyle D_{2}(t)$ with the
respective cumulative number of reported cases and deaths obtained from the
dataset.

(a) Infectious cases in Group 1

(b) Deaths in Group 1

(c) Infectious cases in Group 2

(d) Deaths in Group 2

Figure 2. Fit of the SEIRD system (3) to the reported data. In (2(a)) we plot
the cumulative number of reported infectious cases (black dots) and the
function $\displaystyle C_{1}(t)$ (blue solid line) for the Group 1
(neighboring counties). In (2(b)), the reported deaths (black dots) and the
function $\displaystyle D_{1}(t)$ (red solid line) for the Group 1. In (2(c))
we plot the cumulative number of reported infectious cases (black dots) and
the function $\displaystyle C_{2}(t)$ (blue solid line) for the Group 2 (New
York County). In figure (2(d)), the reported deaths (black dots) and the
function $\displaystyle D_{2}(t)$ (red solid line) for the Group 2.

In order to obtain the estimate in Theorem 3.4, let us define for every
$\displaystyle t=1,...,30$ the distances

$\displaystyle\epsilon_{D_{1}}(t)=\max_{s=1,\ldots,t}|D_{1}(s)-\widetilde{D}_{1}(s)|\;,\;\;\epsilon_{D_{2}}(t)=\max_{s=1,\ldots,t}|D_{2}(s)-\widetilde{D}_{2}(s)|$

and

$\displaystyle\epsilon_{S_{1}}(t)=\max_{s=1,\ldots,t}|S_{1}(s)-\widetilde{S}_{1}(s)|\;,$

where $\displaystyle\widetilde{D}_{1},\widetilde{D}_{2}$ and
$\displaystyle\widetilde{S}_{1}$ are the functions defined at the beginning of
this example.

Therefore, from Remarks 3.10 and 3.11, we obtain that
$\displaystyle\widetilde{C}_{2}$ must satisfy

$\displaystyle\max_{s=1,\ldots,t}|C_{2}(s)-\widetilde{C}_{2}(s)|\leq K_{C}(t)$

where

$\displaystyle
K_{C}(t):=\max_{s=1,\ldots,t}\left\\{{S}_{2}(s)\left(1-\left(1-\frac{\varepsilon_{S_{1}}(s)}{S_{1}(s)}\right)e^{-\left(|V_{1}|\varepsilon_{D_{1}}(s)+|V_{2}|\varepsilon_{D_{2}}(s)\right)}\right)\right\\}.$

The Figure 3 presents the functions $\displaystyle C_{1}(t)$ and
$\displaystyle C_{2}(t)$, the reported cases functions
$\displaystyle\widetilde{C}_{1}(t)$ and $\displaystyle\widetilde{C}_{2}(t)$
and the range of estimated distance $\displaystyle K_{C}(t)$.

For the last day of the period, $\displaystyle t=30$, the values of the
estimated distances for the cumulative function of infectious individuals is
$\displaystyle K_{C}(30)=1713.09$, which represents $\displaystyle 20.78\%$ of
$\displaystyle C_{2}(30)$. In comparison, the fitted curve $\displaystyle
C_{1}(t)$ has an accumulated error from the reported data by 2890.93 cases,
which represents $\displaystyle 9.96\%$ of $\displaystyle C_{1}(30)$.

Figure 3. In this figure we plot the functions cumulative number of infectious
cases $\displaystyle C_{1}(t)$ (light blue solid line) and $\displaystyle
C_{2}(t)$ (dark blue solid line), the reported number of infectious cases for
Group 1 (circle markers) and for Group 2 (square markers), and the estimated
distance range function $\displaystyle K_{C}(t)$ (blue area).

It is important to note that although the black square markers (Group 2 data
on number of infectious cases) in Figure 3 were not used in the fitting
process for the curves, they lie within the estimated interval predicted by
the Theorem 3.4. Indeed, since the reported number of infectious cases in
Group 2 fall into the estimated distance from $\displaystyle C_{2}(t)$, this
suggests that the evolution of cases follows a two group SEIRD model and the
values in (i)-(iv) were well approximated.

###  4.2. Petrolina and Juazeiro Cities

Petrolina and Juazeiro are two neighboring cities in the northeast of Brazil,
with the first one located in the state of Pernambuco, and the latter located
in the state of Bahia. Although the two populations interact freely, their
public health systems are under different administrations, which caused
Petrolina to have performed twice as much testing for COVID-19 proportionally
as Juazeiro, by May 13 2021, the last date we consider in our estimations.

In this example, we use the data on COVID-19 spread in Petrolina and the
deaths caused by the disease in Juazeiro to estimate the number of cumulative
infectious cases in Juazeiro. This estimate and the range of error given by
Theorem 3.4 indicates that Juazeiro’s cases might be under-reported. It is
worth mentioning that our estimations are made under assumptions (i)-(iv) in
Theorem 3.4 on the parameters and therefore, our conclusions here should not
be taken as a precise portrait of the real situation for the considered
population. Our aim in this example is to indicate a possible way of
estimating under-reporting.

We consider the data of reported cases and deaths for the period of 90 days
from February 13 to May 13 of 2021 with Group 1 as the inhabitants of
Petrolina, and Group 2 as the inhabitants of Juazeiro. The number of
individuals in each group is $\displaystyle N_{1}=353800$ and $\displaystyle
N_{2}=226200$, respectively, with total population given by $\displaystyle
N=580000$. The reader can find this dataset on Petrolina’s and Juazeiro’s city
halls websites, respectively [3] and [2]. To generate the fitted SEIRD model,
we proceeded as in the former example and, for the calculated values
$\displaystyle V_{1}=-1.052314\times 10^{-6}$ and $\displaystyle
V_{2}=3.731528\times 10^{-6}$ the fitted parameters for the model (3) can be
seen in Table 3.

Parameters
---
Symbol | Value | Symbol | Value | Symbol | Value
$\displaystyle\beta_{11}$ | $\displaystyle 3.364642\times 10^{-7}$ | $\displaystyle\mu_{1}$ | 0.003744 | $\displaystyle E_{1}(0)$ | 296.46
$\displaystyle\beta_{12}$ | $\displaystyle 3.174329\times 10^{-7}$ | $\displaystyle\mu_{2}$ | 0.002563 | $\displaystyle E_{2}(0)$ | 296.03
$\displaystyle\beta_{21}$ | $\displaystyle 3.325242\times 10^{-7}$ | $\displaystyle\gamma$ | 0.157060 |  |
$\displaystyle\beta_{22}$ | $\displaystyle 3.269971\times 10^{-7}$ | $\displaystyle\alpha$ | 0.3 |  |
Table 3. Fitted parameters for the system (3) for the data of reported cases
and deaths in Group 1, presented in [3], and reported deaths in Group 2,
presented in [2].

The Figure 4 shows the functions $\displaystyle C_{1}(t)$, $\displaystyle
C_{2}(t)$, $\displaystyle D_{1}(t)$ and $\displaystyle D_{2}(t)$ from the
system (3) for the set of parameters in Table 3 subjected to the initial
conditions

$\displaystyle\begin{matrix}S_{1}(0)=337822.5,&I_{1}(0)=207,&R_{1}(0)=15284,&D_{1}(0)=190,\\\
S_{2}(0)=217047,&I_{2}(0)=143,&R_{2}(0)=8554,&D_{2}(0)=160,\end{matrix}$

and the 7-days moving average functions of the total number of reported cases
and deaths for each group.

As shown in Figure 5, the cumulative number of reported cases for Group 2
(square markers) stays far below the correspondent curve $\displaystyle
C_{2}(t)$ and out of the estimated distance range given by the blue area. In
fact, on the last day of the period $\displaystyle t=90$ (May 13), the model
indicates that Group 2 should have a total of $\displaystyle 15770$ infectious
cases, with an error of $\displaystyle 3.26\%$ ($\displaystyle
K_{C}(90)=513.67$), while only $\displaystyle 14111$ were reported. This
suggests that the reported data for Group 2 are under-reported, or that the
hypothesis of the theorem may not be satisfied in this example, in the sense
that the estimated values in (i)-(iv) were not well approximated.

(a) Infectious cases in Group 1

(b) Deaths in Group 1

(c) Infectious cases in Group 2

(d) Deaths in Group 2

Figure 4. Fit of the SEIRD system (3) to the reported data. In (4(a)) we plot
the cumulative number of reported infectious cases (black dots) and the
function $\displaystyle C_{1}(t)$ (blue solid line) for the Group 1
(Petrolina). In (4(b)) we plot the reported deaths (black dots) and the
function $\displaystyle D_{1}(t)$ (red solid line) for the Group 1. In (4(c))
we plot the cumulative number of reported infectious cases (black dots) and
the function $\displaystyle C_{2}(t)$ (blue solid line) for the Group 2
(Juazeiro). In figure (4(d)), the reported deaths (black dots) and the
function $\displaystyle D_{2}(t)$ (red solid line) for the Group 2. Figure 5.
In this figure we plot the function cumulative number of infectious cases
$\displaystyle C_{1}(t)$ (light blue solid line) and the reported number of
infectious cases (black circle markers) for Group 1, the function cumulative
number of infectious cases $\displaystyle C_{2}(t)$ (dark blue solid line) and
reported number of infectious cases (black square markers) for Group 2, and
the estimated distance range function $\displaystyle K_{C}(t)$ (blue area).

##  5\. Conclusions and comments

In the Section 2 we proved that the final size of the susceptible individuals
is always positive for both groups. The result also shows that these values
can be obtained by a fixed-point problem in $\displaystyle\mathbb{R}^{2}$.

The results and discussions presented in Section 3 and 4 show that it is
possible to use the two group SEIRD model to estimate the dynamics of an
epidemics when the number of cases in one of the groups is not known or may be
unreliable.

In Section 4 we illustrate the results of Theorem 3.4 with two examples of the
spread of COVID-19 disease. The first example was used as a verification case,
in order to indicate that the method of estimating the number of cases is
consistent. Since the data from the New York State is more reliable, it was
expected that the number of reported cases in the New York County would fall
into the range of error given by the Theorem with respect to the calculated
function of cumulative cases. The parameters obtained for the model SEIRD led
to this result and it may be seen in Figure 3. The second application concerns
to two neighboring Brazilian cities, Petrolina and Juazeiro, which are under
different health policies with Petrolina submitting more inhabitants to
testing for COVID-19 than Juazeiro. Therefore, in this example we were
interested in checking if the proportionally low number of reported cases in
the city of Juazeiro is a potential situation of under-reporting. The
implementation of the method for fitting the parameters of the SEIRD model and
the range error estimation given by Theorem 3.4 and showed in Figure 5,
suggests that the number of infectious cases in Juazeiro may be indeed greater
than the ones reported by the health system and present an interval were the
number of cases might actually be. It worth mentioning that we are assuming
that reported data from group 1 is reliable and comes from a SEIRD model,
possibly with noise, even though it is well known that under-reporting
occurred worldwide ([20]) and that the data for group 1 may not be accurate,
specially in the Petrolina-Juazeiro case.

Since the full dataset is not available, it is not possible to be sure of the
results of the method in general. However, combined with other tools, the
method might be useful to guide public health polices, for instance when
reliable data is lacking for part of the considered population.

In Section 4 the definition of the groups is based only on geographic aspects,
but one could apply it in situations where the groups are defined by
socioeconomic classes, gender or age.

It is worth mentioning that in both examples, we are supposing that the data
follows a SEIRD model (possibly with noises) for an unknown set of parameters.
We then fitted model presented in each case with parameters which satisfy the
hypotheses of Theorem 3.4 with respect to the unknown parameters. In
particular, from the available dataset we estimate $\displaystyle
V_{1}=(\widetilde{\beta}_{21}-\widetilde{\beta}_{11})/\widetilde{\mu}_{1},\
i=1,2$, and therefore in the parameter fitting process we only consider the
class parameters for (3) that satisfy these values of $\displaystyle V_{i}$.
I.e., the light blue areas of the Figures 3 and 5 contain all solutions to
models whose parameters satisfy these values of $\displaystyle V_{i}$.

We emphasize that these estimations are here only to illustrate the
application of our theoretical results. Additional arguments would be needed
to validate our assumptions before making assertive statements on the actual
scenario of the considered populations.

Acknowledgments The authors are grateful to Sergio Floquet for the fruitful
discussions on the subject and the help with the parameter fitting process.
The authors would also like to express their sincere thanks to the anonymous
referees for their valuable comments and useful suggestions, which
significantly contributed to improving the quality of the article.

Conflict of interest The authors declare that they have no conflict of
interest.

## References

  * [1] Annual Resident Population Estimates for States and Counties - April 1, 2020 - United States Census Bureau.  https://www2.census.gov/programs-surveys/popest/datasets/2010-2020/counties/totals/co-est2020.csv.
  * [2] Boletim Epidemiologico de Juazeiro: Prefeitura Municipal de Juazeiro.  https://www6.juazeiro.ba.gov.br/category/coronavirus.
  * [3] Coronavírus Boletins Diários: Prefeitura Municipal de Petrolina.  https://petrolina.pe.gov.br/coronavirus/coronavirus-boletins-diarios.
  * [4] Report of the WHO-China Joint Mission on Coronavirus Disease 2019 (COVID-19).  https://www.who.int/docs/default-source/coronaviruse/who-china-joint-mission-on-covid-19-final-report.%****␣SEIRD-3rd_Revision.bbl␣Line␣25␣****pdf.
  * [5] USA facts.  https://usafacts.org/visualizations/coronavirus-covid-19-spread-map.
  * [6] L. Acedo, J.-A. M. no, and J. Diez-Domingo.  Cost analysis of a vaccination strategy for respiratory syncytial virus (RSV) in a network model.  Mathematical and Computer Modelling, 52(7):1016–1022, 2010.
  * [7] V. Andreasen.  The final size of an epidemic and its relation to the basic reproduction number.  Bull. Math. Biol., 73(10):2305–2321, 2011.
  * [8] D. Bernoulli.  Essai d’une nouvelle analyse de la mortalité causée par la petite verole et des avantages de l’inoculation pour la prevenir.  Mem. Math. Phys. Acad. Roy. Sci., pages 1–45, 1760.
  * [9] D. Bernoulli.  Reflexions sur les avantages de l’inoculation.  Mercure de France, June issue:173–190, 1760.
  * [10] E. Bertuzzo, L. Mari, D. Pasetto, S. Miccoli, R. Casagrandi, M. Gatto, and A. Rinaldo.  The geography of COVID-19 spread in italy and implications for the relaxation of confinement measures.  Nature Communications, 11(4264), 2020.
  * [11] F. Brauer and C. Castillo-Chavez.  Mathematical models in population biology and epidemiology, volume 40 of Texts in Applied Mathematics.  Springer, New York, second edition, 2012.
  * [12] K. Chalvet-Monfray, M. Artzrouni, J. Gouteux, P. Auger, and P. Sabatier.  A two-patch model of gambian sleeping sickness: Application to vector control strategies in a village and plantations.  Acta Biotheor, 46:207–222, 1998.
  * [13] S. Clouston, G. Natale, and B. G. Link.  Socioeconomic inequalities in the spread of coronavirus-19 in the United States: A examination of the emergence of social inequalities.  Social science & medicine (1982), 268(113554), 2021.
  * [14] O. Diekmann, H. Heesterbeek, and T. Britton.  Mathematical tools for understanding infectious disease dynamics.  Princeton Series in Theoretical and Computational Biology. Princeton University Press, Princeton, NJ, 2013.
  * [15] M. Gatto, E. Bertuzzo, L. Mari, S. Miccoli, L. Carraro, R. Casagrandi, and A. Rinaldo.  Spread and dynamics of the COVID-19 epidemic in Italy: Effects of emergency containment measures.  Proceedings of the National Academy of Sciences, 117(19):10484–10491, 2020.
  * [16] J. Graunt.  Natural and political observations made upon the bills of mortality.  1662\.
  * [17] Q. Griette and P. Magal.  Clarifying predictions for COVID-19 from testing data: The example of New York State.  Infectious Disease Modelling, 6:273–283, 2021.
  * [18] M. J. Keeling and P. Rohani.  Modeling Infectious Diseases in Humans and Animals.  Princeton University Press, 2008.
  * [19] W. O. Kermack, A. G. McKendrick, and G. T. Walker.  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, 1927.
  * [20] H. Lau, T. Khosrawipour, P. Kocbach, H. Ichii, J. Bania, and V. Khosrawipour.  Evaluating the massive underreporting and undertesting of COVID-19 cases in multiple global epicenters.  Pulmonology, 27(2):110–115, 2021.
  * [21] S. A. Lauer, K. H. Grantz, Q. Bi, F. K. Jones, H. R. Zheng, Q.and Meredith, N. G. Azman, A. S. andReich, and J. Lessler.  The incubation period of coronavirus disease 2019 (COVID-19) from publicly reported confirmed cases: Estimation and application.  Annals of internal medicine, 172(9):577–582, 2020.
  * [22] Q. Lin, S. Zhao, D. Gao, Y. Lou, S. Yang, S. S. Musa, M. H. Wang, Y. Cai, W. Wang, L. Yang, and D. He.  A conceptual model for the coronavirus disease 2019 (COVID-19) outbreak in Wuhan, China with individual reaction and governmental action.  International Journal of Infectious Diseases, 93:211–216, 2020\.
  * [23] Z. Liu, P. Magal, O. Seydi, and G. Webb.  A COVID-19 epidemic model with latency period.  Infectious Disease Modelling, 5:323–337, 2020.
  * [24] Z. Liu, P. Magal, O. Seydi, and 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), 2020.
  * [25] P. Magal, O. Seydi, and G. Webb.  Final size of an epidemic for a two-group SIR model.  SIAM J. Appl. Math., 76(5):2042–2059, 2016.
  * [26] M. Martcheva.  An introduction to mathematical epidemiology, volume 61 of Texts in Applied Mathematics.  Springer, New York, 2015.
  * [27] P. R. Martins-Filho, A. A. de Souza Araújo, L. J. Quintans-Júnior, and V. S. Santos.  COVID-19 fatality rates related to social inequality in Northeast brazil: a neighbourhood-level analysis.  Journal of travel medicine, 27(7).
  * [28] K. M. Mullen, D. Ardia, D. L. Gil, D. Windover, and J. Cline.  DEoptim: An R Package for Global Optimization by Differential evolution.  Journal of Statistical Software, Articles, 40(6):1–26, 2011.
  * [29] J. D. Murray.  Mathematical biology. I, volume 17 of Interdisciplinary Applied Mathematics.  Springer-Verlag, New York, third edition, 2002.
  * [30] D. Prodanov.  Analytical parameter estimation of the SIR epidemic model. Applications to the COVID-19 pandemic.  Entropy, 23(1):Paper No. 59, 20, 2021.
  * [31] K. Roosa, Y. Lee, R. Luo, A. Kirpich, R. Rothenberg, J. Hyman, P. Yan, and G. Chowell.  Real-time forecasts of the COVID-19 epidemic in China from February 5th to February 24th, 2020.  Infectious Disease Modelling, 5:256–263, 2020.
  * [32] R. Ross and H. P. Hudson.  An application of the theory of probabilities to the study of a priori pathometry. part ii.  Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 93(650):212–225, 1917.
  * [33] C. Rothe, M. Schunk, P. Sothmann, G. Bretzel, G. Froeschl, C. Wallrauch, T. Zimmer, V. Thiel, C. Janke, W. Guggemos, M. Seilmaier, C. Drosten, P. Vollmar, K. Zwirglmaier, S. Zange, R. Wölfel, and M. Hoelscher.  Transmission of 2019-nCoV infection from an asymptomatic contact in Germany.  New England J. of Medicine, 382(10):970–971, 2020.
  * [34] K. Soetaert, T. Petzoldt, and R. W. Setzer.  Solving differential equations in R: Package deSolve.  Journal of Statistical Software, Articles, 33(9):1–25, 2010.
  * [35] M. T. Sofonea, B. Reyné, B. Elie, R. Djidjou-Demasse, C. Selinger, Y. Michalakis, and S. Alizon.  Memory is key in capturing COVID-19 epidemiological dynamics.  Epidemics, 35:100459, 2021.
  * [36] J. Sun, X. Chen, Z. Zhang, S. Lai, B. Zhao, H. Liu, S. Wang, W. Huan, R. Zhao, M. T. A. Ng, and Y. Zheng.  Forecasting the long-term trend of COVID-19 epidemic using a dynamic model.  Scientific Reports, 10:21122, 2020.
  * [37] M. Triacca and U. Triacca.  Forecasting the number of confirmed new cases of COVID-19 in Italy for the period from 19 May to 2 June 2020.  Infectious Disease Modelling, 6:362–369, 2021.
  * [38] H. Wang, Z. Wang, Y. Dong, R. Chang, C. Xu, X. Yu, S. Zhang, L. Tsamlag, M. Shang, J. Huang, Y. Wang, S. Xu, G., X. T., Zhang, and Y. Cai.  Phase-adjusted estimation of the number of coronavirus disease 2019 cases in Wuhan, China.  Cell discovery, 6(10), 2020.