A Quantum-Compliant Formulation for Network Epidemic Control

TL;DR

This paper presents a quantum-compatible optimization framework for controlling epidemic spread on networks through mobility bans, balancing healthcare and social costs, and leveraging quantum computing for NP-hard problem solving.

Contribution

It introduces a QUBO formulation for epidemic control via mobility bans, enabling quantum computing approaches to address the NP-hard optimization problem.

Findings

QUBO formulation effectively models the control problem.

Quantum computing can potentially solve the NP-hard problem efficiently.

Framework balances epidemic control with social and economic costs.

Abstract

We deal with controlling the spread of an epidemic disease on a network by isolating one or multiple locations by banning people from leaving them. To this aim, we build on the susceptible-infected-susceptible and the susceptible-infected-removed discrete-time network models, encapsulating a control action that captures mobility bans via removing links from the network. Then, we formulate the problem of optimally devising a control policy based on mobility bans that trades-off the burden on the healthcare system and the social and economic costs associated with interventions. The binary nature of mobility bans hampers the possibility to solve the control problem with standard optimization methods, yielding a NP-hard problem. Here, this is tackled by deriving a Quadratic Unconstrained Binary Optimization (QUBO) formulation of the control problem, and leveraging the growing potentialities of quantum computing to efficiently solve it.

Figures

Tables

Table 1. Table I: Performance obtained with different methods in terms of reduction of the peak (p) and average infected/day (a) with respect to the uncontrolled dynamics. For each scenario, the best performance is highlighted using bold font.

Data	n	SA		TS		GA		QA
		p [%]	a [%]	p [%]	a [%]	p [%]	a [%]	p [%]	a [%]
2020/03/08	21	12.40	3.48	12.40	3.48	12.40	3.48	12.40	2.56
2020/03/10	21	9.10	2.63	9.10	2.63	9.10	2.63	9.10	2.20
2020/10/29	21	12.77	1.93	12.77	1.93	12.77	1.93	12.77	3.00
2021/12/20	21	13.73	1.72	13.73	1.72	13.73	1.72	13.73	10.05
2020/03/08	107	4.93	17.96	4.93	17.96	-1.48	14.31	3.61	16.66
2020/03/10	107	3.37	14.79	3.37	14.79	-2.61	11.99	2.15	13.82
2020/10/29	107	15.31	9.02	15.31	9.02	14.02	7.74	15.29	9.02
2021/12/20	107	20.70	7.54	20.70	7.54	20.31	7.34	20.70	7.54

Table 2. Table II: Average solving time for the different methods. For each scenario, the best performance is highlighted using bold font.

Data	n	SA [$\mathrm{ms}$]	TS [$\mathrm{ms}$]	GA [$\mathrm{ms}$]	QA [$\mathrm{ms}$]
2020/03/08	21	2.07	2112.88	55.92	0.04
2020/03/10	21	1.94	2103.59	49.50	0.04
2020/10/29	21	1.93	2110.27	51.32	0.04
2021/12/20	21	1.89	2119.29	68.24	0.04
2020/03/08	107	784.40	2124.64	508.54	0.04
2020/03/10	107	1072.64	2126.41	440.78	0.04
2020/10/29	107	824.93	2124.39	495.81	0.04
2021/12/20	107	1306.50	2116.03	371.79	0.04

Equations

$$x_i(t+1)=(1-\mu )x_i(t)+\frac{\lambda }{n_i}\left(n_i-x_i(t)\right){\alpha }_i(t),$$

$${\alpha }_i(t)=x_i(t)+{\sum }_{j\in \mathcal{L}}{A}_{ij}{x}_j(t)$$

$$\begin{array}[]{l}\displaystyle x_{i}(t+1)=(1-\mu)x_{i}(t)+\frac{\lambda}{n_{i}}\left(n_{i}-x_{i}(t)-y_{i}(t)\right)\alpha_{i}(t),\\ y_{i}(t+1)=y_{i}(t)+\mu x_{i}(t),\end{array}$$

