# A measure model for the spread of viral infections with mutations

<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

Genetic variations in the COVID-19 virus are one of the main causes of the
COVID-19 pandemic outbreak in 2020 and 2021. In this article, we aim to
introduce a new type of model, a system coupled with ordinary differential
equations (ODEs) and measure differential equation (MDE), stemming from the
classical SIR model for the variants distribution. Specifically, we model the
evolution of susceptible $S$ and removed $R$ populations by ODEs and the
infected $I$ population by a MDE comprised of a probability vector field (PVF)
and a source term. In addition, the ODEs for $S$ and $R$ contains terms that
are related to the measure $I$. We establish analytically the well-posedness
of the coupled ODE-MDE system by using generalized Wasserstein distance. We
give two examples to show that the proposed ODE-MDE model coincides with the
classical SIR model in case of constant or time-dependent parameters as
special cases.

###### Key words and phrases:

measure differential equation, generalized Wasserstein distance, SIR model,
viral infections, mutations.

###### 1991 Mathematics Subject Classification:

Primary: 58F15, 58F17; Secondary: 53C35.

∗ Corresponding author: Benedetto Piccoli

Xiaoqian Gong

School of Mathematical and Statistical Science

Arizona State University

Tempe, AZ, 85281, USA

Benedetto Piccoli∗

Department of Mathematical Sciences and Center for Computational and
Integrative Biology

Rutgers University

Camden, NJ, 08102, USA

##  1\. Introduction

The 2020 COVID-19 pandemic generated renewed interests in epidemiological
models for infectious diseases. Researchers and modelers from different areas
including applied math, physics, public health, engineering and others,
developed many different approaches depending on the final aim, which included
nowcasting and possible scenarios, prediction of the pandemic evolution,
evaluation of lock down and social distancing measures as well as economic
impact.

The reasons for the pandemic outbreak included difficulty in detection of the
virus for asymptomatic, fast spread due to globalization and emergence of
different variants usually associated with a specific country where they were
first observed. The latter include: B.1.1.7 initially detected in UK, B.1.351
detected in South Africa, P.1 detected in Japan for travelers from Brazil, and
B.1.427 and B.1.429 identified in California. The name variant is a convenient
way to represent a family of mutations grouped by genetic similarities. The
total number of detected mutations so far exceeds the two millions, see [12].

Our interest is in introducing a new type of models stemming from the
classical SIR model introduced in the pioneering work of Kormack and
McKendrick [11], coupled with a new type of differential equations for
measures, called Measure Differential Equations (briefly MDE) introduced in
[15] for the variants distribution among infected, which is connected to the
concept of Probability Vector Field (briefly PVF).

Let us first recall that the SIR model was generalized in number of recent
papers and in various directions, such as: 1) using time-dependent parameters
as infectivity rates [5, 19]; 2) adding population compartments for different
stages of the disease [8]; 3) increasing complexity with age-structure and
spatial models [4, 6, 7, 9, 23, 22]; 4) using multiscale approaches and
infinite dimensional systems [1, 2, 10]. The different type of models can be
used for some of the scopes discussed above [13, 20], but the difficulty in
tuning with data and use for prediction is known since long time [14].

Let us go back to our model to include virus variants in the SIR model. The
evolution of susceptible $S$ and removed $R$ populations will still be
detected by an Ordinary Differential Equation (briefly ODE). On the other
side, there are more than two million SARS-COV2 genetic variations, see the
figures on pages 844/845 in [12]. Branching out Each dot represents a virus
isolated from a COVID-19 patient in this family tree of SARS-COV2, which shows
a tiny subset of the more than 2 million viruses sequenced so far. The World
Health Organization currently recognizes four variants of concern and four
variants of interest. We assume that the virus variants are captured by a
continuous variable $\alpha\in\mathbb{R}$ and the infected population is
represented by a Radon measure $I$ on $\mathbb{R}$ with finite mass. In simple
words the value $I(A)$, $A\subset\mathbb{R}$, represents the number of
infected people having a virus variant corresponding to parameters $\alpha\in
A$. The dynamics of $I$ will be captured by an MDE with two components: a
finite-diffusion term, which represents the emergence of variants in infected
people, and a source term, which represents the inflow of susceptible getting
infected and the outflow of infected to removed. On the other side, the ODE
for $S$ and $R$ will contain terms which depend on the measure $I$
corresponding to the inflow and outflow term of the MDE. More precisely, the
original SIR dynamics is modified since the classical infection rate $\beta$
is function of the parameter $\alpha$ identifying the virus variant, and the
same occur for the recovery rate $\nu$. Therefore, the resulting dynamics is a
systems of fully coupled ODE-MDE.

To deal with the introduced model, we resort to recent results on MDEs [15]
and MDE with sources [18]. We first recall some basic tools for measures,
including the Wasserstein distance and the generalized Wasserstein distance
(since we deal with measures with variable mass). Then we provide a definition
of solution and Lipschitz semigroup of solutions for the coupled ODE-MDE
system, where the MDE is comprised of a PVF and a source term. The main result
of the paper is the existence of a Lipschitz semigroup under suitable
assumption. To state the assumptions we have to deal with a space
$\mathbb{R}^{m}\times\mathcal{M}$, where $\mathcal{M}$ indicates the space of
Radon measures with finite mass and compact support endowed with the
generalized Wasserstein distance. The conditions for the existence of the
Lipschitz semigroup are the natural generalization of the conditions for ODE
and MDE and include: sub-linear growth of supports for the MDE, Lipschitz
continuity of the vector field of the ODE and the PVF of the MDE, and,
finally, uniform boundedness and Lipschitz continuity of the source of the
MDE. Notice that the uniform boundedness of the source is chosen to simplify
some proofs, but can be relaxed by sub-linear growth conditions. A Lipschitz
semigroup can then be selected, to achieve uniqueness, by prescribing small-
time evolution of finite sums of Dirac masses. The ODE-MDE system representing
the SIR model with virus variants is then shown to satisfy the assumption for
the existence of a Lipschitz semigroup of solutions.

The paper is organized as follows: We first recall basic definitions and
results on generalized Wasserstein distance and measure differential equations
in Section $2$. Then we provide existence and uniqueness results for a system
of coupled ODE-MDE in Section $3$. In the last Section, we introduce a measure
model for viral infections with mutations consisting of one MDE for infected
coupled with a system of two differential equations for susceptible and
removed. Two examples show how the model coincides with the original SIR model
in case of constant parameters (not depending on mutation) and includes SIR
model with time-dependent parameters as a special case.

##  2\. Basic definition and results

We use the symbol $|\cdot|$ to indicate the Euclidean norm, and for every
$R>0$, $B(0,R)$ for the ball of radius $R$ centered at the origin. The symbol
$T\mathbb{R}^{n}$ indicates the tangent bundle of $\mathbb{R}^{n}$,
$\pi_{1}:T\mathbb{R}^{n}\to\mathbb{R}^{n}$ is the base-projection
$\pi_{1}(x,v)=x$ and
$\pi_{13}:(T\mathbb{R}^{n})^{2}\mapsto(\mathbb{R}^{n})^{2}$ is given by
$\pi_{13}(x,v,y,w)=(x,y)$. For every $A\subset\mathbb{R}^{n}$, $\chi_{A}$
indicates the characteristic function of the set $A$ and
${{\mathcal{C}}}^{\infty}_{c}(\mathbb{R}^{n})$ the space of compactly
supported smooth functions on $\mathbb{R}^{n}$.

Given $(X,d)$ Polish space (complete separable metric space) we indicate by
$\mathcal{M}(X)$ the set of positive Borel measures with finite mass and
compact support, by $\mathcal{P}(X)$ the set of probability measures and by
$\mathcal{P}_{c}(X)$ the set of probability measures with compact support on
$X$. For $\mu\in\mathcal{M}$ we denote by $|\mu|=\mu(X)$ the total mass and by
$Supp(\mu)$ its support. Given $(X_{1},d_{1})$, $(X_{2},d_{2})$ Polish spaces,
$\mu\in\mathcal{P}(X_{1})$, $\phi:X_{1}\to X_{2}$ measurable and Borel set
$A\subset X_{2}$, we set $\phi\\#\mu\in\mathcal{P}(X_{2})$ by
$\phi\\#\mu(A)=\mu(\phi^{-1}(A))=\mu(\\{x\in X_{1}:\phi(x)\in A\\})$. Consider
a measure $\mu\in\mathcal{M}(X_{1})$, a family of measures
$\nu_{x}\in\mathcal{M}(X_{2})$, with $x\in X_{1}$, and a function
$\varphi\colon X_{1}\times X_{2}\to\mathbb{R}$ such that $v\to\varphi(x,v)\in
L^{1}(d\nu_{x})$ for $\mu$-almost every $x$ and
$x\to\int_{X_{2}}\varphi(x,v)d\nu_{x}(v)\in L^{1}(d\mu)$. Then we define
$\int_{X_{1}\times X_{2}}\varphi(x,v)\
d(\mu\otimes\nu_{x})=\int_{X_{1}}\int_{X_{2}}\varphi(x,v)d\nu_{x}(v)\
d\mu(x)$. For $\mu$, $\nu\in\mathcal{M}(X)$, we denote by $P(\mu,\nu)$ the set
of transference plans from $\mu$ to $\nu$, i.e. the set of probability
measures on $X\times X$ with marginals equal to $\mu$ and $\nu$ respectively.
The cost of a transference plan $\tau\in P(\mu,\nu)$ is defined as
$J(\tau)=\int_{X^{2}}d(x,y)\,d\tau(x,y)$ and the Monge-Kantorovich or optimal
transport problem amounts to find a cost minimizer. The value of the attained
minimum is called the Wasserstein distance between $\mu$ and $\nu$:

$W^{X}(\mu,\nu)=\inf_{\tau\in P(\mu,\nu)}J(\tau).$

In general, a Wasserstein distance $W^{p}$ can be defined for $p\geq 1$ by
setting
$J(\tau)=\left(\int_{X^{2}}d(x,y)^{p}\,d\tau(x,y)\right)^{\frac{1}{p}}$. We
indicate by $P^{opt}(\mu,\nu)$ the set of optimal transference plans from
$\mu$ to $\nu$.

We will consider measures with time-varying total mass, thus we will consider
the generalized Wasserstein distance:

######  Definition 2.1 (The generalized Wasserstein distance).

Let $\mu,\nu\in\mathcal{M}(\mathbb{R})$ be two measures. We define the
functional

$W^{g}(\mu,\nu):=\inf_{\tilde{\mu},\tilde{\nu}\in\mathcal{M},\,|\tilde{\mu}|=|\tilde{\nu}|}|\mu-\tilde{\mu}|+|\nu-\tilde{\nu}|+W(\tilde{\mu},\tilde{\nu}).$
(1)

The generalized Wasserstein distance is thus obtained combining an $L^{1}$ or
total variation cost for removing/adding mass and transportation cost for the
rest of the mass. Various properties of the generalized Wasserstein distance
can be found in [16, 17].

Lipschitz-type conditions for evolution of measures with time-varying mass
will be defined in terms of the following operator.

######  Definition 2.2 (The operator $\mathcal{W}^{g}$).

Consider $V_{i}\in\mathcal{M}(T\mathbb{R}^{n})$ and $\tilde{V}_{i}$ satisfying
$\tilde{V}_{i}\leq V_{i},i=1,2$. Let $\mu_{i}=\pi_{1}\\#V_{i}$ and
$\tilde{\mu}_{i}=\pi_{1}\\#\tilde{V}_{i},i=1,2$ be the projections on the base
space. We define the following non-negative operator

$\displaystyle\mathcal{W}^{g}(V_{1},V_{2})\coloneqq$
$\displaystyle\inf\left\\{\int_{T\mathbb{R}^{n}\times T\mathbb{R}^{n}}|v-w|\
dT(x,v,y,w)\colon\right.$ (2) $\displaystyle\left.T\in
P(\tilde{V}_{1},\tilde{V}_{2})\text{ with },\tilde{V}_{i}\leq
V_{i},i=1,2,\pi_{13}\\#T\in P^{opt}(\tilde{\mu}_{1},\tilde{\mu}_{2})\right.$
$\displaystyle\left.\text{ where }(\tilde{\mu}_{1},\tilde{\mu}_{2})\text{ is a
minimizer in Definition }\ref{def: g-wass}\right\\}.$

###  2.1. Measure differential equations

In this section we recall the basic definitions for the theory of MDEs
introduced in [15]. The original results were provided for probability
measures but they hold for measures with constant finite mass with obvious
modifications.

These type of equations are based on a generalization of the concept of vector
field to measures as follows.

######  Definition 2.3.

A Probability Vector Field (briefly PVF) on $\mathcal{M}(\mathbb{R}^{n})$ is a
map $V:\mathcal{M}(\mathbb{R}^{n})\to\mathcal{P}(T\mathbb{R}^{n})$ such that
$\pi_{1}\\#V[\mu]=\mu$.

A Measure Differential Equation (MDE) corresponding to a PVF $V$ is defined
by:

$\dot{\mu}=V[\mu].$ (3)

The mass of $\mu$ over a set $A\subset\mathbb{R}^{n}$ is dispersed along the
velocites of the support of $V[\mu]$ restricted to $A\times\cup_{x\in
A}T_{x}\mathbb{R}^{n}$.
Given an MDE and $\mu_{0}\in\mathcal{M}(\mathbb{R}^{n})$ we define the Cauchy
problem:

$\dot{\mu}=V[\mu],\qquad\mu(0)=\mu_{0}.$ (4)

A solution to (4) is defined using weak solutions as follows.

######  Definition 2.4.

A solution to (4) is a map $\mu:[0,T]\to\mathcal{M}(\mathbb{R}^{n})$ such that
$\mu(0)=\mu_{0}$, $|\mu(t)|$ is constant and the following holds. For every
$f\in{{\mathcal{C}}}^{\infty}_{c}(\mathbb{R}^{n})$, the integral
$\int_{T\mathbb{R}^{n}}(\nabla f(x)\cdot v)\ dV[\mu(s)](x,v)$ is defined for
almost every $s$, the map $s\to\int_{T\mathbb{R}^{n}}(\nabla f(x)\cdot v)\
dV[\mu(s)](x,v)$ belongs to $L^{1}([0,T])$, and the map $t\to\int f\,d\mu(t)$
is absolutely continuous and for almost every $t\in[0,T]$ it satisfies:

$\frac{d}{dt}\int_{\mathbb{R}^{n}}f(x)\,d\mu(t)(x)=\int_{T\mathbb{R}^{n}}(\nabla
f(x)\cdot v)\ dV[\mu(t)](x,v).$ (5)

The existence of solutions to (4) holds true under the following assumptions:

  * (H:bound)

$V$ is support sublinear, i.e. there exists $C>0$ such that for every
$\mu\in\mathcal{M}(X)$ it holds:

$\sup_{(x,v)\in Supp(V[\mu])}|v|\leq C\left(1+\sup_{x\in
Supp(\mu)}|x|\right).$

  * (H:cont)

the map $V:\mathcal{M}(\mathbb{R}^{n})\to\mathcal{M}(T\mathbb{R}^{n})$ is
continuous (for the topology given by the Wasserstein metrics
$W^{\mathbb{R}^{n}}$ and $W^{T\mathbb{R}^{n}}$.)

The existence of solutions is achieved via approximate solutions called
Lattice Approximate Solutions (briefly LAS). To define in details LAS, we need
some additional notation.

For $N\in{\mathbb{N}}$ let $\Delta_{N}=\frac{1}{N}$ be the time step size,
$\Delta^{v}_{N}=\frac{1}{N}$ the velocity step size and
$\Delta^{x}_{N}=\Delta^{v}_{N}\Delta_{N}=\frac{1}{N^{2}}$ the space step size.
The discretization in position and velocity are given by the points $x_{i}$ of
${\mathbb{Z}}^{n}/(N^{2})\cap[-N,N]^{n}$ and $v_{j}$ of
${\mathbb{Z}}^{n}/N\cap[-N,N]^{n}$. Every $\mu\in\mathcal{M}(\mathbb{R}^{n})$
can be approximated by Dirac deltas using the operator:

${\mathcal{A}}^{x}_{N}(\mu)=\sum_{i}m^{x}_{i}(\mu)\delta_{x_{i}}$ (6)

