Mathematical analysis of a hybrid model: Impacts of individual behaviors on the spreading of an epidemic

1.   Introduction

The intricate and fascinating relations between the individual and collective behaviors, occurring within a population subject to an evolution problem, have been widely studied for several decades in the field of sociology and economics. Beyond the inescapable textbook of Turner & Killian [37], which is devoted to the study of social collective behaviors, these intricate relations within gatherings, demonstrations and riots have been analyzed, for example in [29] or [38]; self-organization of complex structures, admitting individual or collective social organization, have been studied in [33]; competition between individual and collective behaviors acting in segregation processes have been explored in [18]; more recently, individual roles in collective dynamics or in consensus emergence have been investigated in [26] or in [12]; finally, Galam [17] and Helbing [20] propose to regroup these various questions in a more general framework. However, it is observed that the classical mathematical modelling approaches (especially those using ordinary differential equations) show limitations for reproducing the complex interactions between individual behaviors and collective dynamics. A very simple reason can account for this lack: although ordinary differential equations provide a powerful tool for describing the evolution of population dynamics, it is impossible to follow the trajectory of a single individual within the aggregate flow of such equations. Next, the observation of such real-world phenomena reveals the interleaving of at least two populational scales: the macroscopic scale can be used to describe collective behaviors, whereas the microscopic scale can fit with individuals decisions. Hence, coupling both the macroscopic and microscopic approaches is an attempt to overcome the scaling dilemma. Furthermore, none of both microscopic and macroscopic modelling approaches is better than the other; each one has advantages and disadvantages. However, the differential equations approach offers a formal study potential, which allows to establish qualitative properties of the trajectories determined by the resulting model, at a theoretical level. Thus it appears as a necessity, for the one who aims to advance on this modelling question, to couple ordinary differential equations with a microscopic process. It is precisely the purpose of this work.

Therefore, in this paper, we study a class of hybrid problems, constructed to model the complex features of population dynamic problems, in which the microscopic individual behaviors and the macroscopic collective dynamics are closely intertwined. We consider a general hybrid model, obtained by coupling a system of ordinary differential equations and an agent-based process, which act simultaneously along a common timeline. The complex network structure heavily underlies our hybrid model. Indeed, the subsequent system of ordinary differential equations is embedded into a geographical network structure, so as to reproduce the spatial background of the population dynamics, which possibly present heterogeneous patterns and emergent properties; in parallel, the social interactions occurring between individuals are supported by a social network, which can be partly randomly generated. The resulting mathematical model is analyzed both at a theoretical level and with a numerical approach. The way our model is designed allows its application to a large number of evolution problems arising, for example, in sociology, economics, geography and epidemiology. In this paper, we apply our general model to the current COVID-19 pandemic, whose complex dynamics have been intensely studied.

Besides the general framework of the proposed hybrid models, our approach aims to fill important limitations of previous studies on the topic of epidemics. Indeed, the macroscopic dynamics of epidemics have been intensely studied, using differential equations, very often on the basis of the famous Kermack & McKendrick model [25] (see, for instance, [22], [23], [35] or [36]). Partial differential equations, such as reaction-diffusion equations, have also been used in order to analyze the spatial propagation of epidemics as travelling waves (see for instance [13], [24] or [30]). Recently, epidemiological problems have benefited from the advances on the study of complex networks, which have been used in order to analyze the geographical spreading of various infectious diseases with a different approach (see for instance [7], [9], [27] or [34]). In the mean time, but separately, agent-based models have been used in order to study the microscopic dynamics of epidemics at the individual scale (see for instance [5] or [21]); in particular, agent-based stochastic-process compartmental models have been studied in [2], [3] or very recently in [39]. Afterwards, a first attempt to compare both macroscopic and microscopic approaches is presented in [14], with the aim to establish the importance to take into account the role of individual behaviors. The necessity of coupling differential equations with agent-based models is finally highlighted in ecology (see e.g., [16], [28], [32]) and in epidemiology (see e.g., [1], [4]). However, these studies often lack for a theoretical mathematical validation; for example, the hybrid mathematical model which is studied in the latter paper [4] has not been validated in a theoretical framework and the qualitative properties of its solutions have not been analyzed. Hence, in this article, we do an in-depth study of such hybrid models and fill these temporary deficiencies.

In this paper, we pursue two main objectives. First, we intend to propose a sufficiently large framework which can be applied to a great number of complex problems arising in social and human sciences, and not only in epidemiology, with a rigorous mathematical approach, which supports the numerical simulations and guarantees their relevancy. For that, we construct an abstract class of hybrid models by coupling along a common timeline a macroscopic and continuous process, determined by a system of differential equations, with a microscopic and discrete process, which can be derived from an agent-based model. Under reasonable assumptions, we prove that the resulting hybrid model is well-posed, in the sense that it admits global solutions (see Theorem 2.2), which depend smoothly on a variation of its parameters (see Theorem 2.3). We also prove that the system can exhibit particular solutions with irregular oscillations, (see Theorem 2.4), that generalizes a previous statement presented in [34]. Secondly, we apply our general framework to the current COVID-19 pandemic, so as to prove that the dynamics, spreading and multiple waves of this devastating pandemic cannot be explained only at a macroscopic scale. We use a deterministic ODE compartmental model earlier presented in [35], which is designed for reproducing the specificity of transmission of SARS-CoV-2 virus in a susceptible population, and an agent-based process which focuses on refractory behaviors of citizens to policy strategies. We present several numerical simulations of the hybrid model which highlight the effect of such opposition behaviors on the aggregate flow of the epidemic.

Our paper is organized as follows. In Section 2, we show how to construct the general hybrid model, by coupling a system of differential equations and an agent-based model; we prove that the resulting mathematical problem is well-posed. Next, in Section 3, we apply our general framework to the current COVID-19 pandemic, for which the macroscopic dynamic is described by a deterministic ODE compartmental model, whereas the microscopic opposition behaviors and decision strategies are integrated into an agent-based protocol. The stability analysis of the model is investigated in 4, and a selection of relevant scenarios are finally presented in Section 5 with a numerical approach.

2.   Hybrid models in a general context

In this section, we provide a theoretical mathematical approach to a class of abstract hybrid problems. To construct the hybrid problems, we couple a system of ordinary differential equations with a discrete process, along a common timeline. After, we prove the well-posedness of the resulting mathematical problem and, under reasonable assumptions, the existence of irregular oscillations.

2.1.   Construction of a class of abstract hybrid models

Let us consider a population of individuals and assume that this population is subject to a complex evolution process which cannot be described at a single scale. Thus, we construct a hybrid model by coupling a system of ordinary differential equations and a discrete process, which can be derived from an agent-based model.