$$u_{i}=\left\{\begin{array}[]{ll}1&\text{if location $i$ is isolated,}\\ 0&\text{otherwise.}\end{array}\right.$$

$${\alpha }_i(\mathbf{u},t)=x_i(t)+(1-u_i){\sum }_{j\in \mathcal{L}}{A}_{ij}{x}_j(t).$$

$$f(\mathbf{x},\mathbf{u})={\sum }_{i\in \mathcal{L}}{\sum }_{t=1}^T{x}_i(t)+\gamma {\sum }_{i\in \mathcal{L}}{n}_i{u}_i,$$

$$\begin{array}[]{rl}\boldsymbol{u^{*}}=&\arg\min_{\boldsymbol{x},\boldsymbol{u}}f(\boldsymbol{x},\boldsymbol{u})\\ \text{s.t.}&Eq.\leavevmode\nobreak\ (\ref{eq:sis})\text{ or }Eq.\leavevmode\nobreak\ (\ref{eq:sir})\text{, and }Eq.\leavevmode\nobreak\ (\ref{eq:alpha_control}),\end{array}$$

$${\mathbf{z}}^{\ast }=\arg \underset{\mathbf{u}\in \{0,1{\}}^M}{\min }{\sum }_{i\in \mathcal{L}}{P}_i{z}_i+{\sum }_{i\in \mathcal{L},j\in \mathcal{L}\setminus \{i\}}{Q}_{ij}{z}_i{z}_j$$

$$\begin{array}[]{l}P_{i}=\displaystyle\lambda\Big{(}1\hskip-1.13791pt-\hskip-1.13791pt\frac{\bar{x}_{i}}{n_{i}}\Big{)}\sum_{j\in\mathcal{L}}A_{ij}\bar{x}_{j}+\lambda^{2}(1\hskip-1.13791pt-\hskip-1.13791pt\frac{\bar{x}_{i}}{n_{i}}\Big{)}^{2}\Big{(}\sum_{j\in\mathcal{L}}A_{ij}\bar{x}_{j}\Big{)}^{2}\\ \,\,\,\displaystyle+\lambda^{2}\Big{[}1-\frac{\bar{x}_{i}}{n_{i}}\Big{(}1-\mu+\lambda\Big{(}1-\frac{\bar{x}_{i}}{n_{i}}\Big{)}\Big{]}\Big{(}1-\frac{\bar{x}_{i}}{n_{i}}\Big{)}\sum_{j\in\mathcal{L}}A_{ij}\bar{x}_{j}\\ \,\,\,\displaystyle+\lambda\Big{(}1-\frac{\bar{x}_{i}}{n_{i}}\Big{)}\Big{[}1-\mu+\lambda\Big{(}1-\frac{\bar{x}_{i}}{n_{i}}\Big{)}\Big{]}\bar{x}_{i}\sum_{j\in\mathcal{L}}A_{ij}\bar{x}_{j}-\gamma n_{i}\\ \,\,\,\displaystyle+\lambda\Big{[}1-\frac{\bar{x}_{i}}{n_{i}}\Big{(}1-\mu+\lambda\Big{(}1-\frac{\bar{x}_{i}}{n_{i}}\Big{)}\Big{)}-\frac{\lambda}{n_{i}}\Big{(}1-\frac{\bar{x}_{i}}{n_{i}}\Big{)}\\ \qquad\displaystyle\cdot\sum\nolimits_{j\in\mathcal{L}}A_{ij}\bar{x}_{j}\Big{]}\sum\nolimits_{j\in\mathcal{L}}A_{ij}\Big{[}1-\mu+\lambda\Big{(}1-\frac{\bar{x}_{j}}{n_{j}}\Big{)}\Big{]}\bar{x}_{j}\end{array}$$

$$\begin{array}[]{l}Q_{ij}=\displaystyle\lambda^{2}A_{ij}\Big{[}1-\frac{\bar{x}_{i}}{n_{i}}\Big{(}1-\mu+\lambda\Big{(}1-\frac{\bar{x}_{i}}{n_{i}}\Big{)}\Big{)}\\ \displaystyle-\frac{\lambda}{n_{i}}\Big{(}1-\frac{\bar{x}_{i}}{n_{i}}\Big{)}\hskip-2.84544pt\sum\nolimits_{k\in\mathcal{L}}\hskip-8.5359ptA_{ik}\bar{x}_{k}\Big{]}\hskip-2.84544pt\Big{(}1-\frac{\bar{x}_{j}}{n_{j}}\Big{)}\sum\nolimits_{k\in\mathcal{L}}\hskip-8.5359ptA_{jk}\bar{x}_{k}.\end{array}$$

$$x_{i}(1)=(1-\mu)\bar{x}_{i}+\lambda\Big{(}1-\frac{\bar{x}_{i}}{n_{i}}\Big{)}\Big{(}\bar{x}_{i}+(1-u_{i})\sum\nolimits_{j\in\mathcal{L}}A_{ij}\bar{x}_{j}\Big{)}.$$

$$\begin{array}[]{l}x_{i}(2)=\displaystyle(1-\mu)x_{i}(1)+\lambda\Big{(}1-\frac{x_{i}(1)}{n_{i}}\Big{)}\\ \qquad\qquad\cdot\displaystyle\Big{[}x_{i}(1)+(1-u_{i})\sum\nolimits_{j\in\mathcal{L}}A_{ij}x_{j}(1)\Big{]}.\end{array}$$

$$\begin{array}[]{l}x_{i}(2)=\displaystyle(1-\mu)^{2}\bar{x}_{i}+(1-\mu)\lambda\Big{(}1-\frac{x_{i}}{n_{i}}\Big{)}\bar{x}_{i}\\ +\displaystyle(1-\mu)\lambda(1-u_{i})\Big{(}1-\frac{x_{i}}{n_{i}}\Big{)}\sum\nolimits_{j\in\mathcal{L}}A_{ij}\bar{x}_{j}\\ +\displaystyle\lambda\Big{[}1\hskip-1.13791pt-\hskip-1.13791pt\frac{\bar{x}_{i}}{n_{i}}\Big{(}1\hskip-1.13791pt-\hskip-1.13791pt\mu\hskip-1.13791pt+\hskip-1.13791pt\lambda\Big{(}1\hskip-1.13791pt-\hskip-1.13791pt\frac{x_{i}}{n_{i}}\Big{)}\Big{)}\Big{]}\Big{(}1\hskip-1.13791pt-\hskip-1.13791pt\mu\hskip-1.13791pt+\hskip-1.13791pt\lambda\Big{(}1\hskip-1.13791pt-\hskip-1.13791pt\frac{x_{i}}{n_{i}}\Big{)}\Big{)}\bar{x}_{i}\\ +\displaystyle\lambda^{2}\Big{[}1-\frac{\bar{x}_{i}}{n_{i}}\Big{(}1-\mu+\lambda\Big{(}1-\frac{x_{i}}{n_{i}}\Big{)}\Big{)}\Big{]}(1-u_{i})\Big{(}1-\frac{x_{i}}{n_{i}}\Big{)}\\ \displaystyle\cdot\sum\nolimits_{j\in\mathcal{L}}A_{ij}\bar{x}_{j}+(1-u_{i})\lambda^{2}\Big{(}1-\frac{x_{i}}{n_{i}}\Big{)}^{2}\Big{(}\sum\nolimits_{j\in\mathcal{L}}A_{ij}\bar{x}_{j}\Big{)}^{2}\\ +\displaystyle(1-u_{i})\lambda\Big{(}1-\frac{\bar{x}_{i}}{n_{i}}\Big{)}\Big{(}1-\mu+\lambda\Big{(}1-\frac{x_{i}}{n_{i}}\Big{)}\Big{)}\bar{x}_{i}\sum_{j\in\mathcal{L}}A_{ij}\bar{x}_{j}\\ +\displaystyle(1-u_{i})\lambda\Big{[}1-\frac{\bar{x}_{i}}{n_{i}}\Big{(}1-\mu+\lambda\Big{(}1\hskip-1.13791pt-\hskip-1.13791pt\frac{\bar{x}_{i}}{n_{i}}\Big{)}\Big{)}+\lambda\Big{(}1\hskip-1.13791pt-\hskip-1.13791pt\frac{\bar{x}_{i}}{n_{i}}\Big{)}\\ \displaystyle\cdot\sum\nolimits_{j\in\mathcal{L}}A_{ij}\bar{x}_{j}\Big{]}\Big{[}\sum\nolimits_{j\in\mathcal{L}}A_{ij}\Big{(}1-\mu+\lambda\Big{(}1-\frac{\bar{x}_{j}}{n_{j}}\Big{)}\Big{)}\bar{x}_{j}\Big{]}+\\ \displaystyle\lambda^{2}\hskip-1.13791pt\Big{[}1\hskip-1.13791pt-\hskip-1.13791pt\frac{\bar{x}_{i}}{n_{i}}\Big{(}1\hskip-1.13791pt-\hskip-1.13791pt\mu\hskip-1.13791pt+\hskip-1.13791pt\lambda\Big{(}1\hskip-1.13791pt-\hskip-1.13791pt\frac{\bar{x}_{i}}{n_{i}}\Big{)}\hskip-1.13791pt\Big{)}\hskip-1.13791pt+\hskip-1.13791pt\lambda\Big{(}1\hskip-1.13791pt-\hskip-1.13791pt\frac{\bar{x}_{i}}{n_{i}}\Big{)}\sum\nolimits_{j\in\mathcal{L}}\hskip-1.13791ptA_{ij}\bar{x}_{j}\hskip-1.13791pt\Big{]}\\ \displaystyle\cdot\sum\nolimits_{j\in\mathcal{L}}A_{ij}(1-u_{i})(1-u_{j})\Big{(}1-\frac{\bar{x}_{j}}{n_{j}}\Big{)}\sum\nolimits_{k\in\mathcal{L}}A_{jk}\bar{x}_{k},\end{array}$$

$$\begin{array}[]{l}P_{i}=\displaystyle\lambda\Big{(}1-\frac{\bar{x}_{i}+\bar{y}_{i}}{n_{i}}\Big{)}\sum\nolimits_{j\in\mathcal{L}}A_{ij}\bar{x}_{j}-\gamma n_{i}\\ \,\,\,\displaystyle+\lambda^{2}(1-\frac{\bar{x}_{i}+\bar{y}_{i}}{n_{i}}\Big{)}^{2}\Big{(}\sum\nolimits_{j\in\mathcal{L}}A_{ij}\bar{x}_{j}\Big{)}^{2}\\ \,\,\,\displaystyle+\lambda^{2}\Big{[}1-\frac{\bar{y}_{i}}{n_{i}}-\frac{\bar{x}_{i}}{n_{i}}\Big{(}1-2\mu+\lambda\Big{(}1-\frac{\bar{x}_{i}+\bar{y}_{i}}{n_{i}}\Big{)}\Big{]}\\ \displaystyle\quad\cdot\Big{(}1-\frac{\bar{x}_{i}+\bar{y}_{i}}{n_{i}}\Big{)}\hskip-2.84544pt\sum\nolimits_{j\in\mathcal{L}}\hskip-5.69046ptA_{ij}\bar{x}_{j}+\lambda\Big{(}1-\frac{\bar{x}_{i}+\bar{y}_{i}}{n_{i}}\Big{)}\\ \displaystyle\quad\cdot\Big{[}1-2\mu+\lambda\Big{(}1-\frac{\bar{x}_{i}+\bar{y}_{i}}{n_{i}}\Big{)}\Big{]}\bar{x}_{i}\sum\nolimits_{j\in\mathcal{L}}A_{ij}\bar{x}_{j}\\ \,\,\,\displaystyle+\lambda\Big{[}1-\frac{\bar{y}_{i}}{n_{i}}-\frac{\bar{x}_{i}}{n_{i}}\Big{(}1-2\mu+\lambda\Big{(}1-\frac{\bar{x}_{i}+\bar{y}_{i}}{n_{i}}\Big{)}\Big{)}\\ \,\,\,\displaystyle-\frac{\lambda}{n_{i}}\Big{(}1-\frac{\bar{x}_{i}+\bar{y}_{i}}{n_{i}}\Big{)}A_{ij}\bar{x}_{j}\Big{]}\\ \displaystyle\quad\cdot\sum\nolimits_{j\in\mathcal{L}}A_{ij}\Big{[}1-\mu+\lambda\Big{(}1-\frac{\bar{x}_{j}+\bar{y}_{j}}{n_{j}}\Big{)}\Big{]}\bar{x}_{j}\end{array}$$

$$\begin{array}[]{l}Q_{ij}=\displaystyle\lambda^{2}A_{ij}\Big{[}1-\frac{\bar{y}_{i}}{n_{i}}-\frac{\bar{x}_{i}}{n_{i}}\Big{(}1-2\mu+\lambda\Big{(}1-\frac{\bar{x}_{i}}{n_{i}}-\frac{\bar{y}_{i}}{n_{i}}\Big{)}\Big{)}\\ \,\,\,\displaystyle-\frac{\lambda}{n_{i}}\Big{(}1-\frac{\bar{x}_{i}+\bar{y}_{i}}{n_{i}}\Big{)}\sum_{k\in\mathcal{L}}A_{ik}\bar{x}_{k}\Big{]}\Big{(}1-\frac{\bar{x}_{j}+\bar{y}_{i}}{n_{j}}\Big{)}\sum_{k\in\mathcal{L}}A_{jk}\bar{x}_{k}.\end{array}$$

Full text

I Introduction

The recent COVID-19 health crisis has demonstrated how mathematical modeling and control-theoretic techniques can be powerful assets in the design, evaluation, and optimization of interventions policies during an epidemic outbreak [1, 2, 3]. Of particular interest are network epidemic models, thanks to their ability to capture the complex and heterogeneous patterns of contagion across a geographic region [4, 5, 6, 7, 8].

In this context, different questions have been investigated, including how to distribute drugs [9, 10, 11, 12, 13], plan a vaccination campaign [14, 15], and how guide a collective behavioral response [16, 17, 18, 19, 20]. Among this broad array of questions, a problem of paramount importance is to design intervention policies to control an epidemic outbreak, when no pharmaceutical interventions are available. In this scenario, control actions must rely on non-pharmaceutical interventions, including implementing lockdowns and travel bans.

Controlling network epidemic processes via non-pharmaceutical interventions is typically a complex problem from the computational point of view. In fact, non-pharmaceutical intervention policies are typically of discrete nature. For instance, deciding whether banning mobility between two locations in a geographic network or not is a binary decision variable. Similar, lockdown can be implemented at different (but finite) levels. Hence, the mathematical formulation of optimal control problems for such network epidemic models typically yields binary (or integer) optimization problems that are difficult to be solved in an exact manner with classical optimization methods [21]. Indeed, several heuristics have been developed, especially when one is interested in guaranteeing long-term eradication of the disease [22], but only a minority of the approaches focus on the transient behavior of the epidemic outbreak, limiting their applicability in real-life applications.

This calls for the exploration of novel optimization techniques to solve optimal control problems for the transient behavior of network epidemic processes. Traditional optimization approaches, while powerful, often face scalability limitations and computational bottlenecks when applied to such problems, especially under time constraints or in large-scale networks. In these regards, a promising research field that is recently emerging for the solution of complex and previously intractable optimization problems is quantum computing [23, 24, 25]. Quantum algorithms, such as the Quantum Approximate Optimization Algorithm (QAOA), Grover Adaptive Search (GAS), or quantum annealing (QA) [26, 27, 28], offer new computational paradigms that exploit quantum mechanical principles like superposition and entanglement to explore large solution spaces more efficiently than classical methods. In this work, we consider quantum annealers as the target quantum hardware, as current quantum circuit-based devices are too noisy to reliably execute optimization algorithms. Moreover, simulating such devices on classical hardware imposes severe limitations on the maximum problem size that can be explored.

A common feature of these approaches is their ability to efficiently solve problems that are formulated as Quadratic Unconstrained Binary Optimization (QUBO) problems [28]. The QUBO formulation acts as a unifying framework, allowing a wide class of optimization problems —including those arising in nonlinear model predictive control [29], scheduling [30], and resource allocation [31]— to be mapped onto hardware-compliant models that are solvable by quantum or quantum-inspired devices.

Motivated by these emerging techniques, we take a step towards their integration in the problem of controlling network epidemics. We consider the two most fundamental models of epidemic progression on networks [7]: susceptible–infected–susceptible (SIS) and susceptible–infected–removed (SIR) models, which capture diseases that do not provide immunity after recovery —e.g., most sexually transmitted diseases (STIs) [32]— and diseases that provide permanent (or long-lasting) immunity, respectively. Then, we encapsulate a control action into the model, consisting in isolating one or multiple nodes and we formulate an optimal control problem for the transient behavior of network epidemic processes in order to design the control action in an effective way, trading-off healthcare costs associated to the spread of the disease and the social and economic impact of node isolation.

Besides formulating the control problem, our contribution is twofold. First, we derive a QUBO formulation for our optimal control problem. Second, leveraging this formulation, we explore the potentiality of quantum computing in solving the control problem. In particular, we consider a case study based on the spread of an infectious disease in Italy, with realistic model parameters and network structure. We demonstrate how the control problem, formulated as a QUBO problem, can be efficiently addressed using a quantum annealer [26, 27], outperforming classical methods in terms of computational time, with comparable solutions quality.

II Controlled Network Epidemic Models

We denote by $\mathbb{R}$, $\mathbb{R}_{\geq 0}$, $\mathbb{R}_{>0}$, $\mathbb{Z}_{\geq 0}$, and $\mathbb{Z}_{>0}$ the real, real nonnegative, strictly positive real, positive integer, and strictly positive integer numbers, respectively. Given $n,m\in\mathbb{Z}_{>0}$, bold lowercase font denote a vector $\boldsymbol{x}\in\mathbb{R}^{n}$, with $x_{i}$ its $i$th entry; and bold capital font denote a matrix $\boldsymbol{A}\in\mathbb{R}^{n\times m}$, with $A_{ij}$ the $j$th entry of its $i$th row. The all-$1$ vector and the identity matrix are denoted by $\boldsymbol{1}$ and $\boldsymbol{I}$, respectively.

II-A Network population model

We consider $n$ individuals, partitioned into $M$ locations, denoted by $\mathcal{L}:=\{1,\dots,M\}$. In location $i$, the total population is equal to $n_{i}$ individuals. In this paper, we will focus on large populations. Hence, for the sake of simplicity, it is reasonable to approximate the population as a continuum, i.e., $n_{i}\in\mathbb{R}_{>0}$. Individuals from different locations can come in contact due to their mobility. To capture this feature, we introduce a weighted graph $\mathcal{G}=(\mathcal{L},\boldsymbol{A})$, where the weight matrix $\boldsymbol{A}\in\mathbb{R}_{\geq 0}^{M\times M}$ measures the level of interactions between individuals in different locations: $A_{ij}$ represents the amount of such interaction (relative to the interactions within each location). While not strictly necessary, it is reasonable to assume that $A_{ij}\in[0,1]$. Furthermore, we assume that $\mathcal{G}$ has no self-loops, yielding $A_{ii}=0$ for all $i\in\mathcal{L}$.

II-B Network SIS and SIR epidemic models

We consider two fundamental models of epidemic progression: the SIS and the SIR models [5, 7]. In both models, individuals can be characterized by two possible health conditions: susceptible (S) to the disease, or infected (I) with the disease. The two models differ in whether individuals acquire immunity after recovery. In the SIS model, no immunity is acquired, and recovered individuals become immediately susceptible to the disease again (e.g., many STIs [32]). In the SIR model, instead, individuals become permanently immune after recovery. This is a good proxy also for scenarios where immunity wanes at a slower time-scale than the one of an epidemic wave (e.g., COVID-19). For the SIR model, a third health state, termed removed (R), is used to represent this health condition.

For each location $i\in\mathcal{L}$, let us denote by $x_{i}(t)\in[0,n_{i}]$ the number of infected individuals at time $t\in\mathbb{Z}_{\geq 0}$ in location $i\in\mathcal{L}$. In the network SIS model, the number of susceptible individuals is equal to $n_{i}-x_{i}(t)$. Hence, the health state of the whole system is fully determined by the state vector $\boldsymbol{x}(t)=[x_{1}(t),\dots,x_{M}(t)]^{\top}$. For a generic location $i$, the number of infected individuals $x_{i}(t)$ is updated at each discrete time step according to two contrasting mechanisms: i) recovery, which is regulated by the recovery rate $\mu\in\mathbb{R}_{\geq 0}$, whereby, at each discrete time step, a fraction $\mu$ of the infected individuals recover; ii) contagion, whereby susceptible individuals who interact with infected individuals (in the same location or in other locations) become infected with infection rate $\lambda\in\mathbb{R}_{\geq 0}$. These two contrasting mechanisms yield the following recursive equation:

[TABLE]

for all $i\in\mathcal{L}$, where

[TABLE]

quantifies the infection force in location $i$ at time $t$. Eq. (2) comprises two terms: the first captures the contribution of infected individuals in the same location, the second term captures the impact of mobility to different locations.

The value of the parameters depends on the time-step. The recovery rate $\mu$ can be interpreted as the fraction of infected individuals who recover in a time step or (equivalently) the inverse of the mean duration of the disease (in time-steps).

Remark 1**.**

The infection rate $\lambda$ is related to the well-known concept of the basic reproduction number, whereby one can express $\lambda\propto R_{0}\mu$, where the proportionality coefficient is the largest eigenvalue of matrix $\boldsymbol{A}+\boldsymbol{I}$ [7]. Hence, both parameters $\mu$ and $\lambda$ scale with the duration of a time-step.

A consequence is that, by setting a sufficiently small time-step, it is always possible to guarantee that Eq. (1) is well-defined, as proved in the following statement.

Proposition 1**.**

If \lambda\leq\max_{i\in\mathcal{L}}\big{(}1+\sum_{j\in\mathcal{L}}A_{ij}\frac{n_{j}}{n_{i}}\big{)}^{-1}, then $\mathcal{D}=\prod_{i\in\mathcal{L}}[0,n_{i}]$ is positively invariant under Eq. (1).

Proof.

Consider a generic $i\in\mathcal{L}$ and $t\geq 0$. From Eq. (1), if $x_{i}(t)\in[0,n_{i}]$, then $x_{i}(t+1)\geq(1-\mu)x_{i}(t)\geq 0$. On the other hand, since \lambda\leq\big{(}1+\sum_{j\in\mathcal{L}}A_{ij}\frac{n_{j}}{n_{i}}\big{)}^{-1}, then $\alpha_{i}(t)\leq n_{i}+\sum_{j\in\mathcal{L}}A_{ij}n_{j}\leq\frac{n_{i}}{\lambda}$, implying $x_{i}(t+1)\leq(1-\mu)x_{i}(t)+(n_{i}-x_{i}(t))\leq n_{i}$, yielding the claim. ∎

In the SIR model, for each location, a second state variable should be defined to account for the number of removed individuals, $y_{i}(t)\in[0,n_{i}]$. In this case, individuals become removed after recovery and they cannot contract the disease again. Hence, the number of susceptible individuals is equal to $n_{i}-x_{i}(t)-y_{i}(t)$. In this scenario, the state of the system is characterized by the two vectors $\boldsymbol{x}(t)$ and $\boldsymbol{y}(t)=[y_{1}(t),\dots,y_{M}(t)]^{\top}$, whose components are updated according to the following dynamics:

[TABLE]

for all $i\in\mathcal{L}$, where the first equation differs from Eq. (1) only by the definition of susceptible population, which is now equal to $n_{i}-x_{i}(t)-y_{i}(t)$, and by the presence of an additional equation to account for recovered individuals who become removed. The observations in Remark 1 and Proposition 1 hold true also for Eq. (3), implying that a thoughtful choice of the time-step guarantees that the discrete-time network SIR model is well-defined.

II-C Control

In this paper, we focus on epidemic diseases for which no vaccines or effective pharmaceutical interventions are available. This is the case, e.g., of many STIs or for newly emerged diseases, as it was for the first wave of the COVID-19 pandemic. In such a scenario, the only intervention policies that can be enacted are non-pharmaceutical. Besides incentivizing the use of personal protection equipments (e.g., condoms for STIs or face masks for air-borne diseases), policy makers can enact intervention policies based on isolating one or multiple locations, not allowing the population of a location to leave that location. These interventions were implemented, e.g., during the first waves of the COVID-19 pandemic in many countries.

Technically, these interventions can be encapsulated into the network epidemic model by defining a set of binary control variables $\boldsymbol{u}=[u_{1},\dots,u_{M}]^{\top}\in\{0,1\}^{M}$, each entry associated with a location, such that

[TABLE]

Then, for a controlled network epidemic model, the infection force in Eq. (2) is replaced by the following expression, which accounts for isolation of one or multiple locations:

[TABLE]

In other words, if location $i$ is isolated, individuals in location $i$ are not allowed to leave their location, and so they can have interactions (and thus new contagions) only within their location. Since the control input can only reduce the infection force, the result in Proposition 1 remains valid, guaranteeing well-definedness of the controlled dynamics.

Remark 2**.**

In a more general setting, the control action in Eq. (4) can be time-varying, i.e., having $\boldsymbol{u}(t)$. However, as we shall see in the following, the time-varying nature of interventions can be incorporated into the framework by defining an optimal control problem over a fixed time-window, within a rolling horizon scheme, ultimately inducing a nonlinear model-predictive control scheme [33].

III Problem Statement

Designing optimal control policies using the interventions described in Eq. (4) is nontrivial due to the inherent nonlinearity of the dynamical system and the discrete nature of the control actions, which naturally induce a binary optimization problem, which is often an instance of an NP-hard problem [34]. Indeed, if one is interested only in the long-term behavior, it is known that the local stability of the disease-free equilibrium of the system is determined by the largest eigenvalue (in modulus) of the matrix $\boldsymbol{I}+\operatorname{\text{diag}}(1-\boldsymbol{u})\boldsymbol{A}$. Hence, optimizing the long-term behavior of the system can be mapped into the graph-theoretic problem of minimizing the spectral radius of a (weighted) adjacency matrix by edge removal, which is NP-hard [21, 35].

However, even the heuristics obtained with these graph-theoretic methods might be sub-optimal in the transient phase, and thus can be ineffective in curbing the epidemic curve during an outbreak. Moreover, it is important to keep in mind that enacting intervention policies that restrict mobility might be optimal from an healthcare viewpoint, but could potentially lead to tragic economic and social consequences, as it was observed during the COVID-19 pandemic. For this reason, in this paper, we investigate the problem of controlling the transient of an epidemic process, with a trade-off between the beneficial impact to the healthcare system and the economic and social losses.

We consider a time-window $T\in\mathbb{Z}_{>0}$ and we define the following cost function:

[TABLE]

which accounts for the healthcare cost of an epidemic outbreak (proportional to the total number of infections) and the cost associated with implementing the interventions (proportional to the total population that is impacted by the control action). The parameter $\gamma\geq 0$ can be used to weight the importance of the second contribution: larger values of $\gamma$ are associated with giving more importance to the social and economic aspects in the trade-off with healthcare consequences of the outbreak. At this point, we formalize the optimal control problem, as follows.

Problem 1**.**

The optimal control policy is the solution

[TABLE]

where the first constraint is Eq. (1) or Eq. (3), depending on the characteristics of the disease under investigation.

IV QUBO Formulation

We study Problem 1 and, for $T=2$, we derive a QUBO formulation, which is key for solving it in an efficient way using different techniques, including quantum computing [34, 28]. We start observing that, given the initial conditions $\boldsymbol{\bar{x}}=[x_{1}(0),\dots,x_{M}(0)]$ (and also $\boldsymbol{\bar{y}}$ for the SIR model), the evolution of $\boldsymbol{x}$ is fully determined by the two constraints in Eq. (1) and Eq. (5) and thus, ultimately, by the control $\boldsymbol{u}$. Hence, Eq. (6) can be written in terms of $\boldsymbol{u}$, which is the only decision variable of the optimization problem.

In general, being Eq. (1) and Eq. (3) affine in the control input $\boldsymbol{u}$, its recursive use to express $\boldsymbol{x}(t)$ would yield a polynomial of order $t$. Hence, the intuition would suggest that the cost function in Eq. (6) can be written as a $T$th order polynomial. Here, we consider the scenario of $T=2$, for which we can express Problem 1 as a QUBO problem.

Theorem 1**.**

The solution $\boldsymbol{u^{*}}$ of Problem 1 for the SIS model in Eq. (1) with $T=2$ is equal to $\boldsymbol{u^{*}}={\bf 1}-\boldsymbol{z^{*}}$, where $\boldsymbol{z}^{*}$ is the solution of the QUBO problem:

[TABLE]

with

[TABLE]

Proof.

The result is obtained by using Eq. (1) recursively. In fact, for $t=1$, Eq. (1) becomes

[TABLE]

Then, for $t=2$, Eq. (1) becomes

[TABLE]

The expression of $x_{i}(1)$ in Eq. (10) is a first-degree polynomial in $\boldsymbol{u}$. Inserting it into Eq. (11), we obtain

[TABLE]

which is a second-degree polynomial in the control variable vector $\boldsymbol{u}$. Note that in the computations we use the fact that, being $u_{i}\in\{0,1\}$, then $u_{i}^{2}=u_{i}$. Then, inserting Eq. (10) and Eq. (12) into Eq. (6), we obtain a second-degree polynomial expression in $\boldsymbol{u}$ for the cost function. Finally, by performing the change of variables $z_{i}=1-u_{i}$, we obtain the expression in Eq. (9a), modulo some constant terms that have no impact on the result of the $\arg\min$ problem. Hence, the solution $\boldsymbol{z}^{*}$ of Eq. (8) coincides with the solution of Problem 1, with the change of variables $\boldsymbol{u}={\bf 1}-\boldsymbol{z}$. ∎

The same approach can be used for the SIR model in Eq. (3), for which the equivalent QUBO formulation in Eq. (8) can be used to solve Problem 1, with slightly different expressions for the coefficients $P_{i}$ and $Q_{ij}$ due to the change in the dynamics, which are reported in the following, with the proof omitted due to space constraints.

Theorem 2**.**

The solution $\boldsymbol{u^{*}}$ of Problem 1 for the SIR model in Eq. (3) with $T=2$ is equal to $\boldsymbol{u^{*}}={\bf 1}-\boldsymbol{z^{*}}$, where $\boldsymbol{z}^{*}$ is the solution of the QUBO problem in Eq. (8) with

[TABLE]

Remark 3**.**

The technique used in Theorem 1 to cast Problem 1 as a QUBO problem can be directly used only for $T=2$. For $T>2$, the same approach would lead to a higher-degree polynomial. Different techniques may be used in these scenarios to deal with the problem, including exact methods through the addition of auxiliary variables to replace the higher-order terms, and approximated methods based on the recursive use of Taylor expansions.

V Numerical Results

We demonstrate the effectiveness of our proposed approach on a real-world case study, based on an epidemic outbreak spreading in Italy. Specifically, we consider a network where each node represents an Italian administrative division, and links (with their weights) represent individual mobility and travel between these divisions. Using publicly available data from the Italian National Institute of Statistics [36], we are able to build networks with different granularity considering $n=21$ regions or $n=107$ provinces (second-level administrative divisions). The weight associated with links is set to be proportional to the fraction of population that commutes between two administrative divisions from [36].

We consider the SIR epidemic model in Eq. (3), calibrated on the COVID-19 pandemic. Consistently, we set the parameters following [37] and using Remark 1, making sure that the conditions in Proposition 1 are satisfied. We initialize the simulations using the officially reported cases in each administrative division in the desired date, which varies across the simulations to explore the impact of the proposed control techniques in different phases of the outbreak [38].

In particular, for each one of the two levels of granularity (21 regions and 107 provinces), we selected four different initial dates, corresponding to four different phases of the epidemic outbreak: i) March 8, 2020; ii) March 10, 2020; iii) October 29, 2020; and iv) December 20, 2021. The first two dates coincide with early moments of the epidemic outbreak, which start from a moderately low number of cases. The third date captures a scenario of increasing cases due to a novel variant. The fourth date represents a later stage in which COVID-19 has become endemic, but a new wave was approaching due to higher levels of mobility in correspondence to Winter holidays.