with $m^{x}_{i}(\mu)=\mu(x_{i}+Q)$ amd $Q=([0,\frac{1}{N^{2}}[)^{n}$. The PVF
can also be approximated using the operator:

${\mathcal{A}}^{v}_{N}(V[\mu])=\sum_{i}\sum_{j}m^{v}_{ij}(V[\mu])\
\delta_{(x_{i},v_{j})}$ (7)

where $m^{v}_{ij}(V[\mu])=V[\mu](\\{(x_{i},v):v\in v_{j}+Q^{\prime}\\})$, and
$Q^{\prime}=([0,\frac{1}{N}[)^{n}$.

The two operators are good approximations for the Wasserstein distances in the
following sense:

######  Lemma 2.5.

Given $\mu\in\mathcal{P}_{c}(\mathbb{R}^{n})$, for $N$ sufficiently big the
following holds:

$W({\mathcal{A}}^{x}_{N}(\mu),\mu)\leq{\sqrt{n}}\,\Delta^{x}_{N},\qquad
W^{T\mathbb{R}^{n}}({\mathcal{A}}^{v}_{N}(V[\mu]),V[\mu])\leq{\sqrt{n}}\,\Delta^{v}_{N}.$

The definition of LAS is as follows.

######  Definition 2.6.

Consider $V$ satisfying (H:bound), (4), $T>0$ and $N\in{\mathbb{N}}$. The LAS
$\mu^{N}:[0,T]\to\mathcal{P}_{c}(\mathbb{R}^{n})$ is defined by recursion.
First $\mu_{0}^{N}={\mathcal{A}}^{x}_{N}(\mu_{0})$ and then:

$\mu^{N}_{\ell+1}=\mu^{N}((\ell+1)\Delta_{N})=\sum_{i}\sum_{j}m^{v}_{ij}(V[\mu^{N}(\ell\Delta_{N})])\
\delta_{x_{i}+\Delta_{N}\,v_{j}}.$ (8)

Notice that $Supp(\mu^{N}_{\ell})$ is contained in the set
${\mathbb{Z}}^{n}/(N^{2})\cap[-N,N]^{n}$, thus
$\mu^{N}_{\ell}=\sum_{i}m^{N,\ell}_{i}\delta_{x_{i}}$ for some
$m^{N,\ell}_{i}\geq 0$. We can define LAS for all times by interpolation:

$\mu^{N}(\ell\Delta_{N}+t)=\sum_{ij}m^{v}_{ij}(V[\mu^{N}(\ell\Delta_{N})])\
\delta_{x_{i}+t\,v_{j}}.$ (9)

It is easy to prove uniform bounds for LAS ([15]):

######  Lemma 2.7.

Given a PVF $V$ satisfying (H:bound), $\mu_{0}$ with $Supp(\mu_{0})\subset
B(0,R)$ and $\ell$ such that $\ell\Delta_{N}\leq T$, the following holds true:

$Supp(\mu^{N}_{\ell})\subset B\left(0,e^{C_{N}T}(R_{N}+1)-1\right),$ (10)

where $C_{N}=C+\frac{\sqrt{n}}{N}$ and $R_{N}=R+\frac{\sqrt{n}}{N^{2}}$.

###  2.2. Finite speed diffusion via MDEs

MDEs can be used to model finite speed diffusion. Let us first illustrate a
simple example of mass splitting, which is also related to the Wasserstein
gradient flow with interaction energy
$\Phi(\mu):=-\int_{\mathbb{R}\times\mathbb{R}}|x-y|d\mu(x)d\mu(y)$, see [3].

######  Example 2.8.

Given $\mu\in\mathcal{M}(\mathbb{R})$ we denote by $B(\mu)$ the barycenter of
$\mu$ (i.e. $B(\mu)=\sup\left\\{x:\mu(]-\infty,x])\leq\frac{1}{2}\right\\}$)
and set $V[\mu]=\mu\otimes\nu_{x}$, with:

$\nu_{x}=\left\\{\begin{array}[]{ll}\delta_{-1}&\textrm{if}\ x<B(\mu)\\\
\delta_{1}&\textrm{if}\ x>B(\mu)\\\
{\frac{1}{\mu(\\{B(\mu)\\})}\left(\eta\delta_{1}+\left(\frac{1}{2}-\mu(]-\infty,B(\mu)[)\right)\delta_{-1}\right)}&\textrm{if}\
x=B(\mu),\mu(\\{B(\mu)\\})>0\end{array}\right.$ (11)

where $\eta=\mu(]-\infty,B(\mu)])-\frac{1}{2}$.

We have the following:

######  Proposition 2.9.

The LAS for the PVF (11) converges to:

$\mu(t)(A)=\mu_{0}((A\cap]-\infty,B(\mu)-t[)+t)+\mu_{0}((A\cap]B(\mu)+t,+\infty[)-t)$
$+{\frac{1}{\mu_{0}(\\{B(\mu_{0})\\})}\left(\eta\delta_{B(\mu_{0})+t}(A)+(\frac{1}{2}-\mu_{0}(]-\infty,B(\mu_{0})[))\delta_{B(\mu_{0})-t}(A)\right).}$

Notice that from Proposition 2.9 we deduce that for the initial datum
$\mu_{0}=\delta_{x_{0}}$ the solution split mass in half, i.e.
$\mu(t)=\frac{1}{2}\delta_{x_{0}+t}+\frac{1}{2}\delta_{x_{0}-t}$.

Example 2.8 can be generalized as follows.

######  Example 2.10.

Consider an increasing map $\varphi:[0,1]\to\mathbb{R}$ and define
$V_{\varphi}[\mu]=\mu\otimes_{x}J_{\varphi}(x)$, where

$J_{\varphi}(x)=\begin{cases}\delta_{\varphi(F_{\mu}(x))}&\mbox{~{}~{}if~{}}F_{\mu}(x^{-})=F_{\mu}(x),\\\
\frac{\varphi\\#\left(\chi_{[F_{\mu}(x^{-}),F_{\mu}(x)]}\lambda\right)}{F_{\mu}(x)-F_{\mu}(x^{-})}&\mbox{~{}~{}otherwise,}\end{cases}$

where $F_{\mu}(x)=\mu(]-\infty,x])$, the cumulative distribution of $\mu$, and
$\lambda$ the Lebesgue measure. In simple words $V[\mu]$ moves the ordered
masses with speed prescribed by $\varphi$.
If $\varphi$ is a diffeomorphism, the solution from $\delta_{0}$ is given by
$g(t,x)\lambda$ with

$g(t,x)=\frac{1}{t\varphi^{\prime}(\varphi^{-1}(\frac{x}{t}))}=\frac{(\varphi^{-1})^{\prime}(\frac{x}{t})}{t}.$

For $\varphi(\alpha)=\alpha-\frac{1}{2}$ we get
$g(t,x)=\frac{1}{t}\chi_{[-\frac{t}{2},\frac{t}{2}]}$. In general
$V_{\varphi}$ gives rise to any $g$, which is solution to the equation
$g_{t}+\frac{x}{t}g_{x}=0$.

##  3\. General theory for coupled ODE-MDE

In this section we provide existence and uniqueness results for systems
consisting of a MDE coupled with an ODE.

######  Definition 3.1.

A coupled ODE-MDE system is a system written as:

$\left\\{\begin{array}[]{c}\dot{x}=g(x,\mu)\\\
\dot{\mu}=V[\mu]+s(\mu,x)\end{array}\right.$ (12)

where $g:\mathbb{R}^{m}\times\mathcal{M}(\mathbb{R}^{n})\to\mathbb{R}^{m}$,
$V:\mathcal{M}(\mathbb{R}^{n})\to\mathcal{M}(T\mathbb{R}^{n})$ is a PVF on
$\mathcal{M}(\mathbb{R}^{n})$ and
$s:\mathcal{M}(\mathbb{R}^{n})\times\mathbb{R}^{m}\to\mathcal{M}(\mathbb{R}^{n})$.

######  Definition 3.2.

A solution to (12) with given initial datum
$(x_{0},\mu_{0})\in\mathbb{R}^{m}\times\mathcal{M}(\mathbb{R}^{n})$ is a
couple $(x,\mu)$, with $x:[0,T]\to\mathbb{R}^{m}$ and
$\mu:[0,T]\to\mathcal{M}(\mathbb{R}^{n})$ such that for each
$f\in{{\mathcal{C}}}^{\infty}_{c}(\mathbb{R}^{n})$ the following holds:

  * i)

The map $t\to x(t)$ is absolutely continuous and satisfies
$x(t)=x_{0}+\int_{0}^{t}g(x(\tau),\mu(\tau))\,d\tau$ for almost every
$t\in[0,T]$.

  * ii)

$\mu(0)=\mu_{0}$;

  * iii)

There exists $C(T)>0$ such that $|\mu(t)|\leq C(T)$;

  * iv)

The integral $\int_{T\mathbb{R}^{n}}(\nabla f(y)\cdot v)\,dV[\mu(t)](y,v)$ is
defined for almost every $t\in[0,T]$;

  * v)

The map $t\to\int_{\mathbb{R}^{n}}f(y)\,ds[\mu(t),x(t)](y)$ belongs to
$L^{1}([0,T])$;

  * vi)

The map $t\to\int_{\mathbb{R}^{n}}f(y)\,d\mu(t)(y)$ is absolutely continuous
and for almost every $t\in[0,T]$ it satisfies:

$\frac{d}{dt}\int_{\mathbb{R}^{n}}f(y)\,d\mu(t)(y)=\int_{T\mathbb{R}^{n}}(\nabla
f(y)\cdot v)\ dV[\mu(t)](y,v)+\int_{\mathbb{R}^{n}}f(y)\,ds[\mu(t),x(t)](y).$
(13)

######  Definition 3.3.

A Lipschitz semigroup for (12) is a map
$S:[0,T]\times\mathbb{R}^{m}\times\mathcal{M}(\mathbb{R}^{n})\to\mathbb{R}^{m}\times\mathcal{M}(\mathbb{R}^{n})$
such that for every
$(x,\mu),(y,\nu)\in\mathbb{R}^{m}\times\mathcal{M}(\mathbb{R}^{n})$ and
$s,t\in[0,T]$ the following holds:

  * i)

$S_{0}(x,\mu)=(x,\mu)$ and $S_{t}\,S_{s}\,(x,\mu)=S_{t+s}\,(x,\mu)$;

  * ii)

The map $t\mapsto S_{t}(x,\mu)$ is a solution to (12) in the sense of
Definition 3.2;

  * iii)

Denote by $S^{i}_{t}$, $i=1,2$, the components of $S_{t}$, so that
$S_{t}=(S_{t}^{1},S_{t}^{2})$ with
$S^{1}_{t}\colon\mathbb{R}^{m}\times\mathcal{M}(\mathbb{R}^{n})\to\mathbb{R}^{m}$
and
$S^{2}_{t}\colon\mathbb{R}^{m}\times\mathcal{M}(\mathbb{R}^{n})\to\mathcal{M}(\mathbb{R}^{n})$.
For every $R,M>0$ there exists $C=C(R,M)>0$ such that if $|x|,|y|\leq M$
$Supp(\mu)\cup Supp(\nu)\subset B(0,R)$ and
$\mu(\mathbb{R}^{n})+\nu(\mathbb{R}^{n})\leq M$ then we have for every
$s,t\in[0,T]$:

$|S^{1}_{t}(x,\mu)|\leq C;$ (14) $Supp(S^{2}_{t}(x,\mu))\subset
B(0,e^{Ct}(R+M+1));$ (15)
$|S^{1}_{t}(x,\mu)-S^{1}_{t}(y,\nu)|+W^{g}(S^{2}_{t}(x,\mu),S^{2}_{t}(y,\nu))\leq
e^{Ct}\left(|x-y|+W^{g}(\mu,\nu)\right);$ (16)
$|S^{1}_{t}(x,\mu)-S^{1}_{s}(x,\mu)|+W^{g}(S^{2}_{t}(x,\mu),S^{2}_{s}(x,\mu))\leq
C\ |t-s|.$ (17)

We first state our assumptions on (12) for existence of a semigroup of
solutions.

######  Theorem 3.4 (Existence of Lipschitz semigroup of solutions to (12)).

Consider the ODE-MDE system (12) and assume the following. The PVF $V$
satisfies the assumption (H:bound) and:

(OM:Lip-g):

for each $R>0$ there exists $L=L(R)>0$ such that if $|x|,|y|\leq R$ and
$\mathrm{supp}(\mu)\cup\mathrm{supp}(\nu)\subset B(0,R)$ then:

$|g(x,\mu)-g(y,\nu)|\leq L\ \left(|x-y|+W^{g}(\mu,\nu)\right);$ (18)

(OM:Lip-V):

for each $R>0$ there exists $K=K(R)>0$ such that if
$\mathrm{supp}(\mu)\cup\mathrm{supp}(\nu)\subset B(0,R)$ then:

$\mathcal{W}^{g}(V[\mu],V[\nu])\leq KW^{g}(\mu,\nu);$ (19)

(OM:bound-s):

there exists $\bar{R}$ such that for all $x\in\mathbb{R}^{m}$ and
$\mu\in\mathcal{M}(\mathbb{R}^{n})$ it holds $\mathrm{supp}(s[\mu,x])\subset
B(0,\bar{R})$ and $|s[\mu,x]|\leq\bar{R}$.

(OM:Lip-s):

there exists $M$ such that for all $x,y\in\mathbb{R}^{m}$ and
$\mu,\nu\in\mathcal{M}(\mathbb{R}^{n})$ it holds

$W^{g}(s[\mu,x],s[\nu,y])\leq M\left(|x-y|+W^{g}(\mu,\nu)\right).$ (20)

Then, there exists a Lipschitz semigroup of solutions to (12).

To prove Theorem 3.4, we define an Euler-LAS scheme to construct approximate
solutions. The idea is to use the standard Euler explicit scheme for the ODE
and the LAS for the MDE. The details are as follows.

######  Definition 3.5.

An Euler-LAS approximate solution to (12) is obtained as follows. For fixed
$N$, set $\Delta_{N}=\frac{1}{N}$ and perform the following steps.

  * Step 1

Initial step: Define $x_{0}^{N}=x_{0}$ and
$\mu_{0}^{N}={\mathcal{A}}^{x}_{N}(\mu_{0})=\sum\limits_{i}m_{i}^{x}(\mu_{0})\delta_{x_{i}}$;

  * Step 2

Inductive Step: Define

$\mu^{N}_{\ell+1}=\sum_{i}\sum_{j}m^{v}_{ij}(V[\mu^{N}_{\ell}])\
\delta_{x_{i}+\Delta_{N}\,v_{j}}+\Delta_{N}{\mathcal{A}}^{x}_{N}(s[\mu^{N}_{\ell},x^{N}_{\ell}]),$
$x^{N}_{\ell+1}=x^{N}_{\ell}+\Delta_{N}g(x^{N}_{\ell},\mu^{N}_{\ell});$

  * Step 3

Interpolated measure: For the intermediate times $\tau\in(0,\Delta_{N})$,

$\mu^{N}_{\ell}(\tau)=\sum_{i}\sum_{j}m^{v}_{ij}(V[\mu^{N}_{\ell}])\
\delta_{x_{i}+\tau\,v_{j}}+\tau{\mathcal{A}}^{x}_{N}(s[\mu^{N}_{\ell},x^{N}_{\ell}]),$
(21) $x^{N}_{\ell}(\tau)=x^{N}_{\ell}+\tau g(x^{N}_{\ell},\mu^{N}_{\ell}).$
(22)

######  Lemma 3.6 (Uniform boundedness).

Under the assumption in Theorem 3.4, for every $N,\ell\in\mathbb{N}^{+}$ the
measure $\mu^{N}_{\ell}$ has uniformly bounded mass and support and
$x^{N}_{\ell}$ is uniformly bounded in norm, where $\mu^{N}_{\ell}$ and
$x_{\ell}^{N}$ are defined as in Definition 3.5.

###### Proof.

First fix a time $T>0$ and initial conditions
$(x_{0},\mu_{0})\in\mathbb{R}^{m}\times\mathcal{M}(\mathbb{R}^{n})$ for the
system (12). Choose $\bar{R}>0$ in (OM:bound-s) being an upper bound of the
maximal support and the total mass of $s[x,\mu]$ for all $x\in\mathbb{R}^{m}$
and $\mu\in\mathcal{M}(\mathbb{R}^{n})$. One can also enlarging $\bar{R}$ such
that $\mathrm{supp}(\mu_{0})\subset B(0,\bar{R})$. Let
$(x^{N}_{\ell},\mu^{N}_{\ell})$ be the sequence of Euler-LAS approximate
solutions defined in Definition 3.5, $N\in\mathbb{N}^{+}$.

First notice that $\mathrm{supp}(\mu_{0})\subset B(0,\bar{R})$ implies that
$\mathrm{supp}(\mu^{N}_{0})\subset B(0,\bar{R}+1)$. In addition, for all
$x\in\mathbb{R}^{m},\mu\in\mathcal{M}(\mathbb{R}^{n})$,
$\mathrm{supp}(s[\mu,x])\subset B(0,\bar{R})$ implies that
$\mathrm{supp}(\mathcal{A}_{N}^{x}(s[\mu,x]))\subset B(0,\bar{R}+1)$.

Furthermore, from (H:bound) we get: for each term $i,j$ in equation (21) it
holds that $(x_{i},v_{j})\in\mathrm{supp}(V[\mu_{\ell}^{N}])$ implies that
$|v_{j}|\leq C(1+\sup\limits_{x\in Supp(\mu^{N}_{\ell})}|x|)$. Therefore,

$\sup_{x\in Supp(\mu^{N}_{\ell}(\tau))}|x|\leq\max\\{\bar{R},\sup_{x\in
Supp(\mu^{N}_{\ell})}|x|\\}+\Delta_{N}\cdot C\left(1+\sup_{x\in
Supp(\mu^{N}_{\ell})}|x|\right),$

Note that $\ell\leq\frac{T}{\Delta N}+1$ and by induction one can derive that
$Supp(\mu^{N}_{\ell})$ is uniformly bounded by

$C_{1}=e^{CT}\max\\{\bar{R},\sup_{x\in Supp(\mu_{0})}|x|\\}.$

Now observe that:

$\displaystyle|\mu^{N}_{\ell+1}|\leq|\mu^{N}_{\ell}|+\Delta_{N}\bar{R},$ (23)
$\displaystyle W^{g}(\mu^{N}_{\ell+1},\mu_{0})$ $\displaystyle\leq
W^{g}(\mu^{N}_{\ell},\mu_{0})+W^{g}(\mu^{N}_{\ell+1},\mu_{\ell}^{N})$ (24)
$\displaystyle\leq
W^{g}(\mu^{N}_{\ell},\mu_{0})+\Delta_{N}(1+C_{1})|\mu^{N}_{\ell}|+\Delta_{N}\bar{R},$

and

$\displaystyle|x^{N}_{\ell+1}-x^{N}_{\ell}|$
$\displaystyle=\left|\Delta_{N}g(x^{N}_{\ell},\mu^{N}_{\ell})\right|$ (25)
$\displaystyle\leq\Delta_{N}\left(|g(x_{0},\mu_{\ell}^{N})|+L|x_{\ell}^{N}-x_{0}|\right)$
$\displaystyle\leq\Delta_{N}\left(|g(x_{0},\mu_{0})|+W^{g}(\mu_{\ell}^{N},\mu_{0})+L|x_{\ell}^{N}-x_{0}|\right).$

From inequality (23) we get that the mass of $\mu^{N}_{\ell}$ is uniformly
bounded. This implies, using inequality (24) that
$W^{g}(\mu^{N}_{\ell},\mu_{0})$ is also bounded. Finally, inequality (25)
implies that $x^{N}_{\ell}$ is bounded in $\mathbb{R}^{m}$. ∎

######  Lemma 3.7 (Existence of weak solution to system (12)).

Under the assumption in Theorem 3.4, there exists a limit curve
$(x,\mu)\colon[0,T]\to\mathbb{R}^{m}\times\mathcal{M}(\mathbb{R}^{n})$ such
that $\mu^{N}\rightharpoonup\mu$ and $x^{N}\to x$ as $N\to\infty$, where
$\mu^{N}$ and $x^{N}$ are defined as in Definition 3.5. Moreover, the limit
curve $(x,\mu)$ is a weak solution to system (12) as defined Definition 3.2.

###### Proof.

By Lemma 3.6, we have that $\mu^{N}_{\ell}$ is pre-compact in
$\mathcal{M}(\mathbb{R}^{n})$ and $x^{N}_{\ell}$ is pre-compact in
$\mathbb{R}^{m}$. By a diagonal argument, we can define a limit curve
$(x,\mu)(t)$ defined on $[0,T]$. We claim that the limit curve
$(x,\mu)\colon[0,T]\to\mathbb{R}^{m}\times\mathcal{M}(\mathbb{R}^{n})$ is a
weak solution to system (12) as defined Definition 3.2.

First observe that $W^{g}(\mu_{0}^{N},\mu_{0})\leq\Delta_{N}^{x}|\mu_{0}|\to
0$ as $N\to\infty$. Thus $\mu(0)=\mu_{0}$. By the definition of $x_{0}^{N}$,
it is clear that $x(0)=x_{0}$. By the definition of $x^{N}$ and the fact that
$\mu^{N}\rightharpoonup\mu$ and $x^{N}\to x$ as $N\to\infty$, we have,
$x(t)=x_{0}+\int_{0}^{t}g(x(\tau),\mu(\tau))\,d\tau$ for a.e. $t\in[0,T]$.

For fixed $f\in\mathcal{C}_{c}^{\infty}(\mathbb{R}^{n})$, and
$\sigma,\tau\in[0,T]$ with $\sigma>\tau$, define the operator $F^{N}$ as
follows:

$\displaystyle
F^{N}(\sigma,\tau)\coloneqq\int_{\mathbb{R}^{n}}f(y)\,d(\mu^{N}(\sigma)-\mu^{N}(\tau))(y)$
$\displaystyle-\int_{\tau}^{\sigma}\,dt\left(\int_{T\mathbb{R}^{n}}(\nabla
f(y)\cdot v)\
dV[\mu^{N}(t)](y,v)+\int_{\mathbb{R}^{n}}f(y)\,ds[\mu^{N}(t),x^{N}(t)](y)\right).$

Definition 3.5 implies that for
$\sigma,\tau\in[\ell\Delta_{N},(\ell+1)\Delta_{N}]$,

$\displaystyle\int_{\mathbb{R}^{n}}f(y)\,d(\mu^{N}(\sigma)-\mu^{N}(\tau))(y)=$
$\displaystyle=$
$\displaystyle\int_{\mathbb{R}^{n}}f(y)\,\left(\sum\limits_{i,j}m_{ij}^{v}(V[\mu_{\ell}^{N}])\,d(\delta_{x_{i}+(\sigma-\ell\Delta_{N})v_{j}}-\delta_{x_{i}+(\tau-\ell\Delta_{N})v_{j}})(y)+\right.$
$\displaystyle+\left.(\sigma-\tau)d\mathcal{A}_{N}^{x}(s[\mu_{\ell}^{N},x_{\ell}^{N}])(y)\right)$

Define $g_{ij}(\alpha)\coloneqq
f(x_{i}+\alpha(\sigma-\ell\Delta_{N})v_{j})-f(x_{i}+\alpha(\tau-\ell\Delta_{N})v_{j})$.
Then
$g_{ij}(1)=\int_{\mathbb{R}^{n}}f(y)\,d(\delta_{x_{i}+(\sigma-\ell\Delta_{N})v_{j}}-\delta_{x_{i}+(\tau-\ell\Delta_{N})v_{j}})(y)$.
In addition, one can bound $|F^{N}(\sigma,\tau)|$ from above by

$\displaystyle\left|\sum\limits_{i,j}g_{ij}(1)m_{ij}^{v}(V[\mu_{\ell}^{N}])-\int_{\tau}^{\sigma}\,dt\int_{T\mathbb{R}^{n}}(\nabla
f(y)\cdot v)\ dV[\mu_{\ell}^{N}(t)](y,v)\right|+$ (26)
$\displaystyle\left|(\sigma-\tau)\int_{\mathbb{R}^{n}}f(y)d\mathcal{A}_{N}^{x}(s[\mu_{\ell}^{N},x_{\ell}^{N}])(y)-\int_{\tau}^{\sigma}\,dt\int_{\mathbb{R}^{n}}f(y)\,ds[\mu_{\ell}^{N}(t),x_{\ell}^{N}(t)](y)\right|$
$\displaystyle\leq$
$\displaystyle|\tau-\sigma|\Delta_{N}\|f\|_{C^{3}}\tilde{C}$

for a suitable constant $\tilde{C}$. For more details, please see [18].

For a general pair $\sigma,\tau\in[0,T]$ with $\tau<\sigma$, assume that $N$
is sufficiently large such that

$(k_{1}-1)\Delta_{N}\leq\tau<k_{1}\Delta_{N}\leq
k_{2}\Delta_{N}<\sigma\leq(k_{2}+1)\Delta_{N}$

for some $k_{1},k_{2}\in\mathbb{N}^{+}\setminus\\{0\\}$. Consider the
convergent subsequence $\mu^{N}\rightharpoonup\mu$, and define

$\displaystyle\mathcal{L}=$
$\displaystyle\lim\limits_{N\to\infty}\frac{1}{\sigma-\tau}F^{N}(\sigma,\tau)=\frac{1}{\sigma-\tau}\left(\int
f(y)\,d(\mu(\sigma)-\mu(\tau))\right.$
$\displaystyle\left.-\int_{\tau}^{\sigma}\,dt\int_{T\mathbb{R}^{n}}(\nabla
f(y)\cdot v)\
dV[\mu(t)](y,v)-\int_{\tau}^{\sigma}\int_{\mathbb{R}^{n}}f(y)\,ds[\mu(t),x(t)](y)\right).$

The second identity follows from the continuity of both $V$ and $s$ with
respect to weak convergence of measures. We claim that
$\lim\limits_{\sigma\to\tau}|\mathcal{L}|=0$. In fact, we have

$F^{N}(\sigma,\tau)=F^{N}(\tau,k_{1}\Delta_{N})+\sum\limits_{k=k_{1}}^{k_{2}-1}F^{N}(k\Delta_{N},(k+1)\Delta_{N})+F^{N}(k_{2}\Delta_{N},\sigma).$
(27)

Combining equations (26) and (27), we get

$\displaystyle\lim\limits_{\sigma\to\tau}|\mathcal{L}|$
$\displaystyle=\lim\limits_{\sigma\to\tau}\left|\lim\limits_{N\to\infty}\frac{1}{\sigma-\tau}F^{N}(\sigma,\tau)\right|$
$\displaystyle=\lim\limits_{\sigma\to\tau}\frac{1}{\tau-\sigma}\lim\limits_{N\to\infty}|F^{N}(\sigma,\tau)|$
$\displaystyle\leq\lim\limits_{\sigma\to\tau}\frac{1}{\tau-\sigma}\lim\limits_{N\to\infty}(|\tau-
k_{1}\Delta_{N}|+\dots+|k_{2}\Delta_{N}-\sigma|)\Delta_{N}\|f\|_{C^{3}}\tilde{C}$
$\displaystyle=0.$

∎

Now we are ready to prove Theorem 3.4.

###### Proof.

We only need to verify equations (16) and (17) in Definition 3.3. Choose two
different sets of initial datum $(x_{0},\mu_{0})$ and $(y_{0},\nu_{0})$ in
$\mathbb{R}^{m}\times\mathcal{M}(\mathbb{R}^{n})$, and build Euler-LAS
approximate solutions $(x^{N},\mu^{N})$ and $(y^{N},\nu^{N})$ to (12)
according to Definition 3.5.

Thanks to the uniform boundedness of the supports of $\mu^{N}$, $\nu^{N}$,
$V[\mu^{N}]$, $V[\nu^{N}]$ and the fact that
$\mathcal{A}_{N}^{x}(\mu_{0}^{N})=\mu^{N}_{0}$, we can assume that $N$ is
sufficiently large such that for all $\ell\in\mathbb{N}$ with
$\ell\leq\frac{T}{\Delta_{N}}$the followings are true:

$\displaystyle\mathcal{A}_{N}^{x}(\mu_{\ell}^{N})$
$\displaystyle=\mu^{N}_{\ell}$
$\displaystyle\mathcal{A}_{N}^{x}(\nu_{\ell}^{N})$
$\displaystyle=\nu^{N}_{\ell}$
$\displaystyle\pi_{1}\\#\mathcal{A}_{N}^{v}(V[\mu_{\ell}^{N}])$
$\displaystyle=\pi_{1}\\#V[\mu_{\ell}^{N}]$
$\displaystyle\pi_{1}\\#\mathcal{A}_{N}^{v}(V[\nu_{\ell}^{N}])$
$\displaystyle=\pi_{1}\\#V[\nu_{\ell}^{N}].$

Let $V_{1}\coloneqq\mathcal{A}_{N}^{v}(V[\mu^{N}_{\ell}])$ and
$V_{2}\coloneqq\mathcal{A}_{N}^{v}(V[\nu^{N}_{\ell}])$. By definitions of
$\mu^{N}$, $\nu^{N}$, the generalized Wasserstein distance $W^{g}$ and the
operator $\mathcal{W}^{g}$, one can estimate recursively as

$W^{g}(\mu^{N}_{\ell+1},\nu^{N}_{\ell+1})\leq
W^{g}(\mu^{N}_{\ell},\nu^{N}_{\ell})+\Delta_{N}\mathcal{W}^{g}(V_{1},V_{2}).$
(28)

Furthermore, we can compare $\mathcal{W}^{g}(V_{1},V_{2})$ with
$\mathcal{W}^{g}(V[\mu^{N}_{\ell}],V[\nu^{N}_{\ell})]$ as

$\mathcal{W}^{g}(V_{1},V_{2})\leq\mathcal{W}^{g}(V[\mu^{N}_{\ell}],V[\nu^{N}_{\ell})]+2\sqrt{n}\Delta_{N}|\mu^{N}_{\ell}|.$
(29)

For more details, please see [18].

Now combine equations (28) and (29), and use the hypothesis (OM:Lip-V) and
Lemma 3.6, we have, there exists $\bar{C}>0$, such that

$W^{g}(\mu^{N}_{\ell+1},\nu^{N}_{\ell+1})\leq(1+K)W^{g}(\mu^{N}_{\ell},\nu^{N}_{\ell})+2\sqrt{n}(\Delta_{N})^{2}\bar{C}.$
(30)

By induction and Lemma 2.5, equation (30) implies that

$W^{g}(\mu^{N}(t),\nu^{N}(t))\leq
e^{Kt}W^{g}(\mu_{0},\nu_{0})+2\sqrt{n}(\Delta_{N})^{2}\bar{C}\frac{e^{Kt}-1}{K}.$
(31)

Now we pass to the limit by using a density argument. Let

$\mathcal{D}=\left\\{\mu_{0}\in\mathcal{M},\text{ s.t.
}\mu_{0}=\sum\limits_{i}m_{i}\delta_{x_{i}},m_{i}\in\mathbb{Q}^{+},x_{i}\in\mathbb{Q}^{n}\right\\}.$

