key: cord-0742273-ol14g74r authors: Zino, Lorenzo; Cao, Ming title: Analysis, Prediction, and Control of Epidemics: A Survey from Scalar to Dynamic Network Models date: 2021-02-27 journal: Ieee Circuits and Systems Magazine DOI: 10.1109/mcas.2021.3118100 sha: 713702c08d65c512fef476e6cfcb5fd817d25da0 doc_id: 742273 cord_uid: ol14g74r During the ongoing COVID-19 pandemic, mathematical models of epidemic spreading have emerged as powerful tools to produce valuable predictions of the evolution of the pandemic, helping public health authorities decide which intervention policies should be implemented. The study of these models -- grounded in the systems theory and often analyzed using control-theoretic tools -- is an extremely important research area for many researchers from different fields, including epidemiology, engineering, physics, mathematics, computer science, sociology, economics, and management. In this survey, we review the history and present the state of the art in the modeling, analysis, and control of epidemic dynamics. We discuss different approaches to epidemic modeling, either deterministic or stochastic, ranging from the first implementations of scalar systems of differential equations to describing the epidemic spreading at the population level, and to more recent models on dynamic networks, which capture the spatial spread and the time-varying nature of human interactions. theory to support the effectiveness of variolation, helping its increasing adoption until the development of the smallpox vaccine in 1796. A latter milestone of mathematical modeling of epidemics is materialized in the studies of the 1849 and 1854 Cholera outbreaks in London by William Farr [2] and John Snow [3] , respectively. However, it is not until the beginning of the 20th century that differential and difference equations started being adopted as tools to model and analyze the spread of epidemic diseases. The very origin of this approach can be found in a paper by William Heaton Hamer [4] , in which the British epidemiologist pioneered the use of a nonlinear formula to model the rate of the contagion process, proportional to the product between the number of susceptible individuals and the number of infectious individuals in the population. Such a modeling approach has been formalized and popularized by William Ogilvy Kermack and Anderson Gray McKendrick. In their seminal paper [5] , and in two subsequent works [6, 7] , the two Scottish researchers laid out the basis for the current mathematical theory of epidemic modeling, namely the susceptible-infected-susceptible (SIS) and the susceptible-infected-removed (SIR) models and the concept of epidemic threshold (and consequently, phase transition). In the SIS model, it is assumed that infected individuals after recovering do not acquire immunity, and thus become again susceptible to the disease. This model became very popular to study recurrent and endemic diseases, including sexually transmitted diseases, such as Gonorrhea [8] . Here, we present the deterministic SIS model in a continuous-time framework, i.e., employing differential equations (as in the original contribution by Kermack and McKendrick [5] ). Note that, discrete-time implementations of the model using difference equations have also been extensively studied. More details on discrete-time models can be found, in e.g., [9] . Formally, a unit mass population is split between two compartments: susceptible and infected. Let s(t) ∈ [0, 1] and x(t) ∈ [0, 1] be the fraction of susceptible and infected individuals in the population at time t ≥ 0, respectively. The health state of the population evolves according to two mechanisms: contagion and recovery, illustrated in Fig. 1a . Susceptible individuals may become infected if they interact with infected individuals, while infected individuals spontaneously recover, becoming again susceptible to the disease. According to the contagion mechanism, the total number of newly infected individuals is a nonlinear expression proportional to the number of susceptible individuals and to the number of infectious individuals in the population (as proposed in [4] ), i.e., equal to where λ > 0 is a positive parameter termed contagion rate. The total number of recoveries in the population is proportional to the number of infected individuals, yielding the term µx(t), where µ > 0 is the recovery rate. Hence, the heath state of the population is described by the two-dimensional state variable [s(t), x(t)] , which evolves according to the following system of ordinary differential equations (ODEs): ṡ(t) = −λs(t)x(t) + µx(t) x(t) = λs(t)x(t) − µx(t) . (2) Note that the two equations are linearly dependent, since the total mass of the population is preserved, that is,ṡ(t) +ẋ(t) = 0. Hence, the system of ODEs in Eq. (2) can be reduced to the single nonlinear ODE:ẋ from which one can immediately observe that the domain [0, 1] is positively invariant. From the analysis of this equation and its equilibrium points, one can conclude that, depending on the model parameters λ and µ, two outcomes may occur: either the origin (disease-free equilibrium) is the unique equilibrium of the system in [0, 1] and it is globally asymptotically stable; or the origin becomes unstable, giving rise to a (almost) globally asymptotically stable endemic equilibrium x * ∈ (0, 1]. These two behaviors are illustrated in Fig. 2 . The phase transition between these two behaviors occur when a certain relation between the model parameters is satisfied, which is called epidemic threshold. In this survey, we opt for expressing such a threshold in terms of the ratio between the contagion rate λ and the recovery rate µ. The following result from [5] fully characterizes the behavior of the SIS model. If the population SIS model in Eq. (3) satisfies λ/µ ≤ 1, then the disease free equilibrium x = 0 is globally asymptotically stable. Otherwise, if λ/µ > 1, the disease-free equilibrium is unstable and Eq. (3) has a (almost) globally asymptotically stable endemic equilibrium, corresponding to x * = 1 − µ/λ. The epidemic threshold is often associated with the basic reproduction number R 0 , that is, the average number of contagions that a single infected person will cause. The concept of basic reproduction number, although already touched upon in the original works on the SIS model [5] , was not formally introduced until the work by George MacDonald in the early 1950s [10] . For the SIS model in Eq. (3), in fact, R 0 = λ/µ. However, in more complex models, e.g., those including additional contagion mechanisms, stochasticity, heterogeneity, and dynamic network structures, the basic reproduction number may not be sufficient to characterize the behavior of the system, which may depend on other factors such as its variability [11] . For this reason, in this survey we prefer to present the results in terms of the epidemic thresholds. In the population SIR model (see Fig. 1b 3) with (a) λ/µ = 0.5 (i.e., below the epidemic threshold) and (b) λ/µ = 2 (i.e., above the epidemic threshold). Below the epidemic threshold, both trajectories converge to the disease free equilibrium; above the epidemic threshold, they both converge to the endemic equilibrium x * = 0.5, denoted by the black dashed line. removed is introduced to account for individuals that recover from the disease and become immune, or die; r(t) ∈ [0, 1] is the fraction of population in the removed state. Hence, the heath state of the system is defined by the three-dimensional state [s(t), x(t) , r(t)] (with s(t) + x(t) + r(t) = 1), which evolves according to the following system of ODEs: where the three equations are linearly dependent, sinceṡ(t) +ẋ(t) +ṙ(t) = 0. Hence, the system of ODEs in Eq. (4) reduces to a planar system. The analysis of such a planar system (see, e.g., [12] ) shows that the ratio λ/µ also plays an important role in the SIR model, modulated by the initial fraction of susceptible individuals s(0). If λs(0)/µ < 1, then the fraction of infected individuals exponentially converges to 0. If λs(0)/µ > 1, then the outbreak monotonically grows until reaching an epidemic peak and, only after the peak is reached, decreasing to 0. Two trajectories of the population SIR model to illustrate these two different behaviors are shown in Fig. 3 . Toward the development of more realistic models, further compartments have been introduced. This allows to account for many realistic phenomena, including latency periods between contagion and infectiveness, multiple stages of the progression of the disease, as well as to model pharmaceutical and nonpharmaceutical interventions, including vaccination and quarantine. More details on this extensive family of mathematical models can be found in [13, 14] . This large variety of models has been widely adopted in the recent years for many theoretical and practical studies of epidemic processes, including seasonal and pandemic influenza [15] , SARS [16] , Ebola [17] , and the ongoing COVID-19 outbreak [18] [19] [20] . Control of epidemics was one of the major motivations for the development of the first mathematical models of epidemics. In fact, Bernoulli proposed his model to support the use of variolation [1] . The idea of exploiting the dynamical system formulation of epidemic Figure 3 : Trajectories of the population SIR model in Eq. (4) with (a) λ/µ = 0.5 (i.e., below the epidemic threshold) and (b) λ/µ = 2 (i.e., above the epidemic threshold). The green, red, and blue curves represent the fraction of susceptible s(t), infected x(t), and removed r(t) individuals, respectively. Below the epidemic threshold, the disease is quickly eradicated and almost all of the population remains susceptible; above the epidemic threshold, one observes an initial growth of the epidemic prevalence and the majority of the population is eventually removed. models to derive and assess control techniques can be traced back to the work on Malaria by Robert Ross [21] , in which the British medical doctor claimed that the reduction of the population of mosquitoes below a certain threshold would suffice to eradicate the disease. The development of compartmental models has laid the foundation of a rigorous mathematical framework, grounded in the theory of dynamical systems, to study the spread of epidemic processes and assess the effectiveness of different control strategies implemented to mitigate and combat them, by employing control-theoretic techniques. However, this class of deterministic compartmental models suffers from some inherent limitations, which may hinder their applicability to real-world scenarios, thus calling for several extensions to be discussed in the rest of this survey. The first limitation is concerned with the use of differential equations to model the dynamical process. This usage is under the assumption of having a large population that can be approximated by a continuous distribution; it also assumes that the presence of noise and uncertainties can be omitted. However, epidemics often start from small local outbreaks, where finite population effects and stochasticity cannot be neglected. To account for these important factors, compartmental models of epidemics in populations have been extended to stochastic frameworks, by using the theory of Markov processes and branching processes. Exhaustive surveys of these models can be found in the following books [13, 22] . The second limitation lies in the fact that this class of deterministic compartmental models provides a description of the epidemic process only at the population level, that is, in terms of the overall fractions of susceptible and infected individuals in the system. In real-world epidemic outbreaks, the pattern of interactions between the members of a population -and, consequently, the contagion process -has a complex structure, which may be influenced by the spatial location of the individuals of a population and the pattern of their interactions, calling for a finer representation of the population structure. In the last few decades, network theory has emerged as a powerful tool to capture such a complexity and formalize it in mathematically tractable models [12, [23] [24] [25] . In the rest of this survey, we will focus on network epidemic models, mostly discussing the SIS model, which is among the simplest and most studied models. In Section 2, we will review some classical results for its deterministic and stochastic implementations on static networks. In Section 3, we will move to the most recent developments and the state of the art for the research on epidemics on dynamic networks. Section 4 is devoted to the discussion of control of deterministic and stochastic network epidemic models, with specific attention on the control of dynamic networks. In Section 5, we present a focused discussion on the current challenge of the ongoing COVID-19 pandemic and on how the mathematical modeling of epidemics can provide powerful tools toward fighting against and mitigating its spreading. Finally, in Section 6, we briefly summarize the content of this survey and outline some promising challenges for the future research. Here, we present the network SIS model and report the main results for its deterministic or stochastic implementations. Then, we briefly mention some important extensions of these results to scenarios that incorporate further features of real-world systems. Towards this end, we first gather some definitions and basic notions on network and graph theories, which will be used throughout the survey. A static network is represented by means of a time-invariant graph with n nodes, denoted by a set of positive integer indices V = {1, . . . , n}. Nodes are connected through a set of directed links E ⊆ V × V, such that (i, j) ∈ E if and only if node i is connected to node j. The intensity of such connections is measured by a nonnegative quantity. Specifically, for any pair of nodes i, j ∈ V, we define the connection strength a ij ≥ 0 that quantifies how strongly node i is connected to node j, and so, a ij > 0 ⇐⇒ (i, j) ∈ E. The connection strengths are collected in a (weighted) adjacency matrix A ∈ R n×n ≥0 . The triple G = (V, E, A) defines the graph. An example is illustrated in Fig. 4 . A network is undirected if the corresponding adjacency matrix is symmetric, that is, A = A ; otherwise it is said to be directed. A network is connected (strongly connected, for directed networks) if its adjacency matrix A is irreducible, that is, if for any pair of nodes i and j, there exists a sequence of nodes v 1 = i, v 2 , . . . , v k = j such that (v , v +1 ) ∈ E, for = 1, . . . , k − 1. A network is unweighted if the adjacency matrix of the corresponding graph A has binary entries, that is, all nonzero entries (corresponding to links) are equal to 1. Given a node i ∈ V, we denote by k i = j a ij its (weighted) degree. Note that if the network is unweighted, then the degree k i is equal to the number of nodes that node i is connected to, which are called neighbors of i. A dynamic network is represented by means of a time-varying graph G(t) = (V, E(t), A(t)), where the time t can be a discrete or a continuous index. The n nodes in the node set V = {1, . . . , n} are time-invariant and connected through a time-varying set of links E(t). The matrix A(t) ∈ R n×n ≥0 is the time-varying (weighted) adjacency matrix and measures the strengths of the connections between nodes at time t. From the first implementation of a deterministic SIS model on a (static) network, proposed by Ana Lajmanovich and James A. Yorke in 1976 to study the spread of gonorrhea in a population [26] , network epidemic models have become an established and successful paradigm to study the spread of epidemic diseases in complex populations [12, [23] [24] [25] . In the network SIS model, each node represents a group of individuals and is characterized by its health state, that is, (s i (t), x i (t)). The health state of node i ∈ V represents the fraction of individuals per health states in the ith group [26] . Hence, s i (t) ∈ [0, 1], x i (t) ∈ [0, 1], and, similar to the population SIS model described in the previous section, it holds true that x i (t) + s i (t) = 1. Hence, the health state of each node can be fully determined by the unique variable x i (t) and the entire state of the population is thus represented by the n-dimensional vector x(t) ∈ [0, 1] n . Similar to the scalar model, the contagion rate in node i has a nonlinear expression, equal to the product of the infection rate λ i , the fraction of susceptible individuals in the ith group s i (t) and the strength of the interactions with other infected individuals in the population m i (t), which depends on the network structure and is equal to as illustrated in the simple example in Fig. 5 . For the sake of simplicity, in this section, we will present our main results under the assumption that the graph G = (V, E, A) is strongly connected. Results for directed networks that are not strongly connected can be found in [27, 28] . Hence, the network SIS model is governed by a set of n ODEs, one for each node of the network, where the equation that determines the evolution of x i (t), for a generic node i ∈ V, is the following:ẋ where λ i and µ i are the infection and recovery rate of node i, respectively. While in most of the literature it is assumed that the epidemic parameters are homogeneous in the population, that is, λ i = λ and µ i = µ, for all i ∈ V, the general formulation in Eq. (6) with heterogeneous parameter has been adopted to model and study different scenarios, for instance, to model diseases that affect different age groups differently, which is key to implement and study targeted interventions and vaccination policies. Note that, after introducing the n-dimensional vectors λ and µ to gather all the λ i and µ i , i ∈ V, respectively, the dynamics can be written in a compact, vector form aṡ where Λ = diag (λ), M = diag (µ), and 1 is an n-dimensional all-1 vector from which is straightforward to observe that the domain [0, 1] n is positively invariant. The following result, initially presented in [26] , extends the results for the population SIS model in Theorem 1 to the network SIS model. Different techniques have been used to prove this extension, including Lyapunov arguments [29] and positive systems theory [28] . Theorem 2. Consider the homogeneous network SIS model in Eq. (6) with λ i = λ and µ i = µ, for all i ∈ V, on a strongly connected network. If where ρ(A) is the spectral radius of the (weighted) adjacency matrix A, then the disease free equilibrium x = 0 is globally asymptotically stable. Otherwise, if λ/µ > 1/ρ(A), the diseasefree equilibrium is unstable and Eq. (7) has a unique (almost) globally asymptotically stable endemic equilibrium x * . Different from the population SIS model, where the endemic equilibrium x * has a closedform expression depending on the model parameters, for the network SIS model such an explicit formula cannot be derived in general. However, iterative algorithms to compute such an equilibrium have been proposed in the literature. See, e.g., [12, 30] . The convergence result in Theorem 2 can be easily extended to heterogeneous SIS models, where nodes have different contagion rates λ i and/or recovery rates µ i , as shown in [26, 28, 31] . In the heterogeneous case, the behavior of the system is determined by the spectral radius of the matrix ΛA − M . Specifically, if the disease free equilibrium x = 0 is globally asymptotically stable. Otherwise, the SIS dynamics in Eq. (7) converge to the unique (almost) globally asymptotically stable endemic equilibrium x * . Note that this expression reduces to Eq. (8) in the homogeneous scenario. A technique to approximate the solution of Eq. (7) about the epidemic thresholds has been proposed in [32] . Similarly, the SIR model and more complex compartmental models have been embedded and studied on network structures. For more details on these implementation and their analysis, we refer to the following review papers [12, 24, 25 ]. In their stochastic implementation, network epidemic models are defined as follows. Each node of the network represents an individual and is characterized by the health state X i (t), which coincides with one of the compartments. In the stochastic network SIS model, we have The nodes' health states are gathered into an n-dimensional vector X(t) ∈ {0, 1} n , which represents the health state of the population. Most of the literature on stochastic epidemic models relies on the assumption that the evolution of the epidemic process X(t) can be represented by a Markov jump process, where the state transitions are triggered by Poisson clocks. Such an assumption, although simplistic, allows to use the rich theory on Markov processes [33] to perform rigorous analyses of the model, determining its asymptotic and transient behavior, as detailed in the following. We assume that the vector X(t) evolves according to a continuous-time Markov process, governed by the contagion and the recovery mechanisms, which act on the health state of each node. The former regulates a node's state transitions from the susceptible state to the infected and is modeled by a Poisson clock with the rate proportional to the sum of the connection strengths of the infected individuals that node i is in contact with and to the contagion rate λ i , that is, The recovery mechanism determines the state transitions from the infected state to the susceptible and depends only on the recovery rate of the individual, that is, it is a Poisson clock with the rate equal to λ R i X(t) = µ i . These two mechanisms yield the following transition probabilities for all i ∈ V, which unequivocally determine the dynamics of the Markov process, as illustrated in Fig. 6 . At this stage, one may observe the presence of strong similarities between the expressions of the transition probabilities of the Markov process and the differential equations that govern the deterministic network SIS model. Indeed, these two models have strong connections. In fact, recently, a different formalization of the deterministic network model has been proposed. In this formalization, each node of the network is a single individual of the population (similar to the stochastic framework) and the corresponding health state x i (t) denotes the probability that node i is infected at time t, that is, [30, 34] . Under the so-called N-intertwined mean-field approach, the following approximation is made: , showing that the deterministic system of ODEs in Eq. (7) approximates the evolution of the expected value of the stochastic process X(t). This approach, even though not exact (since it is based on an approximation that typically does not hold true), is often used to approximate and study the evolution of more complex stochastic network epidemic models. From Eq. (12), we observe that the disease free equilibrium x = 0 is the unique absorbing state of the Markov process and is globally reachable. Hence, different from its deterministic counterpart, in the stochastic SIS model, the disease is always eradicated with probability 1 in finite time [33] . However, in the following one will see that an epidemic threshold is still present in terms of the duration of the transient evolution of the epidemic outbreak, before reaching the disease-free equilibrium. Formally, it has been observed that a phase transition with respect to the eradication time can be established, where the latter is defined by Specifically, a set of results that characterize the expected value of the eradication time E[T ] depending on the model parameters and on the network structure has been established [35] [36] [37] . The following theorem gathers some key results from the cited literature. Theorem 3. Consider the homogeneous stochastic network SIS model in Eq. (12) with λ i = λ and µ i = µ, for all i ∈ V, on a strongly connected network. If where ρ(A) is the spectral radius of the (weighted) adjacency matrix A, then where K 1 , K 2 > 0 are two nonnegative constants that depend on the model parameters and on the network structures. The proofs of these results are quite technical and based on the theory of Markov processes [33] . Briefly, the key idea is that, in the fast extinction regime, the probability that the number of infected individuals in the population increases as determined by Eq. (12) is always less than the probability that it decreases, yielding a drift in the direction of the disease-free equilibrium. Above the epidemic threshold, instead, the inequality is reversed when the process is close to the disease-free equilibrium, implying that the infections tend to rise and large stochastic fluctuations are needed to reach the disease-free equilibrium. Similar techniques have been used to study other stochastic epidemic models on networks, including the SIR model. For more details, we refer to these review papers [25, 38] . The results summarized in Theorem 3 show the existence of a sharp phase transition between a regime where the epidemic is on average quickly eradicated, and a regime where the disease lasts on average for an exponentially long time. However, these results on the average eradication time may fail in characterizing the actual behavior of a single instance of the epidemic process. On the one hand, quick eradication can be guaranteed by directly applying the Markov inequality to Eq. (15), yielding T ≤ (ln n) α with probability converging to 1 as n → ∞, for any α > 1. On the other hand, an accurate analysis of the eradication time has been performed in [39] , where it has been established that, if λ µ > σ ρ(A) , where σ ≥ 1 is a constant that depends on the network structure (more specifically, on the isoperimetric constants of the network), then the eradication time is exponentially large in the number of nodes with probability converging to 1 as n → ∞. Note that, if σ > 1, there is a gap between the two regimes. In some specific cases (e.g., for complete graphs and Erdős-Rényi random graphs), σ may be equal to 1, yielding a sharp phase transition not only in the expected duration of the epidemic disease, but also in its actual duration, with high probability. Recent efforts have been made to relax the Markovianity assumption, assuming different forms for the process X(t). In particular, an SIS model, in which the statistical distribution of the contagion time and/or of the recovery time differ from the exponential distribution associated with Markov processes, was proposed in [40] . Therein, a mean-field approach is used to determine conditions for fast eradication of the disease. Further extensions of this approach can be found in [41, 42] . Without relying on any mean-field approximations, in [43] , a lower bound on the decay rate to the infection-free equilibrium is rigorously computed. Table 1 : Notation for network epidemic models Symbol Description G = (V, E, A) network (node set, node set, adjacency matrix) ρ(A) spectral radius of the adjacency matrix A n number of nodes X i (t) health state of node i (stochastic models) λ i infection rate of node i µ i recovery rate of node i x * endemic equilibrium T eradication time (stochastic models) Equivalences and differences between Markovian and non-Markovian epidemic models have been extensively discussed in [44, 45] . All these works suggest that the non-Markovianity of the mechanisms that govern the epidemic process may have a large impact on the spread of a disease and outline an important avenue of future research in the field of stochastic epidemic models toward shedding lights on how the distributions of infections and recovery times shape the spreading process. This survey focuses mostly on continuous-time epidemic models. However, it is important to mention that the continuous-time formulations of deterministic and stochastic models presented in this survey (both the population and the network models) naturally have discretetime counterparts, where differential equations and Markov processes are replaced by difference equations and Markov chains, respectively. Here, we report the equations for the discrete-time deterministic network SIS model, which will be used in some of the models of epidemics on dynamic networks presented in Section 3. For each node i ∈ V, the health state is updated as follows: where m i (t) is defined in Eq. (5); here, λ i and µ i have to be interpreted as the per-contact infection probability and the per-time-unit recovery probability, respectively. Most of the results discussed in this section concerning the existence of a phase transition between a fast extinction regime and a regime where the epidemics becomes endemic and its dependence on the model parameters and on the network structure can be extended with some minor adjustments to the discrete-time counterparts of the models. For more details on the analysis of discrete-time epidemic models and their main results for deterministic models, we refer to [46] [47] [48] [49] [50] , for which a recent review paper by Parè et al. covers most of the results [51] ; for stochastic models, we refer to [52, 53] . Detailed discussions on the main differences between continuous-time epidemic models and their discrete-time counterparts can be found in [54, 55] . A few years ago, L. Pellis et al. in a perspective paper outlined eight important challenges for network epidemic models [56] . Besides other -more practical -directions, calling for the integration of network computational modeling and epidemiological relevant data, two key challenges were identified, which are of great interest for the engineering community. The first challenge concerns the study of epidemic models on dynamic network structures, leading to the study of a nonlinear time-varying dynamical system. The second focuses on understanding how the network structure (static or dynamic) can be exploited to effectively design intervention policies to stop or mitigate the disease spreading; control-theoretic tools are key to address this second challenge. In the rest of this survey, we focus on these two research directions, presenting the state of the art in terms of key progresses of the last few years and most promising lines of current research. The extension of the classic compartmental models to static networks and the corresponding rigorous analysis have allowed the scientific community to understand how the architecture of human social interactions affects the spread of epidemic diseases in interconnected populations. However, in most of the real-world epidemic outbreaks, the underlying network of social interactions is not static, but dynamically changes, co-evolving with and influenced by the spread of the disease [57] [58] [59] [60] . Several endogenous and exogenous reasons may be adduced to explain and motivate the dynamic evolution of the network structure. First, dynamical changes of the individuals' patterns of interactions may be directly or indirectly caused by seasonal factors, such as school holidays and weather conditions, which may favor or hinder gatherings and social events [58, 61] . Second, social interactions are often characterized by an intermittent behavior, whereby individuals' propensity to generate connections is subject to burstiness, yielding clusters of connections separated by latency periods [62, 63] . Third, the infection events themselves may affect the network structure. In fact, not only infected individuals may reduce their interactions as a consequence of the illness, but also susceptible individualsdriven by awareness and risk perception cognitive mechanisms -may dynamically modify their behavior, reducing or rewiring their interactions to avoid contagion [57, 64, 65] . Finally, nonpharmaceutical intervention policies often entail a dynamical modification of the pattern of human interactions, which may be dynamically reduced through the implementation of lockdown or social distancing policies, or reshaped by travel bans and mobility limitation, as observed during the ongoing COVID-19 pandemic [66] . All these evidences of the dynamic nature of social interactions call for the extension of the epidemic models presented in the previous sections to dynamic networks. Here, we review and discuss some successful modeling paradigms that have been recently developed to capture this dynamic nature of the network of social interactions within analytically tractable mathematical models of epidemic diseases. Through the analysis of these models, we shed lights on how the dynamical nature of human interactions plays a key role in shaping the evolution of the epidemic outbreak. As we shall illustrate in the next section, the insight gained into the epidemic process through these analyses has enabled researchers to propose and assess valuable intervention strategies to control and mitigate the epidemic spreading, taking into account and leveraging the dynamical properties of social interactions. The first class of approaches to deal with dynamic networks has extensively relied on timescale separation techniques. These techniques are based on the assumption that the epidemic process and the network dynamics evolve at different paces, as illustrated in Fig. 7 . If the epidemic process evolves much faster than the network of interactions, the system is in the so called quenched regime. In this regime, static networks are accurate proxies of slowly switching topologies, and the corresponding results presented and discussed in the previous section are thus used to study the evolution of the epidemic outbreak. On the other extreme, we encounter the annealed regime, in which the evolution of the network is assumed to be much faster than the disease spreading process [67] . In this regime, if the following limit existsĀ then we can define the average graph G = (V,Ē,Ā), with (i, j) ∈Ē ⇐⇒ā ij > 0, and the dynamic network can be effectively represented and studied by means of its average graph. Also in the annealed regime, results on static networks are applied to the average graph to study the behavior of the dynamical system. In the physics community, epidemics on annealed networks have been extensively studied, aiming at computing -or approximating -the epidemic threshold of the model on a network with known degree distribution, without going through the explicit computation of the spectral radius of the average adjacency matrix. Specifically, by relying on a degree-based mean-field approach, the epidemic threshold for the SIS model on unweighted networks has been computed in [68] as a function of the degree distribution of the network. Note that such an approximation is exact in the limit of large-scale systems n → ∞. Specifically, let denote the mean and the second moment of the degree distribution of the average network G = (V,Ē,Ā), respectively. Then, the following expression for the epidemic threshold can be obtained for uncorrelated annealed networks (that is, ifā ij ∝k ikj ). Theorem 4. Consider an SIS model on a dynamic network in the uncorrelated annealed regime with λ i = λ and µ i = µ, for all i ∈ V. Let us define the following epidemic threshold: Then, in the limit n → ∞, if λ/µ < σ, the disease-free equilibrium is asymptotically stable; otherwise, if λ/µ > σ, the disease-free equilibrium is unstable. From Theorem 4, we establish that, below the epidemic threshold in Eq. (20) , the epidemic is in the fast extinction regime while, above the epidemic threshold, the epidemic becomes endemic. An important implication of Eq. (20) applies to scale-free networks, whose degrees follow a power-law distribution. In fact, in many real-world applications, complex networked systems are modeled by scale-free networks with the power-law exponent between 2 and 3 [69] . In these scenarios, the expression of σ in Eq. (20) vanishes as n → ∞, implying that for any nonzero contagion rate λ > 0, epidemics always spread on large scale-free networks. More details on this approach and on further extensions of these works to more general networks (including correlated networks) and to more complex epidemic models (including the SIR model) can be found in [24] . The quenched and annealed regimes discussed in the above rely on the assumption that the epidemic process and the network evolve at different paces and, thus, on different timescales. However, the arguments raised at the beginning of this section to motivate the need for dynamic networks provide evidence that such a time-scale separation is often restrictive and unrealistic, since the contagion process and the network evolution are often intertwined and thus often evolve at comparable time-scales. In the last few years, several efforts have been made to overcome the limitation of time-scale separation and propose a theory for epidemics on dynamic networks which allow to model and study the coevolution of the two dynamical processes. In the rest of this section, we present and discuss some of these fascinating paradigms. The use of temporal-switching networks to study epidemics in time-varying systems has been initially proposed in [70] . In its original incarnation, a switching network is generated by repeating a deterministic sequence of T static networks, characterized by the adjacency matrices A 1 , . . . , A T , so that A(t) = A t mod T . Therein, the discrete-time deterministic SIS model has been studied, extending the results found for static networks. Specifically (in the homogeneous scenario), if one defines then, the behavior of the system is determined by the spectral radius ρ(P ). If ρ(P ) < 1, then the disease-free equilibrium is asymptotically stable and the disease is quickly eradicated. Otherwise, the disease-free equilibrium becomes unstable and the disease becomes endemic. Such a framework has been extended to the discrete-time stochastic SIS model in [71] , showing the same epidemic threshold by mapping the time-varying system onto a multi-layer network structure. In [72] , sufficient conditions for stability have been established for a more general scenario, where the networks switch arbitrarily among a set of topologies, possibly according to stochastic mechanisms, such as Markov switching rules where, given a Markov process σ(t) with the state space A 1 , . . . , A T , we set A(t) = A σ(t) . This sufficient condition is expressed in terms of the maximum possible norm of products of matrices P t in the set For a review of the most recent developments of this theory, including the computation of a unified formula for the epidemic threshold, we refer to the following paper by Zhang et al. [73] . The use of temporal-switching networks to study epidemics on dynamic networks has been recently extended to continuous-time processes. Specifically, in [74] , the theory of positive linear switched systems is leveraged to derive conditions for global asymptotic stability of the disease-free equilibrium. Such conditions are obtained by combining the stability analysis of the linearized switched system an appropriate notion of irreducibility for the linearization. Specific results are obtained if the topology evolves stocastically according to a Markov switching rule. This approach is followed in [75] to design Markov switching laws to enforce quick eradication of the disease via geometric programming. In [76] , a continuous-time network SIS model is studied in a scenario where the network topology switches deterministically, at discrete-time steps, following a sequence of adjacency matrices. Therein, the epidemic threshold is computed in terms of the Lie commutator bracket of the adjacency matrices, showing that adjacency matrices that are non-commuting yield a lower epidemic threshold, favoring the epidemic spreading. In [77] , the scenario of continuous-time switching networks is considered. The epidemic threshold is explicitly computed in the case when the adjacency matrix A(t) commutes with the aggregated adjacency matrix up to that time t 0 A(s)ds. Under this condition, the order of the switching matrix has no effect on the dynamics, which is fully determined by the average adjacency matrixĀ, defined in Eq. (18) . If λ/µ < 1/ρ(Ā), then the disease-free equilibrium is asymptotically stable; otherwise, it is unstable and the disease becomes endemic. In [78] the analysis of the continuous-time deterministic heterogeneous SIS model from [31] is extended to slowly changing time-varying systems by leveraging Lyapunov arguments. In that work, the effect of stochastic perturbations of the model is also studied. A major limitation of the theoretical analysis of temporal-switching networks is that it is often assumed that the network switches between different adjacency matrices determined a-priori and that the switches typically take place at deterministic (often equally-spaced) time-instants. In the following, we will present two alternative paradigms that do not rely on this assumption. In all the network models presented so far, interactions were determined by some predetermined connectivity patterns, represented by an adjacency matrix, or a sequence of adjacency matrices. Activity-driven networks (ADNs), originally proposed by N. Perra et al. in [79] suggested to perform a paradigm shift. In ADNs, interactions are not seen as a consequence of a network structure which identifies pre-determined dyadic relationships between pair of nodes, but are rather generated by individuals' properties, according to a stochastic process. In their original incarnation, each individual i is characterized by a unique parameter a i ∈ [0, 1], termed activity, which represents an individual's propensity to generate interactions. Then, starting from t = 0, the dynamic network is generated according to the following algorithm: The formation process of an ADN is depicted in Fig. 8 . Note that, at each discrete-time, the network generated by the ADN process is always undirected, but not necessarily connected. The main strength of this paradigm lies in its simplicity: the temporal nature of the system is encapsulated in a single n-dimensional vector. Such a simple formulation has enabled to perform rigorous analytical studies of the properties of the network generated and of dynamical processes evolving on it. Specifically, in [79] , the following result has been established. Theorem 5. Consider an SIS model on an ADN with λ i = λ and µ i = µ, for all i ∈ V. Let us define the following epidemic threshold: For λ/µ < σ the epidemic is quickly eradicated with probability converging to 1 as the network size n → ∞, while for λ/µ > σ the epidemic becomes endemic with probability converging to 1 as the network size n → ∞. The original formulation of ADNs was proposed in a discrete-time framework, so the discrete-time counterpart of the SIS model in Eq. (17) was studied. A continuous-time formulation of ADNs has been proposed in [80] , where the synchronous activations ruled by the activity-based mechanism are substituted by an asynchronous mechanism, where each node is activated according to a Poisson process with the rate equal to its activity, yielding thus a Markov process. Therein, the continuous-time SIS model is analyzed, giving rise to the same threshold identified in Theorem 5. One of the major advantages of the ADN formulation is its simplicity that, besides enabling rigorous analytical studies, allows to expand the formalism in several directions. In fact, in the last few years, several extensions and generalizations of the ADN paradigm have been proposed. These extensions allow to include several features of real-world complex systems in the model. The analytical tractability of ADN-based models has enabled the rigorous computation of the epidemic threshold for these models and the characterization of their behavior, similar to Theorem 5, shedding light on the role of these real-world phenomena on the spreading of epidemics. These extensions include the presence of preferential connectivity patterns [81] [82] [83] , community structures [84, 85] , heterogeneous propensity to receive connections [86] , memory and burstiness in the link formation process [87] [88] [89] , and high-order relations [90] . A detailed review of the results for classical activity-driven networks and for the explicit results of the epidemic thresholds for these recent extensions can be found in [91] . A different approach, which to a certain extent combines the presence of a connectivity pattern, determined a-priori, and the stochasticity of its evolution, are edge-Markovian dynamic graphs. This paradigm has been proposed in [92] to model stochastic evolution of dynamic networks. In edge-Markovian dynamic graphs, each potential link (edge) of the graph (i.e., each pair (i, j) ∈ V ×V, with i = j) is associated with a two-state Markov chain (independent of the other links), where the two states represent the existence and nonexistence of the link, respectively. Two probabilities p, q ∈ [0, 1] are defined so that, at each discrete time-step, the chain switches from nonexistence to existence with probability p, while the opposite transition happens with probability q, as illustrated in Fig. 9 . In plain words, the network is initialized as a given (static) network. Then, each link that exists at time t disappears at the following time-step with probability q (independent of the others), while nonexisting links at time t appear with probability p (independent of the others). A continuous-time formulation of the model can be obtained by substituting the Markov chains with continuous-time Markov processes [93] . Epidemic processes on edge-Markovian dynamic graphs have been proposed and studied in [93, 94] . Specifically, in [94] , the authors derived the value of the epidemic threshold in terms of the basic reproduction number, which has a complex expression. Edge-Markovian dynamic graphs are amenable of several analytically tractable extensions to overcome the limitation posed by its original formulation. For instance, the Markov chain (or process) underlying the network evolution determines the inter-event time distribution for the appearance and disappearance of the links (geometric in the discrete-time model, exponential in the continuous-time model), which is not consistent with the presence of burstiness and temporal clustering discussed in the above [62, 63] . In [95] , this paradigm has been extended to account for general inter-event time distributions, establishing a sufficient condition for (almost sure) exponential stability of the disease-free equilibrium. For engineering researchers who are familiar with the monitoring, intervention or control of complex systems, it is of great interest to know what research works have been carried out to study how to influence, mitigate, and even stop the epidemic processes, especially those based on the models we have presented in the previous sections. Although much fewer control results have been produced compared to the epidemic modeling activities, it is beneficial to give an overview of the existing results, which will serve to inspire researchers, with or without control theory background, to work in this critically important research area. In what follows, we categorize the corresponding results into control of deterministic and stochastic epidemics respectively, and to underscore the importance, devote a separate subsection for the discussion of the distinct features when the underlying networks are dynamic. Note that an earlier survey [25] has summarized some main results on control epidemics on networks by then, and thus we give special attention, wherever appropriate, to those that appeared in The simplest idea of controlling epidemics processes on networks comes from the intuition that removing infected or high-risk individuals and the links associated with them will slow down the propagation, which in practice translates to quarantine and vaccination policies. Following the discussion in the previous sections, this intuition implies that one can lower the epidemic threshold, e.g., by reducing the spectral radius ρ(A) of the adjacency matrix A. Another intuitive idea is to optimize the distribution of antidote, which in practice translates into modifying the entries of matrix M in Eq. (7). These key control actions are illustrated in Fig. 10 . Naturally, there is always a cost associated with the control actions and consequently optimal or sub-optimal control objectives can be formulated. However, to get a flavor on why such intuitive control problems under cost constraints are difficult to solve, we look at the following direct formulation of the control problem. Although this control problem is straightforward to formulate, it is very difficult to be solved analytically. In fact, it was shown in [96] that Problem 1 is NP-hard. A similar spectral minimization control problem through removing nodes is discussed in [97] . In fact, the difficulty in solving such control problems is rooted in the fact that the formulated optimal control problems are variations of the constrained combinatorial optimization problems, which are in general hard problems. For this reason, various heuristics have been proposed to solve the control problems, and a lot of them have smartly taken advantage of the network structures, e.g., the degree distribution of the nodes, centrality indices, and connectivity patterns [98, 99] , or solving a non-convex quadratically constrained quadratic program [100] . Another intuitive approach lies in tuning the values of the parameters by increasing the recovery rates µ i to minimize the expression in Eq. (9) , which characterizes the epidemic threshold. Such an approach can be interpreted as a resource allocation problem, where limited amounts of antidote can be distributed to the population. Such an intuition is formalized in the following optimization problem. Problem 2. Given a network and a fixed budget B > 0 where f : R n → R + is a cost function associated with the cost of increasing the recovery rate and µ i (µ i ) is the minimum (maximum) admissible recovery rate for node i ∈ V. Even though the objective function is a spectral radius, which in general is non-convex, under reasonable assumptions on the cost function, tools from geometric programming and convex optimization can be leveraged to tackle the problem. Solutions have been proposed in a centralized fashion [101] [102] [103] [104] , and through distributed approaches [105] [106] [107] [108] [109] . Some of these works deal with a more general problem, in which, besides increasing the recovery rate, the controller can also reduce the infection rates λ i modeling, for instance, the distribution of personal protective equipment. A resource allocation problem similar to Problem 2 is studied for an extension of the SIS model, in which complication phenomena of the illness are considered [110, 111] . When considering more complex dynamics than the standard SIS epidemic model, ideas from optimal control have already been applied a decade ago [112] , using a linear-quadratic regulator, and, more recently, in [113] , leveraging the Pontryagin's Maximum Principle. More recently, researchers have identified impossibility results which reveal the possible limitation of feedback control. In [114, 115] , the authors have proved that, utilizing the recovery rate µ i as control input, a large class of distributed controllers cannot guarantee convergence to the disease-free equilibrium. In [116] , a similar result has been proved for more complex dynamics that involve two concurrent epidemic processes. The limitations may become even more profound when examining the effect of optimal control in real-world disease management [117] . Because of the stochastic nature of the models, the related control results are centered around evaluating the control performance in terms of bounding as tightly as possible the epidemic thresholds on different classes of networks. In [118, 119] , the studied control problem is how to distribute a fixed amount of antidote to nodes in the given network that may have special topological features, e.g., scale-free networks. In [118] , two different methods are compared, one based on contact tracing, which augments the recovery rates of all neighbors of an infected node, and the other based on degree-centrality, which augments the recovery rates of all nodes, proportional to their degrees. Surprisingly, it is found that contact tracing may only succeed when the number of infected individuals is small (e.g., in early stages of the epidemic outbreak), since it requires a total amount of antidote B = i µ i that grows super-linearly in the number of contacts; otherwise the degree-centrality based approach outperforms contact tracing, as stated in the following result. Theorem 6. Consider the stochastic network SIS model in Eq. (12) on a generic network. If µ i ≥ λ i k i , then the expected eradication time verifies E[T ] ≤ K ln n, for some constant K > 0. However, from this result, we note that the needed amount of antidote B scales linearly with the sum of the weights of the links in the network, hindering its practical implementation even in sparse networks, where such a sum grows linearly in the number of nodes. In a similar setup, in [119, 120] , another control method is proposed, in which the antidote is dynamically allocated to the nodes. The allocation method in [119] requires that all the antidote B is concentrated at a single node at each time; using martingale theory, the following result is established. Theorem 7. Consider the stochastic network SIS model in Eq. (12) on a generic network with the control policy proposed in [119] . Let us define the maximum degree and the cutwidth of the network as ∆ = max i∈V k i and W = min Briefly, the proposed control technique guarantees fast eradication with a sub-linear amount of antidote, depending on the network topology. In [120] , a fundamental limitation is further established for any dynamic allocation by showing how the network topology influences the possibility of eradication the epidemics under the given budget of the antidote to be allocated. A similar control policy based on dynamical resource allocation has been proposed in [121] . Very recently, the powerful tool of model predictive control (MPC), already used for deterministic epidemic models [122] , has been used to deal with the control of stochastic epidemic processes without relying on mean-field approximation [123] . The use of MPC has also allowed to deal with more complex epidemic models. For instance, in [124] , the authors deal with a model with a pre-symptomatic phase, by utilizing MPC with a robust moment closure technique. Optimal control has also been considered for stochastic SIS model with different assumptions on the information available [125] . When the networks are dynamic, the optimization problem in control may have to face time-varying constraints. Fortunately, some of the geometric program techniques still apply, although the complexity in seeking the solution increases [126] . Optimal control theory can be applied to time-varying systems; however, it is well known that the corresponding stability conditions might be more conservative and more difficult to check [127] . In [75] , the authors deal with a resource allocation problem similar to Problem 2, on temporalswitching networks, under the assumption that the switching is determined by a Markov process. Therein, the problem is solved via geometric programming, finding a solution with a computational complexity that grows super-linearly -but polynomially -in the network size. For some specific cases, control strategies for eradicating the outbreak have been successfully proposed. For instance, in [128] , a distributed control scheme has been designed for a deterministic SIS model on periodic time-varying networks. In this scheme, that recalls the antidote distribution proposed in [118] and summarized in Theorem 6, each node dynamically sets its recovery rate proportional to the sum of the weight of the links to its neighbors at that time, that is, Under this scheme, global asymptotic stability of the disease-free equilibrium is guaranteed, under the assumption that all the instances of the time-varying network are strongly connected [128] . An interesting research line is to incorporate the human behavioral mechanisms that lead to variations of the dynamic network structures [65] . Understanding human behavior will be critical to draft and implement vaccination policies [129] , which affects the dynamics of the networks and at the same time is deeply affected by the dynamics of the networks [129] . Further study may dive in how isolation of infected individuals may adaptively change and reshape the network dynamics and how this can be leveraged to devise effective intervention policies [130, 131] . In particular, in [130] , the authors explicitly compute the following epidemic threshold for a network SIS model on activity-driven networks, depending on the possibility of isolating infectious individuals by decreasing their social activity to a factor p ∈ [0, 1], where p = 1 means that interventions are enacted and p = 0 models a complete isolation of infected individuals: By comparing Eq. (26) with Eq. (22), one can assess the effect of isolation policies on increasing the epidemic threshold. More recently, awareness-based control strategies, which were developed a decade ago [132] , have been extended to study temporal networks such as activity-driven networks [133] [134] [135] . The proposed formalism allows to study scenarios in which aware individuals reduce the risk during their physical interactions [134] or they dynamically rewire themselves to avoid interactions with infected individuals [133] , possibly in combination with other control policies, such as isolation of infected individuals [133, 135] . In the model proposed by Yang et al. in [134] , awareness is modeled as a process that co-evolves with the epidemic spreading on a two-layer network. The epidemic process spreads on a physical contact layer, while awareness spreads on a virtual communication layer. Aware individuals reduce their infection probability, as a consequence of the adoption of self-protective practices. The epidemic threshold is then computed as a function of the individuals' activities on the two layers and the coupling between them. In the activity-driven adaptive-SIS model proposed in [133] , two modifications are made to the ADN algorithm illustrated in Section 3.3. First, in step ii), the activity of an infected individual i, is multiplied by a parameter p i ∈ [0, 1], similar to [130] . Second, in step iii), a link (i, i h ) to a node i h that is infected is added to the edge set with probability π i ∈ [0, 1]. The epidemic threshold is computed, in terms of the decay ratio to the disease-free equilibrium, finding a rather complicated expression that depends on the joint distributions of the activities and the parameters p i and π i . Such an analytic expression is used to devise an optimization problem, to optimally allocate resources into isolation of infected individuals and awareness, solved via geometric programming. Since its inception in December 2019, the COVID-19 outbreak has rapidly spread, becoming a worldwide pandemic with over 108 million reported cases as of February 16, 2021. In response to this unprecedented health crisis, we witnessed an extraordinary mobilization of the scientific community toward better understanding the novel disease and combating its spread. Within this joint effort, the systems and control communities are playing a key role in developing accurate mathematical models to predict the evolution of the pandemic and assess the effect of diverse intervention policies that have been enacted or that may be implemented [136] [137] [138] . A first, important contribution comes from the formulation and analysis of more complex models for the epidemic progression in fully-mixed populations. These models allow to capture the specific features of COVID-19, including a latency period before the symptoms onset, the presence of asymptomatic individuals, and the implementation of intervention policies such as hospitalization of severe cases and testing. Among these models, we should mention the SIDARTHE model proposed in [18] (see Fig. 11 ), which takes into consideration also the imperfect reporting of infected individuals. In this work, the epidemic model is studied in a fully-mixed population framework by seeing it as a positive linear system under feedback, and the stability of the disease-free equilibrium is thus characterized in terms of the H ∞ norm of its transfer function, which interestingly coincides with the basic reproduction number R 0 . Using this model, an open-loop fast switching strategy to control and suppress the spread of the disease, consisting in intermittent lockdown phases, has been proposed and studied in [139] . Once the epidemic model has been tailored to capture the epidemic progression of COVID-19, its implementation on a network structure (extensively discussed in this survey) is key to predict the spatial spread of the disease, as highlighted in [136] . The time-varying nature of human mobility as well as the implementation of intervention policies through different sequential phases have put epidemic models on dynamic network at the forefront of the stage. In this vein, we mention several data-informed analyses of the outbreak in different countries, using models with regional granularity. These studies include Italy [140, 141] , Ontario, Canada [142] , Western Australia [143] , and Kazakhstan [144] . Based on these network models, non-linear MPC has been adopted to understand the impact of intervention policies and plan their optimal implementation. We mention the works in [145] and in [146] , with case studies based on the outbreak in Italy and Germany, respectively. The impact of local, timevarying lockdown measures and mobility restrictions is analyzed in [147] , where feedback Figure 11 : Schematic of the SIDARTHE model, proposed in [18] to capture the epidemic progression of COVID-19. The modeled is obtained from an SIR model, in which the health state "infected" is substituted by five states, representing all the four possible combinations of whether the individual has symptoms and is detected, and a state for individuals seriously ill; two states are used instead of the "removed" state to account for individuals that are healed (H) or extinct (E), respectively. The arrows illustrate all the possible transitions, which are governed by 16 different parameters. control laws coordinated by a centralized controller are used to design an intermittent and differentiated regional intervention scheme that outperforms nationwide measures. Dynamic networks also enable modeling the individual behavioral response to the pandemic in terms of the adoption of self-protective behaviors and social-distancing measures, which has a huge impact on disease spreading [148] . At the moment of writing this survey, effective pharmaceutical treatments for COVID-19 were unfortunately still not available, while the research for a vaccine has recently led to some promising findings, and extensive vaccination campaigns are getting started. In this delicate phase, in which only a very limited amount of the vaccine is available, public health authorities need to carefully plan the distribution of vaccines. The control problems described and discussed in this survey (e.g., how to optimally modify the recovery rates by distributing a fixed amount of antidote) will be precious tools to help inform vaccination campaigns and distribution of pharmaceutical treatments. In this survey, we went through 260 years of progresses in mathematical models of epidemics. Starting from the very first intuition that mathematics, and in particular systems theory, can provide effective tools toward better understanding the spread of infectious diseases, we reviewed the recent developments and the state of the art in the analysis and control of epidemic models on networks. Alongside, we outlined some of the most promising avenues of current and future research, which are of particular interest for engineering researchers. These directions include the analysis of epidemics on dynamic networks, where the process is modeled as a nonlinear time-invariant system or a complex stochastic process, and the application of control-theoretical tools to devise intervention policies to stop or mitigate the spreading process. Particular attention should be devoted to the analysis and control of stochastic epidemic models, which are attracting attention in computational epidemiology, physics, and network science, but are still mostly overlooked in the systems and control community. Besides these research directions, we believe that our scientific community can provide important advances to other open problems of mathematical epidemiology. A key challenge that the scientific community is currently facing with the COVID-19 outbreak is how to integrate the available data within the mathematical models, that is, how to identify the model parameters from clinical and epidemiological data, which are often noisy and incomplete [149] . Different optimization methods have been proposed to deal with this problem. See, e.g., [150] . From a modeling point of view, an important challenge for the ongoing research, which we briefly mentioned in this survey, is the inclusion of behavioral factors in the model formulation. Game theory has emerged as a promising modeling framework to model human decision-making processes and has been utilized, for instance, to capture the adoption of self-protective behaviors [151, 152] and to characterize the decision whether to vaccinate or not [153] [154] [155] . Essai d'une nouvelle analyse de la mortalité causée par la petite vérole et des avantages de l'inoculation pour la prévenir Influence of Elevation on the Fatality of Cholera On the Mode of the Communication of Cholera The Milroy Lectures on Epidemic Disease in England: The Evidence of Variability and of Persistency of Type A contribution to the mathematical theory of epidemics The problem of endemicity Further studies of the problem of endemicity Dynamics and Control of the Transmission of Gonorrhea Some discrete-time SI, SIR, and SIS epidemic models The analysis of equilibrium in malaria The failure of R 0 On the dynamics of deterministic epidemic propagation over networks The Mathematics of Infectious Diseases Mathematical Models in Population Biology and Epidemiology Modeling influenza epidemics and pandemics: insights into the future of swine flu (H1N1) A double epidemic model for the SARS propagation Understanding the dynamics of Ebola epidemics Modelling the COVID-19 epidemic and implementation of populationwide interventions in Italy Early Prediction of the 2019 Novel Coronavirus Outbreak in the Mainland China Based on Simple Mathematical Model Can the COVID-19 Epidemic Be Controlled on the Basis of Daily Test Reports? The Prevention of Malaria The Mathematical Theory of Infectious Diseases and its Applications Modeling Infectious Diseases in Humans and Animals Epidemic processes in complex networks Analysis and Control of Epidemics: A Survey of Spreading Processes on Complex Networks A Deterministic Model for Gonorrhea in a Nonhomogeneous Population Predicting epidemics on directed contact networks Stability of epidemic models over directed graphs: A positive systems approach Epidemiological Models and Lyapunov Functions Virus Spread in Networks Epidemic Spreading in Real Networks: An Eigenvalue Viewpoint Time-dependent solution of the NIMFA equations around the epidemic threshold Markov Chains and Mixing Times. Providence RI Generalized Epidemic Mean-Field Model for Spreading Processes Over Multilayer Complex Networks The effect of network topology on the spread of epidemics Decay towards the overall-healthy state in SIS epidemics on networks Exponential extinction time of the contact process on finite graphs Epidemics and Rumours in Complex Networks. Cambridge Time to Extinction for the SIS Epidemic Model: New Bounds on the Tail Probabilities Susceptible-infected-susceptible epidemics on networks with general infection and cure times Generalization of Pairwise Models to non-Markovian Epidemics on Networks Burst of virus infection and a possibly largest epidemic threshold of non-Markovian susceptible-infected-susceptible processes on networks Stability of SIS Spreading Processes in Networks With Non-Markovian Transmission and Recovery Equivalence between Non-Markovian and Markovian Dynamics in Epidemic Spreading Processes Equivalence and its invalidation between non-Markovian and Markovian spreading dynamics on complex networks The basic reproduction number in some discretetime epidemic models Discrete epidemic models Analysis, Estimation, and Validation of Discrete-Time Epidemic Processes On the stability of the endemic equilibrium of a discrete-time networked epidemic model The Viral State Dynamics of the Discrete-Time NIMFA Epidemic Model Modeling, estimation, and analysis of epidemics over networks: An overview Discrete-time Markov chain approach to contact-based disease spreading in complex networks On the mixing time of the SIS Markov chain model for epidemic spread Discrete-Time vs. Continuous-Time Epidemic Models in Networks Limitations of discrete-time approaches to continuous-time contagion dynamics Eight challenges for network epidemic models Adaptive coevolutionary networks: a review The dynamic nature of contact networks in infectious disease epidemiology Temporal networks Modern temporal network theory: a colloquium Seasonality and the dynamics of infectious diseases The origin of bursts and heavy tails in human dynamics Burstiness and memory in complex systems Modelling the influence of human behaviour on the spread of infectious diseases: a review Behavioural change models for infectious disease transmission: a systematic review Effect of non-pharmaceutical interventions to contain COVID-19 in China Spread of epidemic disease on networks Epidemic spreading in scale-free networks Emergence of scaling in random networks Virus propagation on time-varying networks: Theory and immunization algorithms Analytical computation of the epidemic threshold on temporal networks Epidemic Threshold of an SIS Model in Dynamic Switching Networks Spectral Analysis of Epidemic Thresholds of Temporal Networks Stability criteria for SIS epidemiological models under switching policies Optimal Design of Switched Networks of Positive Linear Systems via Geometric Programming Temporal interactions facilitate endemicity in the susceptible-infected-susceptible epidemic model Epidemic Threshold in Continuous-Time Evolving Networks Epidemic Processes Over Time-Varying Networks Activity driven modeling of time varying networks Continuous-time discrete-distribution theory for activity-driven networks Contrasting effects of strong ties on SIR and SIS processes in temporal networks Contagion processes on the static and activity-driven coupling networks Epidemic Spreading in Temporal and Adaptive Networks with Static Backbone Epidemic spreading in modular time-varying networks A novel framework for community modeling and characterization in directed temporal networks Epidemic spreading on activity-driven networks with attractiveness Burstiness and Aging in Social Temporal Networks Modeling Memory Effects in Activity-Driven Networks Burstiness in activitydriven networks and the epidemic threshold Simplicial Activity Driven Model Toward epidemic thresholds on temporal networks: a review and open questions Flooding Time of Edge-Markovian Evolving Graphs Modelling approaches for simple dynamic networks and applications to disease transmission models Epidemic threshold and control in a dynamic network Stability of Spreading Processes over Time-Varying Large-Scale Networks Decreasing the spectral radius of a graph by link removals Spectral analysis of virus spreading in random geometric networks Immunization of complex networks Distributing Antidote Using PageRank Vectors Optimal link removal for epidemic mitigation: A two-way partitioning approach Optimization of network protection against virus spread Optimal Resource Allocation for Network Protection Against Spreading Processes Optimal resource allocation for control of networked epidemic models Optimal curing policy for epidemic spreading over a community network with heterogeneous population Distributed resource allocation for control of spreading processes Distributed discrete-time optimization algorithms with applications to resource allocation in epidemics control Distributed Algorithm for Suppressing Epidemic Spread in Networks Asynchronous Distributed Matrix Balancing and Application to Suppressing Epidemic Sparse Resource Allocation for Control of Spreading Processes via Convex Optimization Resource allocation for epidemic network under complications Optimal Resource Allocation to Reduce an Epidemic Spread and Its Complication Optimal and robust epidemic response for multiple networks Optimal Patching in Clustered Malware Epidemics Distributed feedback control on the SIS network model: An impossibility result Applications of the Poincaré-Hopf Theorem: Epidemic Models and Lotka-Volterra Systems Analysis and Control of a Continuous-Time Bi-Virus Model Applying optimal control theory to complex epidemiological models to inform real-world disease management How to distribute antidote to control epidemics An efficient curing policy for epidemics on graphs When Is a Network Epidemic Hard to Eliminate? Suppressing Epidemics in Networks Using Priority Planning Dynamic Control of Modern, Network-Based Epidemic Models Dynamic Resource Allocation to Control Epidemic Outbreaks A Model Predictive Control Approach Robust Economic Model Predictive Control of Continuous-Time Epidemic Processes Optimal control and the value of information for a stochastic epidemiological SIS-model Optimal resource allocation for containing epidemics on time-varying networks Mitigation of epidemics in contact networks through optimal contact adaptation Analysis and Distributed Control of Periodic Epidemic Processes Controlling contagion processes in activity driven networks Effect of individual behavior on epidemic spreading in activity driven networks Analysis and control of epidemics in temporal networks with self-excitement and behavioral changes The spread of awareness and its impact on epidemic outbreaks Optimal Containment of Epidemics over Temporal Activity-Driven Networks Suppression of epidemic spreading in time-varying multiplex networks On Assessing Control Actions for Epidemic Models on Temporal Networks The challenges of modeling and forecasting the spread of COVID-19 Modelling COVID-19 How control theory can help us control Covid-19 Post-lockdown abatement of COVID-19 by fast periodic switching Spread and dynamics of the COVID-19 epidemic in Italy: Effects of emergency containment measures Modelling and predicting the effect of social distancing and travel restrictions on COVID-19 spreading Local lockdowns outperform global lockdown on the far side of the COVID-19 epidemic curve Modelling Strong Control Measures for Epidemic Propagation With Networks-A COVID-19 Case Study A Network-Based Stochastic Epidemic Simulator: Controlling COVID-19 with Region-Specific Policies Model predictive control to mitigate the COVID-19 outbreak in a multi-region scenario Robust and optimal predictive control of the COVID-19 outbreak A network model of Italy shows that intermittent regional strategies can alleviate the COVID-19 epidemic Using social and behavioural science to support COVID-19 pandemic response COVID-19 and the difficulty of inferring epidemiological parameters from clinical data A time-varying SIRD model for the COVID-19 contagion in Italy Spontaneous behavioural changes in response to epidemics Game Theory of Social Distancing in Response to an Epidemic Vaccination and the theory of games Designing Virus-Resistant, High-Performance Networks: A Game-Formation Approach Game-Theoretic Vaccination Against Networked SIS Epidemics and Impacts of Human Decision-Making This work was partially supported by the European Research Council (ERC-CoG-771687) and the Netherlands Organization for Scientific Research (NWO-vidi-14134).