# SIS epidemics on open networks: A replacement-based approximation

Renato Vizuete, Paolo Frasca, and Elena Panteley This work was supported by
F.R.S.-FNRS via the _KORNET_ project and via the Incentive Grant for
Scientific Research (MIS) _Learning from Pairwise Comparisons_ , and by the
_RevealFlight_ Concerted Research Action (ARC) of the Fédération Wallonie-
Bruxelles and in part by the Agence Nationale de la Recherche (ANR) via grants
HANDY ANR-18-CE40-0010 and COCOON ANR-22-CE48-0011.R. Vizuete is with ICTEAM
institute, UCLouvain, B-1348, Louvain-la-Neuve, Belgium. R. Vizuete is a FNRS
Postdoctoral Researcher - CR. P. Frasca is with Univ. Grenoble Alpes, CNRS,
Inria, Grenoble INP, GIPSA-lab, F-38000 Grenoble, France. E. Panteley is with
Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire des signaux et
systèmes, 91190, Gif-sur-Yvette, France<EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

In this paper we analyze continuous-time SIS epidemics subject to arrivals and
departures of agents, by using an approximated process based on replacements.
In defining the SIS dynamics in an open network, we consider a stochastic
setting in which arrivals and departures take place according to Poisson
processes with similar rates, and the new value of the infection probability
of an arriving agent is drawn from a continuous distribution. Since the system
size changes with time, we define an approximated process, in which
replacements take place instead of arrivals and departures, and we focus on
the evolution of an aggregate measure of the level of infection. So long as
the reproduction number is less than one, the long-term behavior of this
function measures the impact of the changes of the set of agents in the
epidemic. We derive upper bounds for the expectation and variance of this
function and we include a numerical example to show that the approximated
process is close to the original SIS process.

##  I Introduction

Epidemic models are dynamical systems used to model contagion in a community
of individuals. These types of models have long time been a key class of
network dynamics and their interest has further grown in the last few years,
specially due to the COVID-19 outbreak that impacted the entire world. One of
the most important epidemic models is the SIS (Susceptible-Infected-
Susceptible) compartmental epidemic model, which describes diseases where the
agents, after being recovered, have no immunity against re-infection [1, 2].

Most of the stability results derived using network models depend on the
spectral radius of the adjacency matrix associated with the network, assuming
that the composition of the system is fixed [3]. Nevertheless, in more
realistic environments, individuals do not stay in a fixed place for a long
time and are constantly moving between different locations. It is clear that
mobility is one of the most important factors that need to be considered in
the analysis of modern spreads of diseases [4], and many works begin to
consider this aspect to obtain more accurate results [5, 6]. At the level of
network models, time-varying networks have traditionally been used to model to
some degree the mobility of agents by changing the interactions among them [7,
8, 9, 10]. However, this type of models might fail to represent modern
networks where due to the development of transportation infrastructures,
agents are constantly entering and leaving towns, cities and countries, which
implies not only changes of the connections but also modifications of the
infection probabilities associated with the agents. Other approaches have
considered patchy environments where individuals move between layers and
participate in the local dynamics [11]. Nevertheless, a complete knowledge of
all the agents in the system and the transitions between populations might
demand a lot of information, which limits the applicability of this approach.

The analysis of systems characterized by a time-varying set of agents is the
object of study of _open multi-agent systems_ , where the dimension of the
system may vary with time [12, 13, 14]. At the level of large populations of
individuals (e.g., countries, cities), it is common to have similar rates of
arrivals and departures such that the variations of the size of the
populations are almost negligible. However, even if the total number of
individuals in a specific place could remain almost invariant, it is
unrealistic to assume that exactly the same agents stay in a fixed location
during all the evolution of the epidemic. If we consider a similar rate of
departures and arrivals, which preserve the number of individuals, epidemics
under these conditions can be approximated as a dynamical system subject to
replacements of individuals: a similar approach has been taken for opinion
dynamics by [15, 16].

In this paper, we study a continuous time formulation of an SIS epidemic over
a network of individuals when arrivals and departures of agents take place. We
consider a stochastic setting where the time instants at which the events
occur are determined by a homogeneous Poisson process and the infection
probabilities of arrival agents are drawn from a continuous distribution.
Because the dimension of the system changes with time, the analysis of the
epidemic is performed by using an aggregate function (i.e., scalar value) that
measures the level of infection. Due to the high complexity of the stochastic
process, when the rates of arrivals and departures are similar, we propose an
approximated version of the original process based on replacements. Upper
bounds for the asymptotic values of the expectation and variance of the
aggregate function are derived as functions of the parameters of the epidemic
process.

##  II SIS model on closed networks

We begin by recalling, in this section, the definition and main features of
the classical SIS model on “closed” networks (i.e., networks whose node set
does not change).

We consider a system composed by $n$ agents interacting through a connected
undirected graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ where the set of
agents is given by $\mathcal{V}=\\{1,\ldots,n\\}$ and the set of edges by
$\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$. The SIS epidemic model
over a network is given by [3]:

$\dot{x}_{i}(t)=-\delta x_{i}(t)+\sum_{j=1}^{n}a_{ij}\beta
x_{j}(t)\left(1-x_{i}(t)\right),$

where $x_{i}(t)\in[0,1]$ is the probability of infection of an individual at
time $t$, $\delta$ is the recovery rate, $\beta$ is the infection rate,
$a_{ji}=a_{ij}>0$ if there is an interaction between agents $i$ and $j$, and
$a_{ij}=0$ otherwise.

The model of all the network can be expressed as:

$\dot{x}(t)=\left(\beta A-\delta I\right)x(t)-\beta X(t)Ax(t),$ (1)

where $I$ is the identity matrix, $A=[a_{ij}]$ is the adjacency matrix of the
network, and $X(t)=\mathrm{diag}[x_{1}(t),\ldots,x_{n}(t)]$. Since the graph
is undirected, $A$ is symmetric and its eigenvalues are denoted, in non-
increasing order of magnitude as $\lambda_{1}\geq\cdots\geq\lambda_{n}$.

######  Proposition 1 (Stability [17])

For the SIS epidemic (1), the disease-free equilibrium $x=0$ is globally
asymptotically stable if and only if

$\lambda_{1}\frac{\beta}{\delta}<1.$ (2)