First, we assume that the population can be divided into several disjoint groups $x_1$, $\dots$, $x_n$ ($n \geq 1$) and introduce the notation $X = (x_1, \, \dots, x_n)^\top$. Next, we consider a discrete sequence of times

which tends to infinity. We consider the following abstract hybrid problem ($AHP$)

for $s \geq 0$. Here, $F$ is a function defined in $E \times J$ with values in $\mathbb{R}^n$, where $E$ is an open subset of $\mathbb{R}^n$ and $J$ is an open subset of $\mathbb{R}^d$; $G$ is a function defined in $E \times J$ with values in $J$. The first line $( \mathfrak{IC})$ in system $\text{(AHP)}$ determines the initial condition $(X_0, \lambda_0) \in E \times J$; the second line $( \mathfrak{M}_s)$ is an ordinary differential equation which determines the macroscopic part of the hybrid problem, whereas the third line $( \mathfrak{m}_s)$ is a discrete mapping which determines the microscopic part of the hybrid problem.

The hybrid problem $\text{(AHP)}$ can be applied to a great number of population dynamics problems arising in sociology, economics, geography or epidemiology. In particular, we shall study in the next section an application of the hybrid model $\text{(AHP)}$ to the current COVID-19 pandemic; the macroscopic part $( \mathfrak{M}_s)$ will be given by a deterministic ODE compartmental model for the transmission of SARS-CoV-2, whereas the microscopic part $( \mathfrak{m}_s)$ will follows from an explicit agent-based model.

Definition 2.1. For $(X_0, \lambda_0) \in E \times J$, a global solution of the hybrid problem $\text{(AHP)}$ is a pair denoted $\big(X(t, X_0), (\lambda_s)_{s \geq 0}\big)$ in which $X(t, X_0)$ is a function of $t$ defined and differentiable in $[t_0, +\infty)$, such that $X(t_0, X_0) = X_0$ and satisfying $( \mathfrak{M}_s)$ on $(t_s, t_{s+1}]$ for all $s \geq 0$; $(\lambda_s)_{s \geq 0}$ is a sequence of parameters in $J$ satisfying $( \mathfrak{m}_s)$ for all $s \geq 0$.

For brevity, a solution of the hybrid problem $\text{(AHP)}$ starting from $(X_0, \lambda_0)$ shall be denoted as $X(t, X_0, \lambda_0)$.

2.2.   Well-posedness of the hybrid model $\text{(AHP)}$

In order to establish the well-posedness of the hybrid model $\text{(AHP)}$, we consider the three following assumptions.

Assumption 1. The function $F$ involved in the macroscopic part $( \mathfrak{M}_s)$ of the hybrid problem $\text{(AHP)}$ is $\mathscr{C}^1(E \times J)$.

Assumption 2. There exists a compact set $K \subset E$ such that, for all $(X_0, \lambda_0) \in K \times J$, each local solution $X(t, X_0)$ of the Cauchy problem

defined on $[t_0, t_0 + T]$ with $T > 0$ satisfies

for all $t \in [t_0, t_0 + T]$.

Note that in assumption 2, the compact set $K$ does not depend on the initial condition $(X_0, \lambda_0)$.

Assumption 3. The function $G$ involved in the microscopic part $( \mathfrak{m}_s)$ of the hybrid problem $\text{(AHP)}$ is continuous in $E \times J$.

Our first theorem establishes the existence and uniqueness of global solutions to the hybrid problem $\text{(AHP)}$.

Theorem 2.2. Let the assumptions 1 and 2 hold. Then for all $(X_0, \lambda_0) \in K \times J$, the hybrid problem $\text{(AHP)}$ admits a unique global solution $X(t, X_0, \lambda_0)$ defined on $[t_0, +\infty)$, in the sense of Definition 2.1.

Proof. In order to prove the theorem, we construct a solution to the hybrid problem $\text{(AHP)}$ successively on each time interval $[t_s, t_{s+1}]$, with $s \geq 0$.

Let $(X_0, \lambda_0) \in E \times J$. We first consider the Cauchy problem on the time interval $[t_0, +\infty)$:

By virtue of Assumption 1, it follows from the fundamental existence and uniqueness theorem for ordinary differential equations (see for instance Theorem 3.1 in [19]), that the Cauchy problem (2) admits a unique local solution $\tilde{X}(t, X_0)$ defined on $[t_0, t_0 + T]$ with $T > 0$. Furthermore, Assumption 2 guarantees that this solution is global, that is $T = +\infty$. Thus, we define a unique solution of the hybrid problem $\text{(AHP)}$ on $(t_0, t_1]$ by setting $X(t, X_0, \lambda_0) = \tilde{X}(t, X_0)$, $t_0 < t \leq t_1$, $X_1 = X(t_1, X_0, \lambda_0)$ and $\lambda_1 = G\big(X_1, \lambda_0\big)$.

Afterwards, assuming that $X(t, X_0, \lambda_0)$ has been defined on $[t_0, t_s]$ with $s \geq 1$, we set ${X_s = X(t_s, X_0, \lambda_0)}$, $\lambda_s = G\big(X_s, \lambda_{s-1}\big)$ and consider the Cauchy problem on the time interval $[t_s, +\infty)$:

Repeating the above arguments, we construct a unique solution of the hybrid problem $\text{(AHP)}$ on $[t_0, t_{s+1}]$, for all $s \geq 1$. In this way, we have constructed a unique global solution of the hybrid problem $\text{(AHP)}$ in the sense of Definition 2.1.

The next theorem establishes the continuity of the solution $X(t, X_0, \lambda_0)$ of the hybrid problem $\text{(AHP)}$ with respect to a variation of the initial condition $X_0 \in K$ and of the initial parameter value $\lambda_0 \in J$. This continuity principle validates in a rigorous mathematical framework the construction of the hybrid problem, which in turn allows multiple applications to real-world problems.

Theorem 2.3. Let Assumptions 1, 2 and 3 hold. Then each global solution $X(t, X_0, \lambda_0)$ of the hybrid problem $\text{(AHP)}$ is continuous at $(X_0, \lambda_0)$, uniformly on each finite time interval $[t_0, t_0 + T]$ with $T > 0$, that is, for each $T > 0$ and each $\varepsilon > 0$, there exists $\delta > 0$ such that

for all $t \in [t_0, t_0 + T]$, provided that $\left\Vert{(h, k)}\right\Vert_{ \mathbb{R}^n \times \mathbb{R}^d} < \delta$.