Choose $(x_{0},\mu_{0}),(y_{0},\nu_{0})\in\mathbb{R}^{m}\times\mathcal{D}$,
define the corresponding sequences $(x^{N},\mu^{N})$ and $(y^{N},\nu^{N})$
according to Definition 3.5 such that there exist subsequences
$(x^{N_{k}},\mu^{N_{k}})$ and $(y^{N_{k}},\nu^{N_{k}})$ converges to solutions
$(x,\mu)$ and $(y,\nu)$ to system (12). Repeat this diagonal argument for
initial data in $\mathcal{D}$ and pass to the limit in inequality (31) for
$N\to\infty$, we have

$W^{g}(\mu(t),\nu(t))\leq e^{Kt}W^{g}(\mu_{0},\nu_{0})\text{ for all
}\mu_{0},\nu_{0}\in\mathcal{D}.$

Since the set $\mathcal{D}$ is countable and dense in $\mathcal{M}$, the
Lipschitz continuity with respect to the initial condition can be extended to
$\mathcal{M}$. Finally, using hypothesis (OM:Lip-g), we obtain (16).

To prove equation (17), the uniform Lipschitz continuity of the solution to
system (12) with respect to time, it is enough to notice that the sequence
$\mu^{N}$ are equi-Lipschitz in time with respect to the distance $W^{g}$. For
more details we refer the reader to [18]. ∎

Next we will provide the uniqueness result for solutions to system (12). Due
to the weak concept of solution, uniqueness can be achieved only at the level
of semigroup by prescribing small-time evolution of finite sums of Dirac
masses as in [15]. First of all, we recall the definitions of Dirac germs
compatible with a given PVF and semigroup compatible with a given Dirac germ.

######  Definition 3.8 (Dirac germ compatible with a PVF).

Let $V$ be a given fixed PVF and define
$\mathcal{M}^{D}\coloneqq\\{\mu\in\mathcal{M}(\mathbb{R}^{n})\text{ such that
}\mu=\sum\limits_{l=1}^{m}m_{l}\delta_{x_{l}},m_{l}\in\mathbb{Q}^{+},x_{l}\in\mathbb{Q}^{n}\\}$
the space of discrete measures. A Dirac germ $\gamma$ compatible with $V$ is a
map that assigns to each $\mu\in\mathcal{M}^{D}$ a Lipschitz curve
$\gamma_{\mu}\colon[0,\varepsilon(\mu)]\to\mathcal{M}(\mathbb{R}^{n})$ such
that

  * (i)

$\varepsilon[\mu]>0$ is uniformly positive for measures with uniformly bounded
support.

  * (ii)

$\gamma_{\mu}$ is a solution to equation (3) as defined in Definition 2.4.

######  Definition 3.9 (Compatibility of a semigroup for system (12)).

Fix a given PVF $V$ that satisfies (H: bound), a final time $T>0$, and a Dirac
germ $\gamma$ that is compatible with $V$ as in Definition 3.8. We say that a
semigroup
$S\colon[0,T]\times\mathbb{R}^{m}\times\mathcal{M}(\mathbb{R}^{n})\to\mathbb{R}^{m}\times\mathcal{M}(\mathbb{R}^{n})$
for system (12) is compatible with $\gamma$ if for each $R,M>0$ there exists
$C(R,M)$ such that for fixed $x_{0}\in\mathbb{R}^{m}$ with $|x_{0}|\leq M$,
the space $\mathcal{M}_{R,M}^{D}\coloneqq\\{\mu\in\mathcal{M}^{D}\text{ such
that }Supp(\mu)\in B(0,R),|\mu|\leq M\\}$ satisfies for all
$t\in[0,\inf\limits_{\mu\in\mathcal{M}_{R,M}^{D}}\varepsilon(\mu)],x\in\mathbb{R}^{m}$,
one has,

$W^{g}(S_{t}^{2}(x_{0},\mu),\gamma_{\mu}(t))\leq C(R,M)t^{2}.$ (32)

The following lemma plays an essential rule to prove the uniqueness of
solutions to system (12). The proof is entirely similar to that provided in
[15], thus we skip it.

######  Lemma 3.10.

Let
$S\colon[0,T]\times\mathbb{R}^{m}\times\mathcal{M}(\mathbb{R}^{n})\to\mathbb{R}^{m}\times\mathcal{M}(\mathbb{R}^{n})$
be a Lipschitz semigroup for system (12),
$\mu\colon[0,T]\to\mathcal{M}(\mathbb{R}^{n})$ a curve that is Lipschitz
continuous, and $x_{0}\in\mathbb{R}^{m}$. Then it holds

$W^{g}(S_{t}^{2}(x_{0},\mu(0)),\mu(t))\leq
e^{Ct}\int_{0}^{t}\liminf\limits_{h\to
0^{+}}W^{g}(S_{h}^{2}(x_{0},\mu(s)),\mu(s+h))\,ds$

where $C$ is the Lipschitz constant in equation (16).

We are now ready to prove the uniqueness of a semigroup to system (12)
conmpatible with a Dirac germ.

######  Theorem 3.11 (Uniqueness of Lipschitz semigroup of solutions to (12)).

Under the same assumptions of Theorem 3.4, assume a Dirac germ $\gamma$ as in
Defintion 3.8 is given. Then there exists at most one Lipschitz semigroup for
system (12) that is compatible with $\gamma$ as in Defintion 3.9.

###### Proof.

Fix initial datum
$(x_{0},\mu_{0})\in\mathbb{R}^{m}\times\mathcal{M}^{+}(\mathbb{R}^{n})$ and
$T>0$. Assume that there exist two semigroups $S,\bar{S}$ for system (12) that
are compatible with the given germ $\gamma$. By equation (15), there exists
$R>0$, such that $Supp(S_{t}^{2}(x_{0},\mu_{0}))\cup
Supp(\bar{S}_{t}^{2}(x_{0},\mu_{0})\subset B(0,R)$ for all $t\in[0,T]$. By
Lemma 3.10, we have,

$\displaystyle W^{g}(S_{t}^{2}(x_{0},\mu_{0}),\bar{S}_{t}^{2}(x_{0},\mu_{0}))$
$\displaystyle\leq$ $\displaystyle e^{Ct}\int_{0}^{t}\liminf\limits_{h\to
0^{+}}W^{g}(S_{h}^{2}(x_{0},\bar{S}_{s}^{2}(x_{0},\mu_{0})),\bar{S}_{s+h}^{2}(x_{0},\mu_{0}))\,ds.$
(33)

Fix $s>0$ and let $\nu=\bar{S}_{s}^{2}(x_{0},\mu_{0})$. Observe that
$\mathcal{M}^{D}$ is dense in $\mathcal{M}$ with respect to the topology
induced by $W^{g}$. Thus for every $\varepsilon>0$, there exists
$\bar{\nu}\in\mathcal{M}^{D}$ such that $W^{g}(\nu,\bar{\nu})<\varepsilon$.
Applying equation (32) to both $S_{h}^{2}$ and $\bar{S}_{h}^{2}$, we have,
there exists a constant $C_{1}(R,M)$ such that

$W^{g}(S_{h}^{2}(x_{0},\bar{\nu}),\gamma_{\bar{\nu}}(h))\leq
C_{1}(R,M)h^{2},\quad
W^{g}(\bar{S}_{h}^{2}(x_{0},\bar{\nu}),\gamma_{\bar{\nu}}(h))\leq
C_{1}(R,M)h^{2}.$

Therefore

$\displaystyle
W^{g}(S_{h}^{2}(x_{0},\bar{S}_{s}^{2}(x_{0},\mu_{0})),\bar{S}_{s+h}^{2}(x_{0},\mu_{0}))=W^{g}(S_{h}^{2}(x_{0},\nu),\bar{S}_{h}^{2}(x_{0},\nu))$
$\displaystyle\leq$ $\displaystyle
W^{g}(S_{h}^{2}(x_{0},\nu),S_{h}^{2}(x_{0},\bar{\nu}))+W^{g}(S_{h}^{2}(x_{0},\bar{\nu}),\gamma_{\bar{\nu}}(h))$
$\displaystyle+W^{g}(\gamma_{\bar{\nu}}(h),\bar{S}_{h}^{2}(x_{0},\bar{\nu}))+W^{g}(\bar{S}_{h}^{2}(x_{0},\bar{\nu}),\bar{S}_{h}^{2}(x_{0},\nu))$
$\displaystyle\leq$ $\displaystyle 2(e^{Ch}\varepsilon+C_{1}(R,M)h^{2}).$

Note that since both $s$ and $\varepsilon$ were chosen arbitrarily, for each
$s>0$ it holds

$\liminf\limits_{h\to
0^{+}}W^{g}(S_{h}^{2}(x_{0},\bar{S}_{s}^{2}(x_{0},\mu_{0})),\bar{S}_{s+h}^{2}(x_{0},\mu_{0}))=0.$
(34)

Combining equations (3) and (34), for every $t\in[0,T]$ it holds

$S_{t}^{2}(x_{0},\mu_{0})=\bar{S}_{t}^{2}(x_{0},\mu_{0}).$ (35)

Furthermore, for a.e. $t\in[0,T]$,

$S_{t}^{1}(x_{0},\mu_{0})=x_{0}+\int_{0}^{t}g(S_{\tau}^{1}(x_{0},\mu_{0}),S_{\tau}^{2}(x_{0},\mu_{0}))\,d\tau,$

and

$\bar{S}_{t}^{1}(x_{0},\mu_{0})=x_{0}+\int_{0}^{t}g(\bar{S}_{\tau}^{1}(x_{0},\mu_{0}),\bar{S}_{\tau}^{2}(x_{0},\mu_{0}))\,d\tau.$

Thus by assumption (OM:Lip-g) and equation (35), we conclude that for every
$t\in[0,T]$, $S_{t}^{1}(x_{0},\mu_{0})=\bar{S}_{t}^{1}(x_{0},\mu_{0})$. ∎

##  4\. Model for viral infections with mutations

In this section we introduce a model for viral infections, where mutations are
occurring as viral dynamics in the infected population.

We assume that the population of infected is represented by a measure over a
space of virus mutations. For simplicity we start assuming that the virus
mutations can be parameterized by one parameter $\alpha\in\mathbb{R}$ thus
$I\in\mathcal{M}(\mathbb{R}(\alpha))$. On the other side, the population of
susceptible can be identified, as usual by a single parameter, thus
$S\in\mathbb{R}$. The dynamics of infections is captured by and ODE:

$\dot{S}=-\frac{S}{N}\int_{\mathbb{R}}\beta(\alpha)\ dI(\alpha),$

where $N$ is the total population containing the population of susceptible
$S$, the population of infected $I$ and the population of recovered $R$,
$\beta(\alpha)$ is the infectivity rate, which depends on the virus mutation
identified by the parameter $\alpha$.

We now define a generalization of Example 2.10 to be applied to the measure
with time-varying mass $I(t)$. First set
$G_{I}(x)=\frac{I(]-\infty,x])}{I(\mathbb{R})}$, which is the normalized
cumulative distribution, so that $G_{I}(-\infty)=0$ and $G_{I}(+\infty)=1$.
Fix an increasing map $\varphi:[0,1]\to\mathbb{R}$ and define
$V_{\varphi}[I]=I\otimes_{x}J_{\varphi}(x)$, where