Since in open multi-agent systems the size of the system may change in time
due to decoupled arrivals and departures [12, 13, 14], it is desirable to
perform the analysis with the use of an aggregate function $V(x)$ that
condensates the associated process in a scalar value, independently of the
size of the system. In the case of SIS epidemics, a natural choice is the
Lyapunov function $V:[0,1]^{n}\to\mathbb{R}_{\geq 0}$:

$V(x)=\frac{1}{n}\left|\left|x\right|\right|^{2},$ (3)

where $\left|\left|\cdot\right|\right|$ is the Euclidean norm. This Lyapunov
function satisfies

$\lim_{t\to\infty}V(x(t))=0,$ (4)

when the epidemic is stable. Thanks to normalizing by the number of agents,
$V(x)$ enjoys uniform bounds, independent of the size of the system:

$0\leq V(x)\leq 1.$ (5)

In a closed system, where the evolution of the states of the agents is fully
determined by (1), the dynamics of the aggregate function $V(x(t))$ can be
upper bounded as follows (see [18]):

$\displaystyle\dot{V}(x(t))$ $\displaystyle=\frac{2}{n}x^{T}(t)\dot{x}(t)$
$\displaystyle=\frac{2}{n}x^{T}(t)\left(\beta A-\delta I\right)x(t)-\beta
x^{T}(t)X(t)Ax(t)$ (6) $\displaystyle\leq\frac{2}{n}x^{T}(t)\left(\beta
A-\delta I\right)x(t)$ $\displaystyle\leq 2(\beta\lambda_{1}-\delta)V(x(t)),$
(7)

where we used the fact that $x^{T}(t)X(t)Ax(t)\geq 0$ since $x_{i}\in[0,1]$.
Clearly, the right-hand side of (7) is negative if
$\beta\lambda_{1}-\delta<0$, which corresponds to the stability condition (2).
Finally, by applying the Comparison Lemma, we have that $V(x(t))$ satisfies:

$V(x(t))\leq V_{0}+2(\beta\lambda_{1}-\delta)\int_{0}^{t}V(x(\tau))d\tau,$ (8)

where $V_{0}$ denotes the value of $V(x(t))$ at the time $t=0$. Similar
bounds, in differential and integral form, will be derived for open SIS
epidemics in the rest of this paper.

##  III Open SIS epidemics

In this section, we formulate a SIS epidemic where in addition to the
continuous evolution of the states of the individuals according to (1), the
system also experiences arrivals and departures of individuals. Then, in the
next section we will propose an alternative process, characterized by
replacements in a fixed network, as an approximation in the case of similar
rates of arrivals and departures.

######  Definition 1 (Departure)

We say that an agent $j\in\mathcal{V}(t^{-})$ leaves the epidemic at time $t$
if:

$\mathcal{V}(t^{+})=\mathcal{V}(t^{-})\setminus\\{j\\},$

where $\mathcal{V}(t^{-})$ is the set of agents before the departure and
$\mathcal{V}(t^{+})$ is the set of agents after the departure of agent $j$.
Thus, $\left|\mathcal{V}(t^{+})\right|=\left|\mathcal{V}(t^{-})\right|-1$.

######  Assumption 1 (Departure process)

The departure  instants in the system are determined by a Poisson process
$N_{t}^{(d)}$ with rate $\mu_{d}>0$ for $n>1$ and $\mu_{d}=0$ for $n=1$.
During a departure, an agent $j$ in the current epidemic is chosen uniformly
and removed from the system.

Notice that if there is only one agent in the epidemic, we keep this agent and
we set to 0 the rate of the Poisson process for departures such that only
arrivals can happen. In this way, the departure process becomes a slightly
truncated version of a homogeneous Poisson process.

######  Definition 2 (Arrival)

We say that an agent $j$ joins the epidemic at time $t$ if:

$\mathcal{V}(t^{+})=\mathcal{V}(t^{-})\cup\\{j\\},$

where $\mathcal{V}(t^{-})$ is the set of agents before the arrival and
$\mathcal{V}(t^{+})$ is the set of agents after the arrival of agent $j$.
Thus, $\left|\mathcal{V}(t^{+})\right|=\left|\mathcal{V}(t^{-})\right|+1$.

######  Assumption 2 (Arrival process)

The arrival instants are determined by a homogeneous Poisson process
$N_{t}^{(a)}$ with rate $\mu_{a}>0$. During an arrival, an agent $j$ joins the
epidemic with an infection probability $x_{j}$ determined by a random variable
$\Theta$, which takes values according to a continuous distribution with
support in the interval $[0,1]$, with mean $m$ and variance $\sigma^{2}$. The
new agent generates edges with each current node with probability $p$.

Notice that at any time the graphs are realizations of an Erdős-Rényi graph
and therefore $\lambda_{1}$ grows with the number of agents. We consider that
the recovery rate remains constant $\delta_{n}=\bar{\delta}$ and the infection
rate is given by $\beta_{n}=n^{-1}\bar{\beta}$, which means that, as the graph
grows in size, the healing rate (which depends on each individual) remains
constant, whereas the infection rate decreases. This natural scaling law is
also chosen in [19, 20]. Indeed, on dense graphs this assumption means that in
larger graphs, even though there are more potential interactions, the average
strength of the connections is suitably adjusted: this fact accounts for
natural limitations in the rates of contact between individuals.

Similarly to [19], the objective of the aggregate function $V$ is to quantify
the reaction of the epidemic to arrivals and departures that affect the
evolution of the disease, such that any deviation from the equilibrium $x=0$
is due to the openness of the system. The evolution of $V$ in this open
setting is given by a stochastic differential equation (SDE) of the form:

$dV=f(V(t))dt+g_{a}(V(t))dN_{t}^{(a)}+g_{d}(V(t))dN_{t}^{(d)},$ (9)

where $f$ determines the dynamics of $V$ in continuous time, $g_{a}$
determines the change of $V$ during an arrival and $g_{d}$ determines the
change of $V$ during a departure. This SDE is also known as a _Poisson Counter
driven Stochastic Differential Equation_ [21, 22], and its solution is a
cadlag function associated with a Piecewise-deterministic Markov Process
(PDMP) such that it can also be studied using the theory of PDMP [23].

Each realization of the stochastic process (9) can be completely different,
and the analysis of the asymptotic behavior is not appropriate because the
limit $\lim_{t\to\infty}V(x(t))$ does not exist due to the continuous arrivals
and departures. For this reason, we are interested in the statistical
properties of (9), corresponding to $\mathbb{E}\left[V(x(t))\right]$. However,
notice that the stochastic process (9) is highly complex since the mapping $f$
determined by (6) is time-varying due to $A(t)$, and depends on the two
Poisson process $N_{t}^{(a)}$ and $N_{t}^{(d)}$. In addition, due to the
slightly truncation at $n=1$, the Poisson process $N_{t}^{(d)}$ is not
homogeneous.