Proof. Let us consider $T > 0$, $(X_0, \lambda_0) \in K \times J$ and $\varepsilon > 0$. For each $(h, k) \in \mathbb{R}^n \times \mathbb{R}^d$ such that ${(X_0 + h, \lambda_0 + k) \in K \times J}$, we denote $\tilde{X}_0 = X_0 + h$, $\tilde{\lambda}_0 = \lambda_0 + k$, and $X(t, X_0, \lambda_0)$, $X(t, \tilde{X}_0, \tilde{\lambda}_0)$ the solutions of the hybrid problem $\text{(AHP)}$ stemming from $(X_0, \lambda_0)$ and $(X_0+h, \lambda_0+k)$ respectively. By virtue of Theorem 3.4 in [19], there exists $\delta > 0$ such that

for all $t \in [t_0, t_1]$, provided that $\left\Vert{(h, k)}\right\Vert_{ \mathbb{R}^n \times \mathbb{R}^d} < \delta$. It follows that

where $\tilde{X}_1 = X(t_1, \tilde{X}_0, \tilde{\lambda}_0)$ and $X_1 = X(t_1, X_0, \lambda_0)$, for all $\varepsilon > 0$.

Next we consider $\lambda_1 = G(X_1, \lambda_0)$ and $\tilde{\lambda}_1 = G(\tilde{X}_1, \tilde{\lambda}_0)$. Diminishing $\varepsilon$ and $\delta$ if necessary, by virtue of Assumption 3, we have that $\left\Vert{\lambda_1 - \tilde{\lambda}_1}\right\Vert_{ \mathbb{R}^d} < \xi$, for all $\xi > 0$.

Repeating these arguments a finite number of times, we have that for all $s \geq 1$ such that $t_s \in [0, T]$ and for all $\varepsilon > 0$, there exists $\delta > 0$ such that

for all $t \in [t_0, t_s]$, provided that $\left\Vert{(h, k)}\right\Vert_{ \mathbb{R}^n \times \mathbb{R}^d} < \delta$.

Remark 1. If $G$ is discontinuous, then the solution $X(t, X_0, \lambda_0)$ of the hybrid problem $\text{(AHP)}$ can obviously be very sensitive to a perturbation of the initial parameter value $\lambda_0$. Indeed, let us assume that the differential equation $\dot{X} = F(X, \lambda)$ admits two equilibrium points $\Sigma_1$ and $\Sigma_2$, for two parameter values $\lambda_1$ and $\lambda_2$, respectively. Furthermore, let us assume that $\lambda^*$ is a bifurcation value which separates $\lambda_1$ and $\lambda_2$. If two solutions $X(t, X_0, \lambda_0)$ and $\tilde{X}(t, \tilde{X}_0, \lambda_0)$ of the hybrid problem $\text{(AHP)}$, stemming from two initial conditions $X_0$ and $\tilde{X}_0$, satisfy $\lambda_s < \lambda^* < \tilde{\lambda}_s$ for some $s > 0$, then $X(t, X_0, \lambda_0)$ and $\tilde{X}(t, \tilde{X}_0, \lambda_0)$ will be attracted towards $\Sigma_1$ and $\Sigma_2$, respectively, diverging from each other. Roughly speaking, if the microscopic part $( \mathfrak{m}_s)$ of the hybrid problem $\text{(AHP)}$ is determined by a discontinuous mapping $G$, then its flow can be ripped into several sub-flows which diverge from each other.

2.3.   Existence of particular solutions of $\text{(AHP)}$ exhibiting irregular oscillations

In this section, we establish and interpret an important feature of the hybrid model $\text{(AHP)}$, by proving the existence of particular solutions exhibiting irregular oscillations, under reasonable assumptions which we intend to present now.

Let us suppose that there exist two distinct parameters sets $\Lambda_1 \subset J$ and $\Lambda_2 \subset J$. Those parameters sets possibly satisfy $\Lambda_1 \cap \Lambda_2 = \varnothing$. Assume, moreover, that $\Sigma_1$ is an equilibrium point of the equation

for each $\lambda_1 \in \Lambda_1$, and that $\Sigma_2$ is an equilibrium point of the equation

for each $\lambda_2 \in \Lambda_2$. We introduce the minimum step of the timeline $\lbrace t_s \rbrace$ given by (1) by setting:

The next theorem generalizes a recent result proved in [34].

Theorem 2.4. Suppose that assumption 2 holds. Assume that $\Sigma_1$ is globally asymptotically stable in $W_1 \subset K$ for each $\lambda_1 \in \Lambda_1$, $\Sigma_2$ is globally asymptotically stable in $W_2 \subset K$ for each $\lambda_2 \in \Lambda_2$. Assume moreover that $\Sigma_1 \in W_2$, $\Sigma_2 \in W_1$ and $G(Y, \lambda) \in \Lambda_2$ if $Y$ is near $\Sigma_1$ and $\lambda \in \Lambda_1$, $G(Y, \lambda) \in \Lambda_1$ if $Y$ is near $\Sigma_2$ and $\lambda \in \Lambda_2$.

Then every solution $X(t, X_0, \lambda_0)$ of the hybrid problem starting from $(X_0, \lambda_0) \in W_1 \times \Lambda_1$ admits irregular oscillations, that is, oscillations between a neighborhood $\mathcal{N}_1$ of $\Sigma_1$ and a neighborhood $\mathcal{N}_2$ of $\Sigma_2$, provided the minimum step $\tau$ of the timeline is sufficiently large.

Proof. Let us consider $(X_0, \lambda_0) \in W_1 \times \Lambda_1$. Since $\Sigma_1$ is globally asymptotically stable in $W_1$, the solution of the Cauchy problem

belongs to a neighborhood $\mathcal{N}_1$ of $\Sigma_1$ after a finite time $T > 0$. Moreover, the minimum step $\tau$ of the timeline can be chosen arbitrarily large, thus we can assume that $\tau > T$. In this way, the solution $X(t, X_0, \lambda_0)$ of the hybrid problem $\text{(AHP)}$ belongs to $\mathcal{N}_1$ for $t \in [\tau, \tau+\theta]$ with $\theta > 0$ such that $t_1 \in [\tau, \tau+\theta]$. Hence, the solution $X(t, X_0, \lambda_0)$ of the hybrid problem $\text{(AHP)}$ satisfies $X_1 = X(t_1, X_0, \lambda_0) \in \mathcal{N}_1$.