$J_{\varphi}(x)=\begin{cases}\delta_{\varphi(G_{I}(x))}&\mbox{~{}~{}if~{}}G_{I}(x^{-})=G_{I}(x),\\\
\frac{\varphi\\#\left(\chi_{[G_{I}(x^{-}),G_{I}(x)]}\lambda\right)}{G_{I}(x)-G_{I}(x^{-})}&\mbox{~{}~{}otherwise.}\end{cases}$
(36)

The dynamics for $I$ is given by the MDE:

$\dot{I}=V_{\varphi}[I]+\frac{S}{N}\beta(\alpha)I-\nu(\alpha)I$

where $\nu(\alpha)$ is the recovery rate, also dependent on the virus
mutation. Finally a second ODE described the dynamics of the population of
recovered, $R$, by:

$\dot{R}=\int_{\mathbb{R}}\nu(\alpha)\ dI(\alpha).$

The overall dynamics consists of coupled ODEs and MDE:

$\displaystyle\dot{S}=-\frac{S}{N}\int_{\mathbb{R}}\beta(\alpha)\ dI(\alpha),$
$\displaystyle\dot{I}=V_{\varphi}[I]+\frac{S}{N}\beta(\alpha)I-\nu(\alpha)I,$
$\displaystyle\dot{R}=\int_{\mathbb{R}}\nu(\alpha)\ dI(\alpha).$ (37)

######  Proposition 4.1.

Consider the system (4), with $V_{\varphi}$ given by (36), and assume
$\varphi$, $\beta$ and $\nu$ Lipschitz continuous and bounded. Then there
exists a Lipschitz semigroup of solutions. Moreover there exists unique
solutions compatible with Euler-LAS approximations.

###### Proof.

We first show that assumptions in the statement of Theorem 3.4 are satisfied.
The bounds on $\varphi$, $\beta$ and $\nu$ implies that all solutions to a
Cauchy problem is uniformly bounded in time. Therefore the assumptions of
Theorem 3.4 can be verified on bounded sets of
$\mathbb{R}^{2}\times\mathcal{M}(\mathbb{R})$ (the space of $(S,I,R)$) for the
metric given by the Euclidean norm on $\mathbb{R}^{2}$ and $W^{g}$ on
$\mathcal{M}(\mathbb{R})$. Indeed we show that all assumptions are valid on
the whole space except (OM:Lip-V) and (OM:Lip-s).

Since $\varphi$ is increasing, it is uniformly bounded by
$\max\\{|\varphi(0)|,|\varphi(1)|\\}$, thus (H:bound) is verified for
$V_{\varphi}$.

From the Lipschitz continuity and boundedness of $\beta$ and $\nu$ we get:

$\left|\int_{\mathbb{R}}\beta(\alpha)\
dI_{1}(\alpha)-\int_{\mathbb{R}}\beta(\alpha)\ dI_{2}(\alpha)\right|\leq
L(\beta)W^{g}(I_{1},I_{2})+M(\beta),$

where $L(\beta)$ is the Lipschitz constant of $\beta$ and $M(\beta)>0$ a
uniform bound. Moreover, the same estimate holds for $\nu$. Since the other
terms of the ODEs are Lipschitz, we conclude that (OM:Lip-g) is satisfied.

Now, consider $I_{1},I_{2}\in\mathcal{M}(\mathbb{R})$ and fix
$\tilde{I}_{1}\leq I_{1}$ and $\tilde{I}_{2}\leq I_{2}$ realizing
$W^{g}(I_{1},I_{2})$. Choose $\tilde{T}\in
P(V_{\varphi}(\tilde{I_{1}}),V_{\varphi}(\tilde{I_{2}}))$ such that
$\pi_{13}\\#\tilde{T}=\min\\{G_{\tilde{I}_{1}}(x),G_{\tilde{I}_{2}}(y)\\}\in
P^{opt}(\tilde{I_{1}},\tilde{I_{2}})$ (see Theorem 2.18 and Remark 2.19 of
[21]). Since $|\tilde{I}_{1}|=|\tilde{I}_{2}|$, by definition of $V_{\varphi}$
we get:

$\int|v-w|\ d\tilde{T}(x,v,y,w)=0.$ (38)

Define $\tilde{V}_{1}\leq V_{\varphi}(I_{1})$ such that
$\pi_{1}(\tilde{V}_{1})=\tilde{I}_{1}$ and $\int_{v}v\,d\tilde{V}_{1}(x,v)$ is
minimum among the $V\in\mathcal{M}(T\mathbb{R})$ such that
$\pi_{1}(V)=\tilde{I}_{1}$. Similarly define $\tilde{V}_{2}$. In simple words
we select the minimum speeds. Now choose $T\in P(\tilde{V}_{1},\tilde{V}_{2})$
such that $\pi_{13}\\#T=\min\\{G_{\tilde{I}_{1}}(x),G_{\tilde{I}_{2}}(y)\\}$.
Notice that $\pi_{13}\\#T\in P^{opt}(\tilde{I_{1}},\tilde{I_{2}})$.

By definition:

$\mathcal{W}^{g}(V_{\varphi}(I_{1}),V_{\varphi}(I_{2}))\leq\int|v-w|\,dT(x,v,y,w),$

indeed $\pi_{13}\\#T\in P^{opt}(\tilde{I}_{1},\tilde{I}_{2})$ and
$(\tilde{I}_{1},\tilde{I}_{2})$ realizes $W^{g}(I_{1},I_{2})$.
Now choose any $T_{1}\in P(\tilde{V}_{1},V_{\varphi}(\tilde{I_{1}}))$ such
that $\pi_{13}\\#T_{1}=Id$, the transference plan corresponding to the
identity map $Id(x)=x$. The existence of such $T_{1}$ comes from the fact that
$\pi_{1}(\tilde{V}_{1})=\pi_{1}(V_{\varphi}(\tilde{I_{1}}))=\tilde{I_{1}}$. We
choose similarly $T_{2}\in P(V_{\varphi}(\tilde{I_{2}}),\tilde{V}_{2})$ such
that $\pi_{13}\\#T_{2}=Id$. Notice that $T=T_{1}\circ\tilde{T}\circ T_{2}$,
thus we can write:

$\displaystyle\int|v-w|\,dT(x,v,y,w)\leq$
$\displaystyle\int|v-w|\,dT_{1}(x,v,y,w)+\int|v-w|\,d\tilde{T}(x,v,y,w)+\int|v-w|\,dT_{2}(x,v,y,w)$
$\displaystyle\doteq A_{1}+\tilde{A}+A_{2}.$

From (38) it holds $\tilde{A}=0$. Moreover we have:

$A_{1}\leq\int|v-w|\,d(V_{\varphi}(I_{1})\times
V_{\varphi}(\tilde{I}_{1})),\quad
A_{2}\leq+\int|v-w|d(V_{\varphi}(I_{2})\times V_{\varphi}(\tilde{I}_{2})).$
(39)

Now we notice that, since $\varphi$ is Lipschitz continuous, we have:

$|\varphi(G_{I_{1}}(x))-\varphi(G_{\tilde{I}_{1}}(x))|\leq
L(\varphi)|G_{I_{1}}(x)-G_{\tilde{I}_{1}}(x)|$

where $L(\varphi)$ is a Lipschitz constant for $\varphi$, thus:

$\displaystyle\int|v-w|\,d(V_{\varphi}(I_{1})\times
V_{\varphi}(\tilde{I}_{1}))\leq
L(\varphi)\int|x|\,d(G_{I_{1}}(x)-G_{\tilde{I}_{1}}(x))$ $\displaystyle\leq
L(\varphi)\left(sup_{x\in Supp(I_{1})}|x|\right)|I_{1}-\tilde{I}_{1}|.$ (40)

The same estimate holds for $I_{2}$ and $\tilde{I}_{2}$. Using (39) and (4),
we get that $A_{1}$ and $A_{2}$ are bounded by
$C(|I_{1}-\tilde{I}_{1}|+|I_{2}-\tilde{I}_{2}|)$ for some positive $C>0$
bounded over bounded sets. Moreover,
$W^{g}(I_{1},I_{2})\geq|I_{1}-\tilde{I}_{1}|+|I_{2}-\tilde{I}_{2}|$. Therefore
(OM:Lip-V) holds on bounded sets.

The boundedness of $\beta$ and $\nu$ imply (OM:bound-s).

Notice that the source term is given by $(\beta\frac{S}{N}-\nu)I$, thus:

$\displaystyle
W^{g}(s((S_{1},R_{1}),I_{1}),s((S_{2},R_{2}),I_{2}))\leq|I_{1}||(S_{1},R_{1})-(S_{2},R_{2})|$
$\displaystyle+|(S_{2},R_{2})|W^{g}(I_{1},I_{2}).$

Therefore (OM:Lip-s) is verified on bounded subsets of
$\mathbb{R}^{2}\times\mathcal{M}(\mathbb{R})$. This concludes the proof of the
first statement.

Now notice that all Euler-LAS approximations give rise to a Cauchy sequence in
$\mathbb{R}^{2}\times\mathcal{M}(\mathbb{R})$ for the metric given by the
Euclidean norm on $\mathbb{R}^{2}$ and $W^{g}$ on $\mathcal{M}(\mathbb{R})$.
Therefore, we conclude by completeness of the same metric. ∎

######  Lemma 4.2.

Consider the system (4), with $V_{\varphi}$ given by (36), and assume
$\varphi$, $\beta$ and $\nu$ Lipschitz continuous and bounded. Then the
quantity $S+|I|+R$ is conserved along solutions.

###### Proof.

It is sufficient to notice that the total mass of the positive source for $I$
in (4) coincides with the absolute value of the right-hand side for $S$.
Similarly the total mass of the negative source for $I$ in coincides with the
right-hand side for $R$. ∎

######  Example 4.3.

Consider the model (4) and assume that $\beta$, $\nu$ are constant, i.e.
independent of $\alpha$. Then we have:

$\dot{S}=-\frac{S}{N}\,\beta\,|I|,\quad\dot{R}=\nu\,|I|.$

Since $V_{\varphi}$ does not change the total mass of $I$, defining
$m(t)=|I(t)|$, we get:

$\dot{m}(t)=\frac{S}{N}\,\beta\,m-\nu\,m,$

thus the triplet $(S,m,R)$ satisfy the classical SIR model.

######  Example 4.4.

Consider the model (4) with $V$ as in Example 2.8 and $I(0)=|I(0)|\delta_{0}$,
thus the mass of $I$ is split in two parts traveling with velocity $-1$ and
$+1$ respectively. Assume $\beta(\alpha)=\beta(-\alpha)$ and
$\nu(\alpha)=\nu(-\alpha)$, then we get
$\int_{\mathbb{R}}\beta(\alpha)dI(t)(\alpha)=\beta(t)$ and
$\int_{\mathbb{R}}\nu(\alpha)dI(t)(\alpha)=\nu(t)$. Therefore, defining
$m(t)=|I(t)|$, the solution to (4) satisfies the time-dependent SIR-model:

$\dot{S}=-\frac{S}{N}\,\beta(t)\,m(t),\quad\dot{m}=\frac{S}{N}\,\beta(t)\,m-\nu(t)\,m,\quad\dot{R}=\nu(t)\,m.$

## Acknowledgements

The authors would like to acknowledge the support of the NSF CMMI project #
2033580 “Managing pandemic by managing mobility” in collaboration with Cornell
University and Vanderbilt University, and the support of the Joseph and
Loretta Lopez Chair endowment.

## References

  * [1] S. Anita and V. Capasso,  Reaction-diffusion systems in epidemiology,  2017\.
  * [2] (MR4144366) [10.1142/S0218202520500323]  N. Bellomo, R. Bingham, M. A. J. Chaplain, G. Dosi, G. Forni, D. A. Knopoff, J. Lowengrub, R. Twarock and M. E. Virgillito,  A multi-scale model of virus pandemic: Heterogeneous interactive entities in a globally connected world,  _Math. Models Methods Appl. Sci._ , 30 (2020), 1591–1651, 2006.03915.
  * [3] (MR3348406) [10.1051/cocv/2014032]  G. A. Bonaschi, J. A. Carrillo, M. Di Francesco and M. A. Peletier,  Equivalence of gradient flows and entropy solutions for singular nonlocal interaction equations in 1D,  ESAIM Control Optim. Calc. Var., 21 (2015), 414–441.
  * [4] (MR4269269) [10.1126/science.abc6810]  T. Britton, F. Ball and P. Trapman,  A mathematical model reveals the influence of population heterogeneity on herd immunity to SARS-CoV-2,  Science, 369 (2020), 846–849.
  * [5] (MR4205631) [10.1109/TNSE.2020.3024723]  Y.-C. Chen, P.-E. Lu and C.-S. Chang,  A time-dependent SIR model for COVID-19,  _IEEE Trans. Network Sci. Eng._ , 7 (2020), 3279–3294, 2003.00122.
  * [6] (MR4283943) [10.1007/s00245-020-09660-9]  R. M. Colombo and M. Garavello,  Well posedness and control in a nonlocal SIR model,  _Appl. Math. Optim._ , 84 (2021), 737–771.
  * [7] (MR4139038) [10.1186/s13362-020-00090-4]  R. M. Colombo, M. Garavello, F. Marcellini and E. Rossi,  An age and space structured SIR model describing the COVID-19 pandemic,  Journal of Mathematics in Industry, 10 (2020), Paper No. 22, 20 pp.
  * [8] [10.1038/s41591-020-0883-7]  G. Giordano, F. Blanchini, R. Bruno, P. Colaneri, A. Di Filippo, A. Di Matteo and M. Colaneri,  Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy,  Nature Medicine, 26 (2020), 855–860.
  * [9] V. Kala, K. Guo, E. Swantek, A. Tong, M.. Chyba, Y. Mileyko, C. Gray, T. Lee and A. E. Koniges,  _Pandemics in Hawai’i: 1918 Influenza and COVID-19_ ,  The Ninth International Conference on Global Health Challenges GLOBAL HEALTH 2020, IARIA, 2020.
  * [10] A. Keimer and L. Pflug,  Modeling infectious diseases using integro-differential equations: Optimal control strategies for policy decisions and applications in covid-19,  2020\.
  * [11] [10.1098/rspa.1932.0171]  W. O. Kermack and A. G. McKendrick,  Contributions to the mathematical theory of epidemics. II.-the problem of endemicity,  Proceedings of the Royal Society of London. Series A, containing papers of a mathematical and physical character, 138 (1932), 55–83.
  * [12] [10.1126/science.373.6557.844]  K. Kupferschmidt,  Evolving threat,  Science, 373 (2021), 844–849.
  * [13] (MR4269258) [10.1126/science.abd1668]  C. J. E. Metcalf, D. H. Morris and S. W. Park,  Mathematical models to guide pandemic response,  Science, 369 (2020), 368–369.
  * [14] [10.1093/infdis/jiw375]  K. R. Moran, G. Fairchild, N. Generous, K. Hickmann, D. Osthus, R. Priedhorsky, J. Hyman and S. Y. Del Valle,  Epidemic forecasting is messier than weather forecasting: The role of human behavior and internet data streams in epidemic forecast,  The Journal of Infectious Diseases, 214 (2016), S404–S408.
  * [15] (MR3961299) [10.1007/s00205-019-01379-4]  B. Piccoli,  Measure differential equations,  Archive for Rational Mechanics and Analysis, 233 (2019), 1289–1317.
  * [16] (MR3182483) [10.1007/s00205-013-0669-x]  B. Piccoli and F. Rossi,  Generalized Wasserstein distance and its application to transport equations with source,  Archive for Rational Mechanics and Analysis, 211 (2014), 335–358.
  * [17] (MR3544329) [10.1007/s00205-016-1026-7]  B. Piccoli and F. Rossi,  On properties of the generalized Wasserstein distance,  Archive for Rational Mechanics and Analysis, 222 (2016), 1339–1365.
  * [18] (MR4026977) [10.3934/dcds.2019270]  B. Piccoli and F. Rossi,  Measure dynamics with probability vector fields and sources,  Discrete Contin. Dyn. Syst., 39 (2019), 6207–6230.
  * [19] [10.1126/science.abc5096]  N. W. Ruktanonchai, J. R. Floyd, S. Lai, C. W. Ruktanonchai, A. Sadilek, P. Rente-Lourenco, X. Ben, A. Carioli, J. Gwinn, J. E. Steele, et al.,  Assessing the impact of coordinated COVID-19 exit strategies across europe,  Science, 369 (2020), 1465–1470.
  * [20] A. Vespignani, H. Tian, C. Dye, J. O. Lloyd-Smith, R. M. Eggo, M. Shrestha, S. V. Scarpino, B. Gutierrez, M. U. G. Kraemer, J. Wu, et al.,  Modelling covid-19,  Nature Reviews Physics, 2 (2020), 279–281.
  * [21] (MR1964483) [10.1090/gsm/058]  C. Villani,  _Topics in Optimal Transportation,_ Graduate Studies in Mathematics, 58. American Mathematical Society, Providence, RI, 2003.
  * [22] [10.1126/science.abb8001]  J. Zhang, M. Litvinova, Y. Liang, Y. Wang, W. Wang, S. Zhao, Q. Wu, S. Merler, C. Viboud, A. Vespignani, et al.,  Changes in contact patterns shape the dynamics of the COVID-19 outbreak in China,  Science, 368 (2020), 1481–1486.
  * [23] [10.1016/S1473-3099(20)30230-9]  J. Zhang, M. Litvinova, W. Wang, Y. Wang, X. Deng, X. Chen, M. Li, W. Zheng, L. Yi, X. Chen, et al.,  Evolving epidemiology and transmission dynamics of coronavirus disease 2019 outside Hubei province, China: A descriptive and modelling study,  The Lancet Infectious Diseases, 20 (2020), 793–802.

Received May 2021; revised September 2021; early access March 2022.