For this reason, we look for a more tractable mathematical process that
approximates the behavior of the original process (9).

##  IV Replacements as an approximation

Due to the complexity of the stochastic process, we focus on the particular
case of similar rates of arrivals and departures and we perform an analysis of
the evolution of $V$ considering the following two approximations:

#### Replacements

In large populations, the rates of arrivals and departures are similar such
that even if the individuals are constantly changing, the total number of
agents remains the same in expectation. For this reason, we approximate
arrivals and departures as replacements.

#### Fixed network

The connections of all the agents are given by an Erdős-Rényi graph. This
means that all the graph topologies generated in the stochastic process during
arrivals have the same expectation
$\bar{A}=p(\mathds{1}_{n}\mathds{1}^{T}_{n}-I_{n})$, where $\mathds{1}_{n}$ is
a vector of ones of size $n$. In addition, in large networks, the eigenvalues
of the expected graph are close to those of the random graphs [24]. For this
reason, we approximate the time-varying matrices $A(t)$ with a fixed matrix
$\bar{A}$ of size $n_{0}$, where $n_{0}$ denotes the initial number of agents
$n(t)$ at $t=0$.

These approximations are made formal by the following definitions.

######  Definition 3 (Replacement)

We say that an agent $j\in\left\\{1,\ldots,n_{0}\right\\}$ is replaced at time
$t$ if:

$x^{+}=[x_{1}^{-},\cdots,x_{j-1}^{-},{\color[rgb]{0,0,0}b_{j}},x_{j+1}^{-},\cdots,x_{n}^{-}]^{T},\qquad{\color[rgb]{0,0,0}x_{j}^{-}\neq
b_{j}},$

where $x^{-}=x(t^{-})$ is the state of the system before the replacement,
$x^{+}=x(t^{+})$ is the state of the system after the replacement and $b_{j}$
is the probability of infection of the new agent.

###### Assumption 3

The set of edges $\mathcal{E}$ of the network remains invariant for all time
$t$.

In the following proposition we analyze the variation of the aggregate
function $V(x(t))$ under a replacement event where the replaced agent is
chosen uniformly (i.e., the probability of replacement of each agent is the
same) and the value of the new agent after the replacement is determined by a
random variable $\Theta$ following a continuous distribution with support in
the interval $[0,1]$, with mean $m$ and variance $\sigma^{2}$.

######  Proposition 2 (Replacement)

During the replacement of an agent, the aggregate function $V(x(t))$ defined
in (3) satisfies:

$\mathbb{E}\left[V(x(t^{+}))\right]\\!-\\!\mathbb{E}\left[V(x(t^{-}))\right]\\!=\\!\frac{1}{n_{0}}\left(\sigma^{2}+m^{2}-\mathbb{E}\left[V(x(t^{-}))\right]\right).$
(10)

###### Proof:

Due to the space constraints, we avoid including the time dependence of $x(t)$
and the state dependence of $V(x(t))$. This abuse of notation will be made
multiple times in the rest of the paper. Let us assume that an agent $i$ is
replaced. We begin by computing the conditional expectation of
$V^{+}=V(x(t^{+}))$ given $x^{-}$ and the value of the replaced agent
$\Theta$. Then, it holds

$\displaystyle\mathbb{E}\left[V^{+}|x^{-},\Theta\right]$
$\displaystyle=\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\left(\frac{1}{n_{0}}\sum_{j\neq
i}(x_{j}^{2})^{-}+\frac{1}{n_{0}}\Theta^{2}\right)$
$\displaystyle=\frac{1}{n_{0}^{2}}\sum_{i=1}^{n_{0}}\left(\left|\left|x^{-}\right|\right|^{2}-(x_{i}^{2})^{-}+\Theta^{2}\right)$
$\displaystyle=\frac{1}{n_{0}}\left|\left|x^{-}\right|\right|^{2}-\frac{1}{n_{0}^{2}}\left|\left|x^{-}\right|\right|^{2}+\frac{\Theta^{2}}{n_{0}}$
$\displaystyle=V^{-}-\frac{1}{n_{0}}V^{-}+\frac{\Theta^{2}}{n_{0}},$

where $V^{-}$ denotes $V(x(t^{-}))$. Then, by computing the total expectation,
we complete the proof. ∎

######  Remark 1 (Other systems)

The proof of Proposition 2 did not use equation (1), only the definition of
$V$. Therefore, the statement remains valid for other systems characterized by
replacements of agents.

Along this work, we make the following assumption about the occurrence of
replacements in the evolution of the epidemic.

######  Assumption 4 (Replacement process)

The replacements instants are determined by a homogeneous Poisson process with
rate $\mu>0$. During a replacement, an agent
$j\in\left\\{1,\ldots,n_{0}\right\\}$ is chosen uniformly and it is assigned a
new infection probability $b_{j}$ determined by a random variable $\Theta$,
which takes values according to a continuous distribution with support in the
interval $[0,1]$, with mean $m$ and variance $\sigma^{2}$.

###  IV-A Pure replacements ($\dot{x}(t)=0$)

Since the replacements process is independent of the evolution of the epidemic
given by (1), we analyze the behavior of the aggregate function $V$ during
replacements independently of the continuous time dynamics by considering
$\dot{x}(t)=0$. Notice that the jumps of the state $x(t)$ corresponds to the
jump of the aggregate function $V(x(t))$. In this case, the value of $V(x(t))$
is given by:

$\displaystyle V(x(t))$ $\displaystyle=V_{0}+\sum_{k=1}^{N_{t}}\Delta
V_{k}=V_{0}+\sum_{k=1}^{N_{t}}\left(V(x(T_{k}^{+}))-V(x(T_{k}^{-}))\right),$

where $\Delta V_{k}$ is the size of the jump of the function $V(x(t))$ at the
$T_{k}$ jump time of the Poisson process $N_{t}$ and

$V(x(T_{k}^{-}))=\lim_{t\uparrow
T_{k}}V(x(t))\qquad\\!\\!\\!\text{and}\qquad\\!\\!\\!V(x(T_{k}^{+}))=\lim_{t\downarrow
T_{k}}V(x(t)).$