Next, it is assumed that $\Sigma_1$ belongs to the basin of attraction $W_2$ of $\Sigma_2$; moreover, since $\lambda_0 \in \Lambda_1$, we have $G(Y, \lambda_0) \in \Lambda_2$ if $Y$ is near $\Sigma_1$, thus $\lambda_1 = G(X_1, \lambda_0) \in \Lambda_2$. Since $\Sigma_2$ is globally asymptotically stable in $W_2$, we can repeat the above arguments to conclude that the solution of the hybrid problem $X(t, X_0, \lambda_0)$ reaches a neighborhood $\mathcal{N}_2$ of $\Sigma_2$ after a finite time $t_1 + T$ with $T > 0$, thus satisfies $X(t_2, X_0, \lambda_0) \in \mathcal{N}_2$. Since $\Sigma_2 \in W_1$, it follows that after a finite time $t_2 \in [t_1, t_1+T]$, the solution $X(t, X_0, \lambda_0)$ of the hybrid problem $\text{(AHP)}$ is back in the basin of attraction of $\Sigma_1$.

Since $G(Y, \lambda_2) \in \Lambda_1$ if $Y$ is near $\Sigma_1$, we are brought back to an initial condition $(X_2, \lambda_2) \in W_1 \times \Lambda_1$ with $X_2 = X(t_2, X_0, \lambda_0)$ and $\lambda_2 = G(X_1, \lambda_1)$. Therefore, the solution $X(t, X_0, \lambda_0)$ of the hybrid problem $\text{(AHP)}$ oscillates between $\mathcal{N}_1$ and $\mathcal{N}_2$, which completes the proof.

Remark 2. Obviously, the conclusion of Theorem 2.4 still holds if the initial condition $(X_0, \lambda_0)$ belongs to $W_2 \times \Lambda_2$. Furthermore, if the equation $\dot{X} = F(X, \lambda)$ admits more than two equilibrium points, then the solution $X(t, X_0, \lambda_0)$ of the hybrid problem $\text{(AHP)}$ can exhibit more complex dynamics, with chaotic oscillations between more than two neighborhoods of equilibrium points.

Remark 3. In the sequel, our main application of the hybrid problem $\text{(AHP)}$ is devoted to epidemic models; these epidemic models often admit a disease-free equilibrium and an endemic equilibrium, thus two equilibrium points $\Sigma_1$ and $\Sigma_2$, for two distinct parameter sets which are determined by the value of a basic reproduction number $R_0$. In this context, the existence of pseudo-periodic solutions can be interpreted as the risk of multiple epidemic waves. Assuming that the microscopic part $( \mathfrak{m}_s)$ is partly determined by the decisions of policy makers, the condition "$G(Y, \lambda) \in \Lambda_2$ if $\lambda \in \Lambda_1$ and $Y$ is near $\Sigma_1$" means that the decisions switch the dynamics of the epidemic from one equilibrium to the other. Finally, the condition "$\tau$ is sufficiently large" can be interpreted as a slowness in the decision processes.

3.   Hybrid model applied to COVID-19

In this section, we present an important application of the hybrid framework, constructing an instance of the hybrid problem $\text{(AHP)}$, designed for studying the evolution of COVID-19, whose dynamics cannot be described at a single scale. The macroscopic part $( \mathfrak{M}_s)$ is determined by a deterministic compartmental $SAIRP$ model for the transmission of SARS-CoV-2 with opposition behaviors, whereas the microscopic part $( \mathfrak{m}_s)$ follows from an agent-based model. Both microscopic and macroscopic parts of the problem are supported by a complex network structure.

3.1.   $SAIRP$ model with opposition behaviors

Let us consider a population of individuals affected by the COVID-19 pandemic, caused by the infection of SARS-CoV-2. In order to integrate the impact of refractory behaviors on the dynamics of the epidemic, we propose to improve the deterministic compartmental $SAIRP$ mathematical model, recently analyzed in [34]. This epidemiological model takes into account the transmission of SARS-CoV-2 virus in a population of susceptible individuals.

The followings assumptions are in agreement with the ones from [35]. The population is subdivided into five distinct classes: susceptible individuals ($S$); asymptomatic infected individuals ($A$); active infected individuals ($I$); removed (including recovered and COVID-19 induced deaths) ($R$); protected individuals ($P$). The total population, $N(t) = S(t) + A(t) + I(t) + R(t) + P(t)$, with $t \in [0, T]$ representing the time (in days) and $T > 0$, has a variable size where the recruitment rate $\Lambda$, and the natural death rate $\mu > 0$, are assumed to be constant. The susceptible individuals $S$ become infected by contact with active infected $I$ and asymptomatic infected individuals $A$, at a rate of infection $\beta \frac{\left( \theta A + I \right)}{N}$, where $\theta$ represents a modification parameter for the infectiousness of the asymptomatic infected individuals $A$. Only a fraction $q$ of asymptomatic infected individuals $A$ develop symptoms and are detected, at a rate $v$. Active infected individuals $I$ are transferred to the recovered/removed individuals $R$, at a rate $\delta$, by recovery from the disease or by COVID-19 induced death. A fraction $p$, with $0 < p < 1$, is protected from infection, by the application of non-pharmaceutical interventions, such as, physical distancing, limited size of indoor and outdoor gatherings, teleworking, regular cleaning of frequently-touched surfaces and appropriate ventilation of indoor spaces, mask use and hand washing (see e.g. [8,15]), that prevent from being exposed to the infection, and is transferred to the class of protected individuals $P$, at a rate $\phi$. A fraction $f$ of protected individuals $P$ returns to the susceptible class $S$, at a rate $w$.

In addition to what was proposed in [34], we consider a fraction $u \in [0, 1]$ of individuals who adopt refractory or opposition behaviors to the protection strategy. In this way, the transfer term from the susceptible class $S$ to the protected class $P$ is given by $\pm \phi p(1-u)S(t)$, and the infection term is modified in correspondence. In what follows, for the sake of simplification, we use the notation $\nu = v q$ and $\omega = w f$.

Overall, the $SAIRP$ model with refractory behaviors is determined by the following system of ordinary differential equations:

The latter deterministic compartmental model can be written

with $x = (S, A, I, R, P)^\top \in \mathbb{R}^5$,

and $f(x, \alpha) = \big(f_j(x, \alpha)\big)_{1 \leq j \leq 5}^\top$ with

The equilibrium points and their stability analysis are studied in Section 4.

3.2.   Geographical network modeling the spatial distribution of the population

In order to take into account the geographical distribution of the population which is affected by the epidemic, we propose to embed the latter $SAIRP$ model into a complex network structure. Thus we assume that the individuals are spatially distributed into a finite number of regions $D_1$, $\dots$, $D_m$ with $m \geq 1$; some of these regions are interconnected and individuals present a spatial mobility from one region to another (an example of such a geographical network is depicted in Figure 6 below). Thus we introduce a matrix

of geographical connectivity, by setting $L_{i, k} = 1$ if the region $D_k$ is connected to the region $D_i$ ($i \neq k$); $L_{i, k} = 0$ else; $L_{i, i} = - \sum_{k \neq i} L_{k, i}$ for $i = k$. In this way, the matrix $L$ is a matrix whose sum of coefficients of each column is null.