In our simulations, we implement the proposed control technique in a rolling horizon fashion, for a total simulation time-window of 30 weeks. Specifically, starting from the initial date of the simulation setting (which is denoted as $t=0$), we solve the control problem at time $t$ using the QUBO formulation in Theorem 2 over a time-window $T=2$, and we implement the optimal control policy $\boldsymbol{u^{*}}$, solution of Problem 1 for a single time step. Then, at time $t+1$, the same procedure is iterated, until the end of the simulation time-window is reached.

The QUBO problem is solved using four different methods, including simulated annealing (SA) and tabu search (TS), both executed using D-Wave Solver, a genetic algorithm (GA) implemented via the pymoo optimization library, and quantum annealing (QA). These methods were selected to provide a representative comparison between classical metaheuristics and emerging quantum technologies [39]. The results of our simulations are reported in Fig. 1, which show the evolution of the epidemic curve over time under each control policy, compared with the uncontrolled baseline. The plots highlight how all methods significantly reduce the number of infected individuals over time.

Table I quantifies these improvements by reporting, for each solver, the relative reduction in the infection peak (denoted by p) and average number of infected/day (denoted by a) with respect to the uncontrolled dynamics; QA achieves performance comparable with those of the state-of-the-art classical methods in most scenarios. Importantly, such results are obtained with a remarkably smaller computational effort, as can be observed from Table II. More details on the computational efficiency of the proposed method and further numerical simulations can be found in [40].

VI Conclusion

In this article, we addressed the problem of optimally controlling the spread of epidemic diseases on networks through mobility bans. Building on classical discrete-time SIS and SIR models, we formulated a control problem to trade-off healthcare impact and socioeconomic costs of a mobility ban policy. We proved that this problem can be cast as a Quadratic Unconstrained Binary Optimization (QUBO) problem, allowing its solution through quantum annealing. We implemented our controller in a rolling-horizon scheme, and we benchmarked the performance of four solvers —including a quantum annealer— on realistic epidemic scenarios calibrated on COVID-19 spreading in Italy, proving the effectiveness and computational efficiency of our approach.

The preliminary results presented in this paper pave the way for several promising avenues of future research. First, our approach was developed under the assumption two-step optimization time-window. In Remark 3, we observed that the use of recursive Taylor expansions is a possible approach to overcome this limitation. Future research should explore this research line. Second, while this paper focuses on control actions in terms of banning mobility between nodes, other types of interventions can be enacted during an epidemic outbreak. Extending our mathematical framework to introduce different control actions and compare their effectiveness is a key direction for future research. Third, the effectiveness of quantum computing in solving this optimization problem suggests that similar framework can be developed also in different contexts of controlling complex network systems.