The aggregate function satisfies $V(x(t^{-}))=V(x(t^{+}))$ for almost all $t$,
except in a countable number of time jumps. Therefore, it can be rewritten as
[25]:

$\displaystyle V(x(t))$
$\displaystyle=V_{0}+\sum_{k=1}^{N_{t}}\left(V(x(T_{k}^{+}))-V(x(T_{k}^{-}))\right)$
$\displaystyle=V_{0}+\int_{0}^{t}(V(x(\tau^{+}))-V(x(\tau^{-})))dN_{\tau}.$

A homogeneous Poisson process satisfies $\mathbb{E}\left[N_{t}\right]=\mu t$,
which gives us
$\mathbb{E}\left[dN_{t}\right]=d\mathbb{E}\left[N_{t}\right]=\mu dt$ [25].
Then, it holds

$\mathbb{E}\left[V(x(t))\right]=V_{0}+\int_{0}^{t}\mathbb{E}\left[V(x(\tau^{+}))-V(x(\tau^{-}))\right]\mu
d\tau,$

which is the integral version of the ODE:

$\frac{d}{dt}\mathbb{E}\left[V(x(t))\right]=\mu\mathbb{E}\left[V(x(t^{+}))-V(x(t^{-}))\right].$
(11)

######  Proposition 3 (Pure replacement process)

Under Assumption 4 and inter-event dynamics $\dot{x}=0$, the aggregate
function $V$ defined in (3) satisfies:

$\mathbb{E}\left[V(x(t))\right]=\left(V_{0}-\sigma^{2}-m^{2}\right)e^{-\frac{\mu}{n_{0}}t}+(\sigma^{2}+m^{2}).$
(12)

###### Proof:

The result is obtained by applying (10) in (11)

$\frac{d}{dt}\mathbb{E}\left[V(x(t))\right]=-\frac{\mu}{n_{0}}\mathbb{E}\left[V(x(t))\right]+\frac{\mu}{n_{0}}(\sigma^{2}+m^{2}),$

and solving the ordinary differential equation (ODE). ∎

Proposition 3 shows that in a system subject only to replacements, the
expected value of the aggregate function is bounded and its asymptotic value
is:

$\lim_{t\to\infty}\mathbb{E}\left[V(x(t))\right]=\sigma^{2}+m^{2},$ (13)

which corresponds to the second moment of the distribution used to generate
the new values of the replaced agents. The rate of the Poisson process $\mu$
only has an influence on the rate of convergence of $V(x(t))$ such that a
large value of $\mu$ will guarantee a fast convergence. The asymptotic value
of $V(x(t))$ is independent of the rate $\mu$.

##  V First moment

Now, we consider the behavior of the aggregate function in continuous time
subject to replacements at time instants determined by a Poisson process.
Since the solution to (6) is unique and the jumps are determined by a
homogeneous Poisson process, the aggregate function $V(x(t))$ is given by
[21]:

$V(x(t))=V_{0}+\int_{0}^{t}\frac{2}{n_{0}}\big{(}x^{T}(\tau)(\beta\bar{A}\\\
-\delta I)x(\tau)-\beta
x^{T}(\tau)X(\tau)\bar{A}x(\tau)\big{)}d\tau+\int_{0}^{t}\Delta
V_{k}dN_{\tau},$ (14)

which can be expressed as the SDE

$\\!\\!\\!\\!\\!\\!dV(x(t))=\frac{2}{n_{0}}\left(x^{T}(t)\left(\beta\bar{A}-\delta
I\right)x(t)-\beta x^{T}(t)X(t)\bar{A}x(t)\right)dt\\\ +\Delta V_{k}dN_{t}.$
(15)

Notice that by denoting

$\omega(x(t))=\frac{2}{n_{0}}\left(x^{T}(t)\left(\beta\bar{A}-\delta
I\right)x(t)-\beta x^{T}(t)X(t)\bar{A}x(t)\right),$

(14) can be expressed as

$V(x(t))\\!=\\!V_{0}+\\!\\!{\int_{0}^{T_{1}}\\!\\!\omega(x(\tau))d\tau+\\!\\!\int_{T_{1}}^{T_{2}}\\!\\!\omega(x(\tau))d\tau+\cdots}+\\!\sum_{k=1}^{N_{t}}\\!\Delta
V_{k},$

which shows that $V(x(t))$ jumps at the time instants $T_{k}$ determined by
the Poisson process, corresponding to the sum, and between these time
intervals, $V(x(t))$ evolves according to (6), corresponding to the integrals.

######  Theorem 1 (Expectation)

Consider a SIS epidemic satisfying the stability condition (2). Under
Assumptions 3 and 4, the aggregate function $V(x(t))$ defined in (3)
satisfies:

$\limsup_{t\to\infty}\mathbb{E}\left[V(x(t))\right]\leq\frac{\mu(\sigma^{2}+m^{2})}{\mu+2n_{0}\bar{\delta}-2\bar{\beta}p(n_{0}-1)}.$
(16)

###### Proof:

We take the expectation of both sides of (15) and we obtain:

$d\mathbb{E}\left[V\right]=\mathbb{E}\left[\frac{2}{n_{0}}\left(x^{T}\left(\beta\bar{A}-\delta
I\right)x-\beta x^{T}X\bar{A}x\right)\right]dt\\\
+\frac{\mu}{n_{0}}\left(\sigma^{2}+m^{2}-\mathbb{E}\left[V\right]\right)dt,$

where we applied the Dominated Convergence Theorem to interchange the
expectation and the derivative
$\mathbb{E}\left[dV\right]=d\mathbb{E}\left[V\right]$ and we used the fact
that $\mathbb{E}\left[\Delta
V\right]=\frac{1}{n_{0}}\left(\sigma^{2}+m^{2}-\mathbb{E}\left[V\right]\right)$
according to Proposition 2. Then we obtain the following ODE:

$\displaystyle\dfrac{d}{dt}\mathbb{E}\left[V\right]$
$\displaystyle=\mathbb{E}\left[\frac{2}{n_{0}}\left(x^{T}\left(\beta\bar{A}-\delta
I\right)x-\beta x^{T}X\bar{A}x\right)\right]$
$\displaystyle\;\;\;\;+\frac{\mu}{n_{0}}\left(\sigma^{2}+m^{2}-\mathbb{E}\left[V\right]\right)$
$\displaystyle\leq
2(\beta\bar{\lambda}_{1}-\delta)\mathbb{E}\left[V\right]+\frac{\mu}{n_{0}}\left(\sigma^{2}+m^{2}-\mathbb{E}\left[V\right]\right)$
(17)
$\displaystyle=\left(\frac{2n_{0}\beta\bar{\lambda}_{1}-2n_{0}\delta-\mu}{n_{0}}\right)\mathbb{E}\left[V\right]+\frac{\mu(\sigma^{2}+m^{2})}{n_{0}},$
(18)