Next, we divide the population into a finite number of regional sub-populations

where $x_i = (x_{i, j})_{1 \leq j \leq 5}^\top = (S_i, A_i, I_i, R_i, P_i)^\top \in \mathbb{R}^5$ denotes the population of the region $D_i$, $1 \leq i \leq m$.

In this way, the dynamics of the epidemic is modeled at the macroscopic scale by a complex network of ordinary differential equations, which is written

Here, $f = (f_j)_{1 \leq j \leq 5}$ is the function defined by (5) in $\mathbb{R}^5 \times \mathbb{R}^{10}$, with $1 \leq j \leq 5$, and $\alpha_i \in \mathbb{R}^{10}$ is the vector parameter given by (4). The vector-valued function $f = (f_j)_{1 \leq j \leq 5}$ models the dynamic of the epidemic in each region. Since the dynamic is likely to differ from one region to another, the parameters $\alpha_i$ can be distinct for two different regions. The parameters $\sigma_j$ are non-negative coefficients and model the rate of mobility of the regional sub-populations $x_{i, j}$ of type $j$, $1 \leq i \leq m$.

The problem (7) can be rewritten into the short form

with $X = (x_i)_{1 \leq i \leq m} \in \mathbb{R}^{5m}$, $F$ is defined by

and

thus, it is an instance of the macroscopic part $( \mathfrak{M}_s)$ of the hybrid problem $\text{(AHP)}$ (see Section 2), with $n = 5m$ and $d = 10m + 5 + m^2$.

3.3.   Transition from a system of ordinary differential equations towards an agent-based model

After having determined the macroscopic part $( \mathfrak{M}_s)$ of the hybrid model $\text{(AHP)}$, our aim is now to construct the microscopic part $( \mathfrak{m}_s)$. However, two elementary obstacles should be overcome. Indeed, it is worth emphasizing that the differential equation approach is an aggregate approach in which it is impossible to follow the trajectory of a single individual. Furthermore, stemming from an integer number of individuals, the output of the differential equation leads to real numbers which are delicate to interpret. This modeling question has been studied in [4] by prioritizing an important constraint in their case (conservation of population) and therefore by adopting rough approximations, implying a loss of information.

Here, we consider an agent-based process in which agents are determined by non-integer individuals. Let us describe our approach. Assume we have solved system (8) on a finite time interval $[t_s, t_{s+1}]$, with $s \geq 0$, and denote by $X(t, X_s, \lambda_s)$ the corresponding orbit. For each sub-population $x_{ij}$ of type $j$, $1 \leq j \leq 5$, in each region $D_i$, $1 \leq i \leq m$, we evaluate $x_{ij}(t_{s+1})$ and we compute the floor value $N_{ij}(t_{s+1}) = \lfloor x_{ij}(t_{s+1}) \rfloor$. We then consider a finite set

and we call its elements non-integer individuals or simply agents. As constructed, the set $\mathfrak{A}_{ij}$ admits $N_{ij}$ agents, which model the individuals of the sub-population $x_{ij}$ ($1 \leq i \leq m$, $1 \leq j \leq 5$). In parallel, the size of the population at time $t = t_{s+1}$ is stored in the vector $X_{s+1} = X(t_{s+1}, X_s, \lambda_s)$, which is the initial condition for the next instance of the macroscopic part $( \mathfrak{M}_s)$ on the time interval $[t_{s+1}, t_{s+2}]$. The set $\mathfrak{A}_{ij}$ is called the group of agents $\mathfrak{a}_{ij}^g$ in the sub-population of type $j$ in region $D_i$, $1 \leq g \leq N_{ij}$. We then generate a social network over the groups $\big( \mathfrak{A}_{ij})_{1 \leq j \leq 5}$ of each region $D_i$, $1 \leq i \leq m$, by running a graph generation algorithm of edges over the finite set

Here, we choose to generate a Newman–Watts–Strogatz small-world graph, since it is well recognized to reproduce important aspects of the structure of social interactions [31]. However, other graph generation algorithms could be considered. In this social network, each agent has a finite number of neighbors; each agent can observe the types of its neighbors and make decisions, as illustrated in Figure 2. The microscopic part $( \mathfrak{m}_s)$ of the hybrid model $\text{(AHP)}$ is finally determined by a decision protocol which we intend to describe in the following section.

3.4.   Agent-based model approach for integrating the collective behaviors in response to the epidemic

Let us finally describe the agent-based model which is performed at each time step $t_s$ of the time line (1). Assume that the system (8) has been solved on a time interval $[t_s, t_{s+1}]$ with the parameter values

and that the social network of agents $\mathfrak{A}_{ij}$ has been generated at time step $t_{s+1}$.

We assume that agents model citizens or decision makers and focus on two types of actions.

● Action 1. In each region $D_i$, $1 \leq i \leq m$, decision makers evaluate the rate $\rho_I(D_i, t_{s+1})$ of infected individuals; this rate is given by

- If the rate $\rho_I(D_i, t_{s+1})$ overcomes a given threshold $\mathcal{T}_1$, that is $\rho_I(R_i, t_{s+1}) > \mathcal{T}_1$, then decision makers decide to increase the fraction $p_i$ of individuals in region $D_i$ which are protected. This increase is determined by the mapping

where $d_1$ is a decision parameter such that $0 < d_1 < 1$.

- Else, that is, if $\rho_I(D_i, t_{s+1}) \leq \mathcal{T}_1$, then the fraction of protected individuals is maintained, that is

- If at least one of the rates $\rho_I(D_i, t_{s+1})$, $1 \leq i \leq m$, overcomes a second threshold $\mathcal{T}_2 > \mathcal{T}_1$, that is $\rho_I(D_i, t_{s+1}) > \mathcal{T}_2$ for at least one $i \in \lbrace 1, \dots, m \rbrace$, then decision makers decide to confine the individuals in their region, which is determined by setting

- Else, that is, if $\rho_I(D_i, t_{s+1}) \leq \mathcal{T}_2$ for all $i \in \lbrace 1, \dots, m \rbrace$, then decision makers reestablish the last positive value of $\sigma_j$, $1 \leq j \leq 5$.

● Action 2. In each region $D_i$, $1 \leq i \leq m$, agents $\mathfrak{a}_{ij}^g$ of each type $j$, $1 \leq g \leq N_{ij}$, observe the types of their neighbors: $S_i$, $A_i$, $I_i$, $R_i$ or $P_i$. Among these neighbors, each agent $\mathfrak{a}$ evaluates the number $\mathfrak{N}(I, \mathfrak{a}, t_{s+1})$ of infected neighbors. These numbers of infected neighbors are collected into the rate $\rho_ \mathfrak{N}(D_i, t_{s+1})$ of infected neighbors:

where $\mathfrak{A}_i$ denotes the set of agents in region $D_i$, and $N_i$ the number of agents in region $i$.

- If the rate of infected neighbors overcomes a given threshold $\mathcal{T}_3$, that is $\rho_ \mathfrak{N}(D_i, t_{s+1}) > \mathcal{T}_3$, then citizens decide to be in opposition with the protection strategy, which is determined by the mapping

where $d_2$ is an opposition parameter such that $0 < d_2 < 1$. These opposition behaviors mean that policy makers loose the confidence of citizens, who are led to doubt on the protection strategy which is observed as inefficient.

- Else, that is, if $\rho_ \mathfrak{N}(D_i, t_{s+1}) \leq \mathcal{T}_3$, then citizens accept to decrease their level of opposition, which is determined by the mapping

The above protocol defines a discrete mapping

which determines the microscopic part $( \mathfrak{m}_s)$ of a hybrid model. Overall, with the system of ordinary differential equations (8) and the discrete mapping (13) derived from the above agent-based model, we have constructed a hybrid model of the form $\text{(AHP)}$. In the sequel, our aim is to study the dynamics of the COVID-19 hybrid model (8)-(13). We already remark that the continuity assumption 1 on $F$ is obviously satisfied; we shall demonstrate that assumption 2 is also always fulfilled, whereas assumption 3 is likely not to be satisfied.

Remark 4. The two actions protocol determined by (11) and (12) could easily be adapted to other decision strategies. For instance, if the rate $\rho_I(D_i, t_{s+1})$ of infected individuals in some region $D_i$ ($1 \leq i \leq m$) is less than the threshold $\mathcal{T}_1$, then it can be decided that the fraction of protected individuals is decreased, which would be modeled by setting

Similarly, a great number of decision strategies and behavioral changes can be integrated in the protocol, which shows the wide potential of application of our model.

Remark 5. We emphasize that the agent-based model determined by (11) and (12) is not an optimal strategy, since the decisions and their acceptance or opposition behaviors of the agents may produce a negative effect on the dynamic of the epidemic. However, the behavioral agent-based model allows us to introduce more reality on the analysis of epidemics dynamics than the just considering the deterministic compartmental model (8) and, for example, to study the impact of individuals behaviors on the spreading of an infection disease.

4.   Analysis of the COVID-19 hybrid model

In this section, we analyze the dynamics of the COVID-19 hybrid model (8)-(13). We prove that the model satisfies the compactness assumption 2, which guarantees that it admits global solutions. Then we study the local and global stability of its equilibrium states.

4.1.   Uniformly invariant region and global solutions

Let us consider the system of ordinary differential equations (8), where $F$ is given by (9). We introduce the compact set $K \subset \mathbb{R}^{5m}$ defined by

where $\Lambda_0 = \sum_{i = 1}^m \Lambda_i$ and $\mu_0 = \min\big\lbrace \mu_i \; ;\; 1 \leq i \leq m \big\rbrace$. The following theorem guarantees that the compact set $K$ is an invariant region for the system of ordinary differential equations (8).

Theorem 4.1. For any $X_0 \in K$ and any $\lambda_0 \in J$, the Cauchy problem given by (8) and $X(0) = X_0$ admits a unique solution denoted by $X(t, \, X_0, \, \lambda_0)$, defined on $[0, \, \infty)$, whose components are non-negative. Furthermore, the compact set $K$ defined by (14) is positively invariant.

Proof. First, the existence and uniqueness of local in time solutions is immediate, since $F$ is polynomial. Next, the non-negativity of the components is guaranteed by the quasi-positivity of the non-linear operator $F$, which means that it satisfies the property

for all $(x_1, \, \dots, \, x_{m}) \in ( \mathbb{R}^+)^{5m}$, $j \in \lbrace 1, \, \dots, \, m \rbrace$ and $\lambda \in J$. By virtue of Proposition A.17 in [36], it follows that the components of any solution $X(t, X_0, \lambda_0)$ stemming from $X_0$ in $K$ remain non-negative in future time.

Afterwards, it is easily seen that the total population in the geographical network, defined by

satisfies

since the matrix of connectivity $L$ defined by (6) is a zero column sum matrix. By virtue of the Gronwall lemma, it follows that

which leads to the desired conclusion.

The latter statement guarantees that the compactness assumption 2 is fulfilled. Hence, Theorem 2.2 applies to the COVID-19 hybrid model (8)-(13). We obtain the following corollary.

Corollary 1. For any $(X_0, \lambda_0) \in K \times J$, the hybrid epidemiological model determined by (8) and (13) admits a unique global solution $X(t, X_0, \lambda_0)$ in the sense of definition 2.1.

Remark 6. The latter corollary guarantees that the COVID-19 hybrid model (8)-(13) admits global solutions. However, the continuity assumption 3 on the function $G$ given by (13) is very likely not to be satisfied. Indeed, even if the graph generation algorithm over the agents set (10) satisfies the continuity requirement (which would mean that two social networks are arbitrarily closed to each other, provided their corresponding agents sets are sufficiently closed), the decision protocol determined by (11) and (12) involves thresholds which can lead to discontinuities of the function $G$. In this case, the solutions of the COVID-19 hybrid model (8)-(13) can exhibit sensitive dynamics, as will be investigated in Section 5. This sensitivity of the system should not be seen as a pathology of the model, since it reproduces the complex dynamics of epidemics which can be very difficult to predict.

4.2.   Equilibrium points and basic reproduction number of the $SAIRP$ model with opposition behaviors

Here we study the equilibrium states of the epidemiological model (3), which determines the local dynamics on each node of the complex network given by (7).

Basic computations show that the model (3) has two equilibrium points:

● disease-free equilibrium, denoted by $\Sigma_0$, given by

● endemic equilibrium, $\Sigma_+$, whenever $R_0 > 1$, given by

with

where the basic reproduction number, $R_0$, is given by

In what follows, we analyse the influence of the parameters $p$ and $u$ on the values of $R_0$. The partial derivatives of $R_0$ with respect to $p$ and $u$ are given by, respectively,

As all parameter values assume positive values and $0 < p, u < 1$ we have $\frac{\partial R_0}{\partial u } > 0$ and $\frac{\partial R_0}{\partial p } < 0$. Meaning that $R_0$ increases with $u$, but decreases as $p$ increases.

In order to observe this variation of $R_0$ for different values of $p$ and $u$, let us assume, fixed values for the parameters (see Tables 1-2):

Then,

Let $0.25 \leq p \leq 0.675$ and $0 \leq u \leq 0.4$, then there holds $0.4147461 \leq R_0(p, u) \leq 4.5124381$, see Figure 3.

4.3.   Stability analysis of the equilibrium points

The Jacobian matrix of system (3), evaluated at the disease-free equilibrium (15), is given by

with

The eigenvalues of the matrix $M\left(\Sigma_0\right)$ are given by $\lambda_1 = \lambda_2 = -\mu$, $\lambda_3 = -(\phi p (1-u) + \mu +\omega)$ and the remaining two, $\lambda_4$ and $\lambda_5$, are the roots of the polynomial $P(\lambda)$ given by

where $B = \frac{-\beta \theta (1-p(1-u)) (\omega+\mu)}{(p(1-u) \phi+\mu+\omega)} +\delta+2 \mu+\nu$ and $C = \frac{\mathcal{D-\mathcal{N}}}{p (1-u)\phi+\mu+\omega}$. The local stability of the disease-free equilibrium $\Sigma_0$ comes by the application of the Routh–Hurwitz criterion, if, and only if, $B > 0$ and $C > 0$. It is easy to show that $C > 0$ whenever $R_0 < 1$. The proof that $B>0$ is positive whenever $R_0$ is analogous to the one made in [34] replacing $p$ by $p(1-u)$. We have just proved the following theorem.

Theorem 4.2. The disease free equilibrium of the $SAIRP$ model (3), $\Sigma_0$ given by (15), is locally asymptotically stable whenever $R_0 < 1$.

To analyze the local stability of the endemic equilibrium, $\Sigma_+$, we evaluate the Jacobian matrix of system (3) at the endemic equilibrium (17), here denoted by $M\left(\Sigma_+\right)$ (we do not write the expression of $M\left(\Sigma_+\right)$ due to the long size of each terms). We follow by computing the characteristic equation, given by,

where $I_5$ represents the identity matrix of dimension 5.

The roots of the equation $Q(t) = 0$ are however hard to analyze analytically. Therefore, here we study the local stability of $\Sigma_+$ considering fixed parameter values (Pfixed). For these parameter values, the characteristic equation $Q(t) = 0$ is given by

with

Applying the Routh–Hurwitz criterion to the fourth order polynomial (19), the endemic equilibrium $\Sigma_+$ is locally asymptotically stable whenever the following inequalities hold:

The conditions $A_0 > 0$ and $A_3 > 0$ hold for all positive $0 < p, u \leq 1$. As for the conditions $A_3 A_2 - A_4 A_1 > 0$ and $C > 0$, we assume the parameter values (Pfixed) and observe that they hold whenever $R_0(u, p) > 1$, for $0.25 \leq p \leq 0.675$ and $0 \leq u \leq 0.4$, see Figure 4.

Remark 7. Note that the global stability of the disease free and endemic equilibrium points $\Sigma_0$ and $\Sigma_+$ can also be established; the proofs are analogous to the ones made in [34], replacing $p$ by $p(1-u)$.

4.4.   Fit the $SAIRP$ model with COVID-19 data in Portugal

The $SAIRP$ model with piecewise constant parameters presented in [34], can be used to fit the real data of active infected individuals with COVID-19, for a fixed period of time and for a single region. Here we consider the number of active infected individuals with COVID-19 in Portugal from the date of the first confirmed case, March 2, 2020, until April 15, 2021 (the data can be found in [11]). We consider a time interval of 410 days which we subdivide into 9 sub-intervals based on the behavior of the epidemic curve and the public health policies imposed by the Portuguese Government. For each sub-interval of time, the parameters $\beta$, $f$ and $p$ take piecewise constant values and we assume $u = 0$. See Table 1 for the time sub-intervals and the corresponding piecewise constants parameter values. The remaining parameter values are fixed and take the values presented in Table 2.

However, the piecewise constant parameters model does not integrate at the microscopic scale the impact of opposition behaviors to policy strategies; moreover, it does not take into account the geographical distribution of individuals. Thus our aim in the next section is to deeper explore the dynamics of the hybrid model (8)-(13), with a numerical approach.

5.   Numerical simulations of the COVID-19 hybrid model

In this section, our aim is to explore numerically the dynamics of the hybrid model (8)-(13), in order to show how and to what extent individual behaviors can influence the dynamic of the epidemic. The hybrid problem has been implemented with the $\texttt{python}$ language, using the library $\texttt{networkx}$, which integrates graph generation algorithms. We consider an abstract geographical network of 5 regions (depicted in Figure 6), which is inspired from the complex network studied in [34]; the connectivity matrix of this geographical network is given by

where we assume a weak and homogeneous connectivity by setting $\varepsilon_{ij} = 0.1$ for all $i, j \in \lbrace 1, \dots, 5\rbrace$. The coupling strengths are given by $\sigma_S = \sigma_A = 0.1$ and $\sigma_I = \sigma_R = \sigma_P = 0$, which means that only susceptible and asymptomatic individuals present a spatial mobility. The values of other parameters are given in Table 3. We have chosen these parameter values so that the endemic equilibrium is stable in the regions 1, 2, 3, whereas the disease free equilibrium is stable in the two other regions. In other words, in absence of couplings in the geographical network, the orbits of the model (3) are attracted to $\Sigma^+$ in three regions, whereas they are attracted to $\Sigma_0$ in two other regions. Each social network is generated by running a Watts-Strogatz small-world algorithm in which each agent has $10$ neighbors at most, and where the probability of rewiring each edge is fixed to $0.02$. In this way, we generate a set of agents $\mathfrak{A}_{ij}$, as detailed in Subsection 3.3, whose structure reproduces the main social interactions in a group of individuals (see [31]). Now we investigate the effects of individual behaviors on the dynamics of the hybrid model, by testing four relevant scenarios, which are presented and discussed below. To this end, we solve the hybrid model (8)-(13) on a fixed timeline and store the outputs of the agent-based process (11)-(12) in data files. The results of the numerical simulations are illustrated in Figure 7: each sub-figure shows the number of infected individuals $I_i(t)$ of one region $D_i$ ($1 \leq i \leq 5$), for the four selected scenarios, and a comparison with the simulation of the "only-macroscopic" model (8). In these four scenarios, we fix the end of the timeline (1) at $T = 50$ and we focus on a variation of the number $N$ of steps along this timeline, the thresholds $\mathcal{T}_1$, $\mathcal{T}_2$, $\mathcal{T}_3$ and the decision parameters $d_1$, $d_2$ involved in the agent-based model given by (11) and (12). Obviously, all parameter values could be modified for testing other relevant scenarios. For instance, the influence of the topology of the geographical network has been discussed (for a non-hybrid model) in [7].

5.1.   Scenario 1: infection under control with fast confinement decision