where $\bar{\lambda}_{1}$ is the largest eigenvalue of $\bar{A}$. By the
Comparison Lemma, we can guarantee that $\mathbb{E}\left[V\right]$ is upper
bounded by the solution of the right-hand side of (18). Then, it holds

$\displaystyle\mathbb{E}\left[V\right]$
$\displaystyle\leq\\!e^{\left(\frac{2n_{0}\beta\bar{\lambda}_{1}-2n_{0}\delta-\mu}{n_{0}}\right)t}\\!\\!\int_{0}^{t}\frac{\mu(\sigma^{2}\\!+\\!m^{2})}{n_{0}}e^{-\left(\frac{2n_{0}\beta\bar{\lambda}_{1}-2n_{0}\delta-\mu}{n_{0}}\right)\tau}d\tau$
$\displaystyle\quad+\mathbb{E}\left[V_{0}\right]e^{\left(\frac{2n_{0}\beta\bar{\lambda}_{1}-2n_{0}\delta-\mu}{n_{0}}\right)t}$
$\displaystyle=\left(\mathbb{E}\left[V_{0}\right]-\frac{\mu(\sigma^{2}+m^{2})}{\mu+2n_{0}\delta-2n_{0}\beta\bar{\lambda}_{1}}\right)\\!e^{\left(\frac{2n_{0}\beta\bar{\lambda}_{1}-2n_{0}\delta-\mu}{n_{0}}\right)t}$
$\displaystyle\quad+\frac{\mu(\sigma^{2}+m^{2})}{\mu+2n_{0}\delta-2n_{0}\beta\bar{\lambda}_{1}}.$
(19)

Finally, the result follows by taking the limit $t\to\infty$ in (19) and using
$\bar{\lambda}_{1}=p(n_{0}-1)$, $\delta=\bar{\delta}$ and
$\beta=n_{0}^{-1}\bar{\beta}$. ∎

Clearly, the bound (16) is a function of the parameters of the SIS epidemic
(i.e., $\bar{\beta}$, $\bar{\delta}$ and $p$) and the parameters of the
replacement process (i.e., $\mu$, $\sigma^{2}$ and $m^{2}$). Similarly to
(13), the bound is proportional to $\sigma^{2}+m^{2}$, but in this case the
rate of the Poisson process $\mu$ also appears as a multiplication factor. The
upper bound for the rate of convergence of the SIS epidemic in a closed system
determined by (7) appears in the denominator of (16), but this value is
increased by the rate of the Poisson process $\mu$. From (19), it can be seen
that a large value of $\mu$ will increase the rate of convergence of
$V(x(t))$, such that the convergence will be fast, but the asymptotic value
will also increase. When $\mu\to\infty$, the bound (16) is given by
$\sigma^{2}+m^{2}$, which coincides with the result of Proposition 3 since the
process will be characterized only by replacements.

Notice that the right-hand side in (17) corrresponds to the expectation of a
SDE of the form:

$dV=f(V)dt+g(V)dN_{t}.$ (20)

Unlike the SDE (9) associated to the original process, the mapping $f$ is not
time-varying and the Poisson process $N_{t}$ is homogeneous.

##  VI Second moment

Although the asymptotic behavior of the expected aggregate function
$\mathbb{E}\left[V(x(t))\right]$ provides useful information of a system
subject to the replacements of agents, it does not give any information about
the deviation of single realizations from the expected value. To this purpose,
it is necessary to analyze the variance of the aggregate function, denoted by
Var$(V(x(t)))$. In this case, we use Itô Lemma formulated for jump processes
of the form (20).

######  Lemma 1 (Itô Lemma [21])

For a continuous differentiable function $\phi:\mathbb{R}\to\mathbb{R}$:

$d\phi(V)=\langle\frac{d\phi}{dV},f(V)\rangle
dt+\left[\phi(V+g(V))-\phi(V)\right]dN_{t},$

where $f(V)$ and $g(V)$ are the mappings in the SDE (20) and $N_{t}$ is a
Poisson process.

We first analyze the change of $\mathbb{E}\left[V^{2}(x(t))\right]$ during a
replacement.

###### Lemma 2

During the replacement of an agent, the aggregate function $V(x(t))$ defined
in (3) satisfies:

$\displaystyle\mathbb{E}\left[V^{2}(x(t^{+}))\right]\\!-\\!\mathbb{E}\left[V^{2}(x(t^{-}))\right]\\!\leq\\!-\frac{2}{n_{0}}\mathbb{E}\left[V^{2}(x(t^{-}))\right]\\!+\\!\frac{\mathbb{E}\left[\Theta^{4}\right]}{n_{0}^{2}}$
$\displaystyle+\left(\frac{2\mathbb{E}\left[\Theta^{2}\right](n_{0}-1)+1}{n_{0}^{2}}\right)\mathbb{E}\left[V(x(t^{-}))\right],$

where $\Theta$ is the new value of a replaced agent according to Assumption 4.

###### Proof:

During the replacement of an agent $x_{j}$, the function $V^{2}$ satisfies:

$\displaystyle(V^{2})^{+}\\!$
$\displaystyle=\left(V^{-}+\frac{\Theta^{2}}{n_{0}}-\frac{(x_{j}^{2})^{-}}{n_{0}}\right)^{2}$
$\displaystyle=(V^{2})^{-}\\!+\\!\frac{\Theta^{4}}{n_{0}^{2}}\\!+\\!\frac{(x_{j}^{4})^{-}}{n_{0}^{2}}\\!+\\!\frac{2\Theta^{2}V^{-}}{n_{0}}\\!-\\!\frac{2(x_{j}^{2})^{-}V^{-}}{n_{0}}\\!-\\!\frac{2\Theta^{2}(x_{j}^{2})^{-}}{n_{0}^{2}},$

where $(V^{2})^{+}$ and $(V^{2})^{-}$ denote $V^{2}(x(t^{+}))$ and
$V^{2}(x(t^{-}))$ respectively. Then, the conditional expectation is given by:

$\displaystyle\mathbb{E}\left[\\!(V^{2})^{+}\\!|x^{-}\\!,\Theta\right]$
$\displaystyle=\frac{1}{n_{0}}\sum_{j=1}^{n_{0}}\Bigg{(}(V^{2})^{-}+\frac{\Theta^{4}}{n_{0}^{2}}+\frac{(x_{j}^{4})^{-}}{n_{0}^{2}}+\frac{2\Theta^{2}}{n_{0}}V^{-}$
$\displaystyle\quad-\frac{2(x_{j}^{2})^{-}}{n_{0}}V^{-}-\frac{2\Theta^{2}(x_{j}^{2})^{-}}{n_{0}^{2}}\Bigg{)}$
$\displaystyle\leq(V^{2})^{-}+\frac{\Theta^{4}}{n_{0}^{2}}+\frac{1}{n_{0}^{3}}\sum_{j=1}^{n_{0}}(x_{j}^{2})^{-}+\frac{2\Theta^{2}}{n_{0}}V^{-}$
$\displaystyle\quad-\frac{2(V^{2})^{-}}{n_{0}}-\frac{2\Theta^{2}}{n_{0}^{2}}V^{-}$
$\displaystyle=\\!\left(\frac{n_{0}\\!-\\!2}{n_{0}}\right)\\!\\!(V^{2})^{-}\\!+\\!\frac{\Theta^{4}}{n_{0}^{2}}\\!+\\!\left(\frac{1}{n_{0}^{2}}\\!+\\!\frac{2\Theta^{2}}{n_{0}}\\!-\\!\frac{2\Theta^{2}}{n_{0}^{2}}\right)\\!V^{-}.$

Finally, we compute the total expectation to get the desired result. ∎

######  Proposition 4 (Second moment)

Consider a SIS epidemic satisfying the stability condition (2). Under
Assumptions 3 and 4, the aggregate function $V(x(t))$ defined in (3)
satisfies:

$\limsup_{t\to\infty}\mathbb{E}\left[V^{2}(x(t))\right]\\!\leq\\!\\\
\frac{\mu\mathbb{E}\\!\\!\left[\Theta^{4}\right]\\!\\!(\mu\\!+\\!2n_{0}\bar{\delta}\\!-\\!2\bar{\beta}p(n_{0}\\!-\\!1))\\!+\\!\mu^{2}\mathbb{E}\\!\\!\left[\Theta^{2}\right]\\!\\!(2\mathbb{E}\\!\\!\left[\Theta^{2}\right]\\!\\!(n_{0}\\!-\\!1)\\!+\\!1)}{2n_{0}(\mu+2n_{0}\bar{\delta}-2\bar{\beta}p(n_{0}-1))^{2}}.$
(21)

###### Proof:

Thanks to Itô Lemma with $\phi(V)=V^{2}$ we obtain:

$\displaystyle dV^{2}$ $\displaystyle\leq\langle
2V,2(\beta\bar{\lambda}_{1}-\delta)V\rangle
dt+\left((V^{2})^{+}-(V^{2})^{-}\right)dN_{t}$
$\displaystyle=4(\beta\bar{\lambda}_{1}-\delta)V^{2}dt+\left((V^{2})^{+}-(V^{2})^{-}\right)dN_{t}.$
(22)

We take the expectation in (22) and we get:

$\displaystyle d\mathbb{E}\left[V^{2}\right]\\!\\!$
$\displaystyle\leq\\!4(\beta\bar{\lambda}_{1}\\!-\\!\delta)\mathbb{E}\left[V^{2}\right]\\!dt\\!+\\!\left(\mathbb{E}\left[(V^{2})^{+}\right]\\!-\\!\mathbb{E}\left[(V^{2})^{-}\right]\right)\\!\mu
dt$
$\displaystyle=2\left(\frac{2n_{0}\beta\bar{\lambda}_{1}-2n_{0}\delta-\mu}{n_{0}}\right)\mathbb{E}\left[V^{2}\right]dt$
$\displaystyle\quad\\!\\!+\mu\\!\\!\left(\frac{2\mathbb{E}\left[\Theta^{2}\right]\\!(n_{0}-1)+1}{n_{0}^{2}}\right)\\!\mathbb{E}\left[V\right]dt+\frac{\mu\mathbb{E}\left[\Theta^{4}\right]}{n_{0}^{2}}dt,$

where we use the Dominated Convergence Theorem to interchange the expectation
and the derivative of $V^{2}$ and we apply Lemma 2. Let us denote
$a_{V}\\!=\\!2\left(\frac{2n_{0}\beta\bar{\lambda}_{1}-2n_{0}\delta-\mu}{n_{0}}\right)$,
$b_{V}\\!=\\!\mu\left(\frac{2\mathbb{E}\left[\Theta^{2}\right](n_{0}-1)+1}{n_{0}^{2}}\right)$
and $c_{V}\\!=\\!\frac{\mu\mathbb{E}\left[\Theta^{4}\right]}{n_{0}^{2}}$, so
that we have

$\frac{d}{dt}\mathbb{E}\left[V^{2}\right]\leq
a_{V}\mathbb{E}\left[V^{2}\right]+b_{V}\mathbb{E}\left[V\right]+c_{V}.$ (23)

The solution of (23) is given by:

$\displaystyle\mathbb{E}\left[V^{2}\right]\\!$ $\displaystyle\leq
e^{a_{V}t}\mathbb{E}\left[V_{0}^{2}\right]+e^{a_{V}t}\int_{0}^{t}e^{-a_{V}\tau}\left(b_{V}\mathbb{E}\left[V(x(\tau))\right]+c_{V}\right)d\tau$
$\displaystyle=e^{a_{V}t}\mathbb{E}\left[V_{0}^{2}\right]\\!+b_{V}e^{a_{V}t}\\!\\!\\!\int_{0}^{t}\\!\\!e^{-a_{V}\tau}\mathbb{E}\left[V\right]d\tau\\!-\\!\frac{c_{V}}{a_{V}}\\!\left(1\\!-\\!e^{a_{V}t}\right).$

We recall that by (19), the function $\mathbb{E}\left[V\right]$ satisfies:

$\mathbb{E}\left[V(x(t))\right]\leq d_{V}e^{\frac{a_{V}}{2}t}+m_{V},$