In the first scenario, we divide the timeline (1) of the hybrid model in $N = 60$ sub-intervals of constant width. The decision parameters and the thresholds involved in the agent-based protocol (11)-(12) are given by

These numerical values have been chosen arbitrarily and could easily be modified. With those parameter values, we observe, as depicted in green in Figure 7, that the numbers of infected individuals $I_i(t)$ decrease to $0$ in each region $D_i$ ($1 \leq i \leq 5$). More specifically, the output of the agent-based process leads to a first decision of confinement (induced by Action 1 in the protocol (11)-(12)) from $t = 0$ until $t = 0.8$ and a second decision of confinement from $t = 4.9$ until $t = 8.2$. These decisions of confinement are favored by a low value of the threshold $\mathcal{T}_2$. In parallel, the weak value of the parameter $d_2$ means minority opposition behaviors in the population. By comparison with the output of the only-macroscopic model given by the system (8) (which is depicted in red in Figure 7), we observe that the decision protocol with fast confinement decision, succeeds in avoiding the persistence of the disease in the whole population. We also note a particular situation in region $5$: the hybrid model does not modify the decrease of infection to $0$, which is guaranteed by a low basic reproduction number ($R_0 \simeq 0.0375$, see Table 3) and also the geographical position of region $5$, which is not connected to regions $1$, $2$, $3$, where the basic reproduction number is greater than $1$.

5.2.   Scenario 2: postponed extinction of the disease by slowness in the decision process

In the second scenario, we divide the timeline (1) in $N = 20$ sub-intervals of constant width. In this way, the first execution of the agent-based process (11)-(12) occurs at $t = 2.5$, whereas it occurs at $t \simeq 0.83$ in the first scenario. This choice roughly models a slowness decision process. Other parameters are still given by

As depicted in blue in Figure 7, we observe that the numbers of infected individuals $I_i(t)$ still decrease to $0$ in each region $D_i$ ($1 \leq i \leq 5$). However, the decrease is postponed with respect to scenario 1. Indeed, the output of the agent-based model (11)-(12) leads to a single decision of confinement, from $t = 4.8$ until $t = 14.3$. Overall, the decision of confinement is made later and its duration is longer than in scenario 1, which leads to a postponed extinction of the disease. It is also interesting to observe the emergence of small waves in region $4$; these small waves seem to be caused by the connections with regions $1$, $2$, $3$, where the decrease of infection is postponed.

5.3.   Scenario 3: postponed extinction of the disease by opposition behaviors

In the third scenario, our aim is to explore the possible conflict between a slow decision process of policy makers (as described by Action 1 in the agent-based protocol (11)-(12)) and a high level of opposition behaviors of citizens (determined by Action 2 in the protocol (11)-(12)). To that aim, we still divide the timeline (1) in $N = 20$ sub-intervals of constant width; the values of $d_1$, $\mathcal{T}_1$ and $\mathcal{T}_2$ are the same as in scenarios 1 and 2, whereas we modify the values of the opposition behaviors parameters by setting

This choice roughly models an increase of opposition behaviors in the population and a lower threshold of acceptation of protection strategies. As depicted in dark yellow in Figure 7, we observe that the infection is still under control, but the extinction of the disease is almost postponed by the opposition behaviors. Moreover, we remark that the small infection waves, which were already observed in region $4$ with scenario $2$, are amplified.

5.4.   Scenario 4: emergence of multiple pandemic waves due to risk negation

In the fourth scenario, we divide the timeline (1) in $N = 60$ sub-intervals of constant width. We come back to minority opposition behaviors as in scenarios $1$ and $2$ by setting $d_2 = 0.01$, $\mathcal{T}_3 = 0.3$, but we modify the decision of confinement parameters by assuming

The value of the threshold $\mathcal{T}_1$ is greater than in other scenarios, which means that the decision of protection are postponed; in parallel, the value of the threshold $\mathcal{T}_2$ is much greater, which means that the control strategy does not consider seriously confinement as a major decision, or even that confinement should be avoided, for example in order to maintain the economic activity, at the risk to exacerbate the level of infection. In this case, we observe, as depicted in black in Figure 7, that the level of infection does not decrease to $0$, except in region $5$, which is not connected to the regions where the basic reproduction number is greater than $1$. We observe the emergence of infection oscillations around a level which corresponds to the persistence of the disease in the non-hybrid model (8). The emergence of such oscillations was predicted by Theorem 2.4. It suggests that the negation of the health risk by decision makers could compromise the control of the epidemic. We also indicate that a small perturbation of the parameters $d_1$ and $\mathcal{T}_1$ can amplify or destroy these oscillations, which might correspond to a bifurcation process in the hybrid model (8)-(13).

6.   Conclusion and future work

In this paper, we have brought an original contribution to the study of complex systems arising in social science, economics and epidemiology, by constructing a class of hybrid models, in which the macroscopic dynamics of a population subject to an evolution problem, were coupled with the microscopic dynamics of individuals. The macroscopic dynamics was modeled by a system of differential equations, embedded in a geographical network structure, whereas the microscopic dynamic was modeled by an agent-based process, which can integrate various individual behaviors. The transition between the macroscopic and the microscopic scales involves the generation of a social network, which reproduces the social interactions occurring in the population. Our hybrid model was studied in an abstract and theoretical framework, by establishing the existence and uniqueness of relevant solutions, their continuous dependence with respect to a variation of its parameters and the possible emergence of pseudo-periodic solutions. We applied our hybrid model to the study of the current COVID-19 pandemic, and performed several numerical simulations, which highlight how microscopic behaviors can have a strong impact on the macroscopic of the epidemic.

As future work, we aim to generalize our hybrid framework, by considering a larger class of dynamical systems modeling the macroscopic dynamic of the population: for example, non-autonomous equations, delay differential equations or reaction-diffusion equations could be considered. Analogously, the agent-based process studied in Section 3, can be subject to various generalizations, which show the wide potential of the hybrid approach. Moreover, it is also a challenge to improve the SAIRP model by considering specific subgroups for the protected class $P$ distinguishing the different protective measures that may be applied to control the epidemic. We also aim to investigate a question which was left open in the present paper: the emergence of unstable oscillations under an increase of opposition behaviors has been observed, and will be deeper investigated in a separate work.

Acknowledgments

This work was partially supported by Portuguese funds through CIDMA, The Center for Research and Development in Mathematics and Applications of University of Aveiro, and the Portuguese Foundation for Science and Technology (FCT–Fundação para a Ciência e a Tecnologia), within project UIDB/04106/2020. Silva is also supported by the FCT Researcher Program CEEC Individual 2018 with reference CEECIND/00564/2018.

The authors are sincerely grateful to the anonymous reviewers for their valuable comments and suggestions which improved the presentation of the paper.