where
$d_{V}\\!=\\!\mathbb{E}\left[V_{0}\right]-\frac{\mu(\sigma^{2}+m^{2})}{\mu+2n_{0}\delta-2n_{0}\beta\bar{\lambda}_{1}}$
and
$m_{V}\\!=\\!\frac{\mu(\sigma^{2}+m^{2})}{\mu+2n_{0}\delta-2n_{0}\beta\bar{\lambda}_{1}}$.
Then $\mathbb{E}\left[V^{2}\right]$ satisfies:

$\displaystyle\mathbb{E}\left[V^{2}\right]$ $\displaystyle\leq
e^{a_{V}t}\mathbb{E}\left[V_{0}^{2}\right]-\frac{c_{V}}{a_{V}}\left(1-e^{a_{V}t}\right)$
$\displaystyle\quad+b_{V}e^{a_{V}t}\int_{0}^{t}e^{-a_{V}\tau}\left(d_{V}e^{\frac{a_{V}}{2}\tau}+m_{V}\right)d\tau$
$\displaystyle=\left(\mathbb{E}\left[V_{0}^{2}\right]+\frac{c_{V}}{a_{V}}+\frac{2b_{V}d_{V}}{a_{V}}+\frac{b_{V}m_{V}}{a_{V}}\right)e^{a_{V}t}$
$\displaystyle\quad-\frac{2b_{V}d_{V}}{a_{V}}e^{\frac{a_{V}}{2}t}-\left(\frac{c_{V}}{a_{V}}+\frac{b_{V}m_{V}}{a_{V}}\right).$
(24)

Finally, since $a_{V}$ is negative, the result follows by taking the limit
$t\to\infty$ in (24) and using $\bar{\lambda}_{1}=p(n_{0}-1)$,
$\delta=\bar{\delta}$ and $\beta=n_{0}^{-1}\bar{\beta}$. ∎

The bound (21) is derived using (16), which corresponds to the asymptotic
behavior of $\mathbb{E}\left[V(x(t))\right]$. Unlike the previous result (16),
the value of (21) is also affected by the fourth moment of $\Theta$, while
(16) only depends on the second moment
($\mathbb{E}\left[\Theta^{2}\right]=\sigma^{2}+m^{2}$).

The previous bound on the second moment can also be used to estimate the
variance

$\displaystyle\limsup_{t\to\infty}\,\mathrm{Var}(V)$
$\displaystyle=\limsup_{t\to\infty}\left(\mathbb{E}\left[V^{2}\right]-\left(\mathbb{E}\left[V\right]\right)^{2}\right)$
$\displaystyle\leq\limsup_{t\to\infty}\mathbb{E}\left[V^{2}\right],$

although the result is likely to be conservative since the effect of
$-\left(\mathbb{E}\left[V(x(t))\right]\right)^{2}$ is neglected.

##  VII Numerical and Simulation Results

To illustrate the main results of this paper, we consider a SIS epidemic over
a graph with $n_{0}=50$ agents generated with a probability $p=0.5$. The
parameters of the epidemic are chosen as $\beta_{n}=\bar{\beta}/n=0.1/n$ and
$\delta_{n}=\bar{\delta}=1.5p\bar{\beta}$, and the Poisson processes for
arrivals, departures and replacements have a rate $\mu_{a}=\mu_{d}=\mu=7$. The
states of the arriving agents are generated with a uniform distribution with
$m=1/2$ and $\sigma^{2}=1/12$. Fig. 1 presents a realization of the stochastic
process. In the top plot, we show the evolution of $V(x(t))$ where the
aggregate function jumps during the occurrence of arrivals (dashed green line)
and departures (dash-dotted red line). In the bottom plot, we can see the
changes of the number of agents during the evolution of the epidemic, which
remains around the initial number $n_{0}=50$.

Figure 1: Sample realization of the aggregate function $V(x(t))$ and the
number of agents $n(t)$ for a graph with $n_{0}=50$ agents and probability
$p=0.5$. The SIS epidemic is considered with $\beta_{n}=0.1/n$ and
$\delta_{n}=0.15p$, and the rate of the Poisson processes is
$\mu_{a}=\mu_{d}=\mu=7$. The arrival and departure instants correspond to the
green and red lines respectively.

In Fig. 2, we present the computations of the moments of $V(x(t))$. The top
plot corresponds to the evolution of $\mathbb{E}\left[V(x(t))\right]$ where
the solid blue line is the original process, the dashed red line is the
approximate process and the dash-dotted yellow line corresponds to the upper
bound (16). Similarly, the middle plot shows the evolution of
$\mathbb{E}\left[V^{2}(x(t))\right]$ and the bottom plot presents the
evolution of Var$(V(x(t)))$. Notice that the behavior of
$\mathbb{E}\left[V(x(t))\right]$, $\mathbb{E}\left[V^{2}(x(t))\right]$ and
$\text{Var}(V(x(t)))$ of the process with arrivals/departures and the process
with replacements are really close.

Figure 2: Evolution of the moments of the aggregate function $V(x(t))$ for a
graph with $n_{0}=50$ agents and probability $p=0.5$. The SIS epidemic is
considered with $\beta=0.1/n$ and $\delta=0.15p$, and the rate of the Poisson
processes is $\mu_{a}=\mu_{d}=\mu=7$. In the top plot, the solid blue line and
dashed red line correspond to the estimation of
$\mathbb{E}\left[V(x(t))\right]$ for the original and approximate process
respectively, while the dash-dotted yellow line is the upper bound (16). In
the middle plot, the solid blue line and dashed red line correspond to the
estimation of $\mathbb{E}\left[V^{2}(x(t))\right]$ for the original and
approximate process respectively, while the dash-dotted yellow line is the
upper bound (21). In the bottom plot, the solid blue line and dashed red line
correspond to the estimation of $\,\mathrm{Var}(V(x(t)))$ for the original and
approximate process respectively, while the dash-dotted yellow line is the
upper bound (21). The simulated values were obtained considering 1000
realizations of the process.

##  VIII Conclusion

In this paper, we analyzed a SIS epidemic in open multi-agent systems by using
replacements to approximate a process with similar rates of arrivals and
departures, and a fixed network to approximate time-varying graphs that are
sampled from the same distribution. We defined an aggregate function to
analyze the epidemic and derived upper bounds for the asymptotic values of the
expectation and variance as functions of the parameters of the infection and
replacement processes. Through simulations we showed that the original process
(i.e., with arrivals and departures) and the approximate process (i.e., with
replacements) are really close, such that replacements can be used as a tool
to analyze real open systems.

For future work, we would like to consider the more general and challenging
case when arrivals and departures have different rates [26], which would not
preserve the initial number of agents in expectation. Finally, it would be
important to extend the results to process where the connections of new agents
are given by a more complex distribution. In this direction, graphons, which
have already been used to study epidemics in large networks [19], appear to be
a promising model to generate network topologies that are consistent across
time [27].

## References

  * [1] H. W. Hethcote, “The mathematics of infectious diseases,” _SIAM review_ , vol. 42, no. 4, pp. 599–653, 2000.
  * [2] I. Z. Kiss, J. C. Miller, and P. L. Simon, _Mathematics of Epidemics on Networks_. Cham, Switzerland: Springer, 2017.
  * [3] C. Nowzari, V. M. Preciado, and G. J. Pappas, “Analysis and control of epidemics: A survey of spreading processes on complex networks,” _IEEE Control Systems Magazine_ , vol. 36, no. 1, pp. 26–46, 2016.
  * [4] T. Alamo, P. Millán, D. G. Reina, V. M. Preciado, and G. Giordano, “Challenges and future directions in pandemic control,” _IEEE Control Systems Letters_ , vol. 6, pp. 722–727, 2022.
  * [5] M. Tizzoni, P. Bajardi, A. Decuyper, G. Kon Kam King, C. M. Schneider, V. Blondel, Z. Smoreda, M. C. González, and V. Colizza, “On the use of human mobility proxies for modeling epidemics,” _PLoS computational biology_ , vol. 10, no. 7, p. e1003716, 2014.
  * [6] S. Hazarie, D. Soriano-Paños, A. Arenas, J. Gómez-Gardeñes, and G. Ghoshal, “Interplay between population density and mobility in determining the spread of epidemics in cities,” _Communications Physics_ , vol. 4, no. 1, pp. 1–10, 2021.
  * [7] M. Ogura and V. M. Preciado, “Stability of spreading processes over time-varying large-scale networks,” _IEEE Transactions on Network Science and Engineering_ , vol. 3, no. 1, pp. 44–57, 2016.
  * [8] P. E. Paré, C. L. Beck, and A. Nedić, “Epidemic processes over time-varying networks,” _IEEE Transactions on Control of Network Systems_ , vol. 5, no. 3, pp. 1322–1334, 2018.
  * [9] J. Leitch, K. A. Alexander, and S. Sengupta, “Toward epidemic thresholds on temporal networks: a review and open questions,” _Applied Network Science_ , vol. 4, pp. 1–21, 2019.
  * [10] L. Zino and M. Cao, “Analysis, prediction, and control of epidemics: A survey from scalar to dynamic network models,” _IEEE Circuits and Systems Magazine_ , vol. 21, no. 4, pp. 4–23, 2021.
  * [11] V. Abhishek and V. Srivastava, “SIS epidemic spreading under multi-layer population dispersal in patchy environments,” _IEEE Transactions on Control of Network Systems_ , 2023.
  * [12] J. M. Hendrickx and S. Martin, “Open multi-agent systems: Gossiping with random arrivals and departures,” in _2017 IEEE 56th Annual Conference on Decision and Control (CDC)_ , 2017, pp. 763–768.
  * [13] M. Franceschelli and P. Frasca, “Stability of open multiagent systems and applications to dynamic consensus,” _IEEE Transactions on Automatic Control_ , vol. 66, no. 5, pp. 2326–2331, 2021.
  * [14] R. Vizuete, C. Monnoyer de Galland, P. Frasca, E. Panteley, and J. M. Hendrickx, “Trends and questions in open multi-agent systems,” in _Hybrid and Networked Dynamical Systems_. Springer, 2024.
  * [15] T. Carletti, D. Fanelli, A. Guarino, F. Bagnoli, and A. Guazzini, “Birth and death in a continuous opinion dynamics model: The consensus case,” _The European Physical Journal B_ , vol. 64, pp. 285–292, 2008.
  * [16] J. Török, G. Iñiguez, T. Yasseri, M. San Miguel, K. Kaski, and J. Kertész, “Opinions, conflicts, and consensus: Modeling social dynamics in a collaborative environment,” _Physical Review Letters_ , vol. 110, no. 8, p. 088701, 2013.
  * [17] A. Lajmanovich and J. A. Yorke, “A deterministic model for gonorrhea in a nonhomogeneous population,” _Mathematical Biosciences_ , vol. 28, no. 3-4, pp. 221–236, 1976.
  * [18] A. Khanafer, T. Başar, and B. Gharesifard, “Stability properties of infection diffusion dynamics over directed networks,” in _53rd IEEE Conference on Decision and Control_ , 2014, pp. 6215–6220.
  * [19] R. Vizuete, P. Frasca, and F. Garin, “Graphon-based sensitivity analysis of SIS epidemics,” _IEEE Control Systems Letters_ , vol. 4, no. 3, pp. 542–547, 2020.
  * [20] S. Gao and P. E. Caines, “Spectral representations of graphons in very large network systems control,” in _2019 IEEE 58th conference on decision and Control (CDC)_. IEEE, 2019, pp. 5068–5075.
  * [21] R. Brockett, “Stochastic control,” _Lecture Notes, Harvard University_ , 2009\.
  * [22] R. Brockett, W. Gong, and Y. Guo, “Stochastic analysis for fluid queueing systems,” in _Proceedings of the 38th IEEE Conference on Decision and Control (Cat. No.99CH36304)_ , vol. 3, 1999, pp. 3077–3082 vol.3.
  * [23] M. Jacobsen, _Point Process Theory and Applications: Marked Point and Piecewise Deterministic Processes_. Birkhäuser Boston, 2006.
  * [24] F. Chung and M. Radcliffe, “On the spectra of general random graphs,” _The Electronic Journal of Combinatorics_ , vol. 18, no. 1, p. 215, 2011.
  * [25] N. Privault, _Stochastic Finance: An Introduction with Market Examples_. CRC Press, 2013.
  * [26] C. Monnoyer de Galland, S. Martin, and J. M. Hendrickx, “Modelling gossip interactions in open multi-agent systems,” _arXiv preprint arXiv:2009.02970_ , 2020.
  * [27] R. Vizuete, “Contributions to open multi-agent systems: consensus, optimization and epidemics,” Ph.D. dissertation, Université Paris-Saclay, 2022.