Complex network model for COVID-19: Human behavior, pseudo-periodic solutions and multiple epidemic waves.

We propose a mathematical model for the transmission dynamics of SARS-CoV-2 in a homogeneously mixing non constant population, and generalize it to a model where the parameters are given by piecewise constant functions. This allows us to model the human behavior and the impact of public health policies on the dynamics of the curve of active infected individuals during a COVID-19 epidemic outbreak. After proving the existence and global asymptotic stability of the disease-free and endemic equilibrium points of the model with constant parameters, we consider a family of Cauchy problems, with piecewise constant parameters, and prove the existence of pseudo-oscillations between a neighborhood of the disease-free equilibrium and a neighborhood of the endemic equilibrium, in a biologically feasible region. In the context of the COVID-19 pandemic, this pseudo-periodic solutions are related to the emergence of epidemic waves. Then, to capture the impact of mobility in the dynamics of COVID-19 epidemics, we propose a complex network with six distinct regions based on COVID-19 real data from Portugal. We perform numerical simulations for the complex network model, where the objective is to determine a topology that minimizes the level of active infected individuals and the existence of topologies that are likely to worsen the level of infection. We claim that this methodology is a tool with enormous potential in the current pandemic context, and can be applied in the management of outbreaks (in regional terms) but also to manage the opening/closing of borders.

Introduction

The history of the pandemics that have already plagued mankind has revealed that there are consistently periods of marked increase in infected individuals, followed by phases in which the numbers are relatively lower. In these cases, there is also a repetition of these oscillations that are called pandemic waves. A second wave of the pandemic poses an imminent threat to society, with an immense cost in terms of human lives and a devastating economic impact, it is therefore, crucial to try to avoid the emergence of pandemic waves.

The coronavirus disease 2019 (COVID-19) is an infectious disease caused by the severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2). The World Health Organization announced the COVID-19 outbreak as a pandemic on 11 March 2020 [29]. At the time of writing, several countries across Europe are seeing a resurgence in COVID-19 cases after successfully controlling the first outbreak and started to take action to face the so-called second wave [12]. Mathematical tools have been important in the analysis and control of COVID-19 pandemic, see e.g. [16], [18], [19], [22] and references cited therein.

In this paper, we propose a mathematical model given by a system of ordinary differential equations based on the model from [23], for the transmission dynamics of SARS-CoV-2 in a homogeneously mixing population, and we generalize it to a non constant population model with piecewise constant parameters. Considering parameters determined by piecewise constant functions allows to model the governmental and public health decisions of political actors, which have a large influence on the behaviors of individuals, which in turn can change the dynamics of the epidemic, see e.g. [1], [3], [11]. Moreover, the piecewise constant parameters also allow to mathematically model the human behavior in the application of non-pharmaceutical interventions (NPI), 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, avoiding close contact and hand washing, see e.g. [8], [13].

Viruses need cells to replicate themselves, as they are not able to do so on their own. Depending on the type of virus, contagion is more or less facilitated depending on the means of transmission. In the case of SARS-COV-2, the transmission is carried out through droplets containing viruses that are sent from the respiratory tract. Accordingly, as the airways are the gateway to the virus, the proximity of people and the dynamics imposed by the movement are a determining factor in the spread of the COVID-19 disease. The proof that the movement of people is a determining factor for contagion is that during the recent periods of confinement the number of infections has decreased very markedly [27].

In this paper, the COVID-19 mathematical model is improved by constructing a complex network of dynamical systems, following the framework presented in [6], [7], in order to take into account the mobilities of individuals, which are also known to play a decisive role in the dynamics of the epidemic. We investigate in which cases such a COVID-19 second wave can occur, by establishing sufficient conditions for pseudo-periodic solutions: see, e.g., [9], [15], [28] for traveling waves results in epidemiological models and [2], [4], [5], [14] for complex network models.

Our model is calibrated in order to fit with the real data of the COVID-19 dynamics in six regions of Portugal mainland, namely Norte, Centro, Lisboa e Vale do Tejo, Alentejo, Algarve and Pinhal Litoral.

Numerical simulations are provided where we explore the effect of the topology of the network on the dynamics of the epidemics (disposal of connections and coupling strength). We identify which type of topology minimizes the level of infection of the epidemic, and which type of topology worsens the number of infected individuals.

This paper is organized as follows. In Section 2, we propose a mathematical model for the transmission dynamics of SARS-COV-2, with constant parameters and variable population size. We show that the model admits two equilibrium points and we analyze their local and global stability. In Section 3, we consider piecewise constant parameters, which allows to model the impact of public health policies and the human behavior in the dynamics of the COVID-19 epidemic. The existence and uniqueness of global solutions of the model with piecewise constant parameters is proved. An important result in the context of the COVID-19 pandemic and the resurgence of epidemic waves, is proved in Section 4, where a sufficient condition is proved to the existence of pseudo-periodic solutions. In Section 5, we construct a complex network of models with piecewise constant parameters, where each node represents one of six regions in Portugal, and where the values of the parameters differ from one region to another. The case study, with COVID-19 real data from Portugal, is analyzed in Section 6, where we calibrate the model to each of the six regions and after, numerical simulations are performed in order to determine the topology that minimizes the average number of active infected individuals in Portugal. The results obtained in this paper are discussed, in Section 7, from a practical point of view and its implications in the management of the COVID-19 pandemic waves. We end the paper in Section 8, with some conclusions and future work.

Model with constant parameters

We propose a compartmental mathematical model, based on [23], where the population is subdivided into five 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, , with representing the time (in days) and , has a variable size where the recruitment rate, Λ, and the natural death rate, , are assumed to be constant. The susceptible individuals S become infected by contact with active infected I and asymptomatic infected A individuals, at a rate of infection , where θ represents a modification parameter for the infectiousness of the asymptomatic infected individuals A. The remaining assumptions follow the ones from [23]. 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 δ, by recovery from the disease or by COVID-19 induced death. A fraction p, with , is protected (not immune) from infection, by the application of non-pharmaceutical interventions (NPI), 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], [13], that prevent from being exposed to the infection, and is transferred to the class of protected individuals P, at a rate ϕ. A fraction m of protected individuals P returns to the susceptible class S, at a rate w. In what follows, for the sake of simplification, we use the notation and . The previous assumptions are described by the following system:

Existence, positivity and boundedness of solutions

The equations of the SAIRP model (1) can be rewritten as with and , where the non-linear operator f is defined in by

In order to prove that the problem determined by (2) is well-posed, we introduce the compact region defined by

The following theorem establishes the existence of global solutions to (2).

For any , the Cauchy problem given by (2) and admits a unique solution, denoted by , defined on , whose components are non-negative. Furthermore, the region Ω defined by (4) is positively invariant.

The existence of a local in time solution to problem (2) starting from follows from the theory of ordinary equations (see, for instance, [20]). The non-negativity of the components is guaranteed by the quasi-positivity of the non-linear operator , which means that it satisfies the property for all , and . By virtue of Proposition A.17 in [24], it follows that the components of any solution stemming from in Ω remain non-negative in future time. Finally, summing the five equations of system (2) leads to from which it is deduced, using Gronwall's lemma, that The latter inequality proves that the region Ω defined by (4) is positively invariant. □

A possible alternative for establishing the existence of a positively invariant region is to apply the method described in [25, Chapter 14]: if the vector field corresponding to the system points into the interior of the region, then the region is seen to be positively invariant.

Note also that the upper bound involved in the definition of the region Ω can be physically interpreted: if the mortality rate μ increases, then the diameter of the region Ω decreases, which roughly speaking means that the population described by the system is smaller. In parallel, if Λ increases, then the region Ω is enlarged, which means that a higher initial recruitment guarantees a numerous population.

Equilibrium points and basic reproduction number

In this section, we study the equilibrium points of system (1) and derive their expressions with respect to the parameters of the model. We compute the basic reproduction number, denoted by , following the approach in [26].

The model (1) has two equilibrium points: where the basic reproduction number, , is given by

disease-free equilibrium, denoted by , given by

endemic equilibrium, , whenever , given by with

The computation of the equilibrium points and the basic reproduction number follows standard arguments, see [26]. □

Local and global stability analysis

In what follows we prove the local and global asymptotic stability of the disease-free equilibrium (DFE), when the basic reproduction number satisfies the inequality . For , we prove that the endemic equilibrium is globally asymptotically stable in a biologically meaningful compact positive invariant region.

local stability of the DFE

The disease-free equilibrium, , is locally asymptotically stable whenever .

For the sake of simplicity, let denote the numerator and denominator, respectively, of the basic reproduction number, , given by (8). The Jacobian matrix of system (1), evaluated at the disease-free equilibrium (5), is given by The eigenvalues of the matrix are given by , and the remaining two, and , are the roots of the polynomial given by where and . Applying the Routh–Hurwitz criterion, we conclude that model (1) is locally stable if, and only if, and . It is easy to show that whenever . The coefficient B is positive when Since all parameters take positive values, we have From , we have and, therefore, From the biological meaning of δ, μ and ν, we have that (observe that their unit is ). Therefore, and that is, Multiplying the previous inequality by we have In other words, Therefore, which ensures that . It follows the local asymptotic stability of the disease-free equilibrium . □

global stability of the DFE

If , then the disease-free equilibrium, , is globally asymptotically stable in Ω.

Since , we can write with . We obtain that which leads to where k is defined by A fortiori, we have: Now we consider the following functional given by where ζ and ξ are defined by and χ is a positive constant which will be determined below. Recall that , and note that by virtue of (10). As constructed, L is a non-negative functional and we have Now we compute the derivative of L along the solutions of the SAIRP model (1) starting in Ω. We have Now we use the relations to obtain First step. Let us examine the terms involving and . We have Choosing ξ such that , the latter equality becomes Elementary computations show that Now we use the standard inequality to conclude that . Thus, we obtain Second step. We examine the terms involving A and I. After simplifications, we have Next, we write By virtue of Young's inequality, we have where ε is a positive coefficient which can be chosen arbitrarily small. It follows that Since , we obtain Similar computations show that We obtain that Now we choose which yields . Since , we obtain, by virtue of (11), Analogously, it is seen that by virtue of (9). After those simplifications, we obtain since . Finally, we choose sufficiently large so that which guarantees that . In other words, the functional L is a Lyapunov function for the flow induced by the SAIRP model (1). The conclusion follows from LaSalle's invariance principle [17]. □

The next theorem establishes the global stability of the endemic equilibrium (EE) defined by (7), in an invariant region . It remains an open question to determine if the endemic equilibrium is globally asymptotically stable in the whole region Ω.

global stability of the EE

The compact region Γ defined by is positively invariant under the flow induced by system (1) . It contains the disease-free equilibrium, , and the endemic equilibrium, , if . Furthermore, if , then the endemic equilibrium is globally asymptotically stable in Γ.

Let us first prove that Γ is positively invariant under the flow induced by system (1). We denote by the solution of system (1) starting from . Summing again the five equations of system (1) leads to where . Consequently, we have Now we have since . It follows that for all , which proves that Γ is positively invariant. Next, we easily verify that and belong to Γ. Now we turn to prove the global stability of in Γ. To that aim, we introduce the functional V defined by with positive coefficients , , , and satisfying where . As constructed, V is a non-negative functional and we have We compute the time derivative of V along the solutions of system (1) starting in Γ. Using the fact that is an equilibrium of system (1), we obtain, after simplifications, that We perform elementary computations, so as to split into four terms, as follows: with First, similar computations as in the proof of the global stability of the disease-free equilibrium show that , since . Next, it is seen that , by virtue of the elementary inequality Afterwards, remarking that , algebraic computations show that with We introduce , and . It is observed that . By virtue of the standard inequality we deduce that . Finally, using (7) and (14), we show that . Gathering the above results shows that . Therefore, we have proved that the functional V is a Lyapunov function for the flow induced by the SAIRP model (1). The conclusion follows once again from LaSalle's invariance principle [17]. □

The theoretical results proved in Theorem 3, Theorem 4 are illustrated in Fig. 1 .

Fig. 1

Phase portraits in the (S, A, I) space, illustrating the global stability of the equilibrium points. (a) Global stability of the disease-free equilibrium (R0 < 1). (b) Global stability of the endemic equilibrium (R0 > 1).

In the critical case , it is easily observed, by virtue of equation (7), that the equilibrium points and coincide. In that case, is not hyperbolic and looses its stability, thus it might belong to a compact invariant set, while not being asymptotically stable.

Model with piecewise constant parameters

The human behavior and the governmental public health decision makers have a strong influence on the dynamics of the epidemics and it appears necessary to take into account these complex effects into the mathematical model, albeit it seems very difficult to accurately measure the impacts of these behaviors on the model parameters, especially on the basic reproduction number . In this section, we model the human behavior and the impact of the decision makers policies, by considering in model (1) parameters determined by piecewise constant functions. We prove the existence and uniqueness of global solutions of the resulting model. We start by subdividing the time line into a finite number of n intervals with disjoint unions, and introduce a piecewise constant function α defined on each time interval as with , and . Next, we consider the sequence of Cauchy problems defined for each initial condition by

We are now in a position to derive the existence and uniqueness result.

For any initial condition , the sequence of Cauchy problems given by (15) admits a unique global solution, denoted again by , whose components are non-negative. Furthermore, the region Ω is positively invariant.

Applying Theorem 1, a finite number of times, directly provides the existence and uniqueness of global solutions to problem (15). □

Note that the solutions of problem (15) are continuous on the time interval , but may not be of class at , . From the modeling point of view, each change of parameters occurring at time () corresponds, for example, to a public announcement of confinement/lift of confinement or prohibition of displacement.

Existence of pseudo-periodic solutions

In this section, we show that piecewise constant parameters can lead to pseudo-periodic solutions. This result has important practical applications in the context of the COVID-19 pandemic. During the first months of the COVID-19 pandemic, the first concern of governments was to decrease the level of infected individuals. This was achieved, in great part, by confinement decisions. However, those confinement decisions had a large impact on the economy, and it was important to decide to lift this ban. Consequently, the risk of a premature relaxation was the second concern of governments, since it was feared to provoke a second wave of infection. The existence of multiple epidemic waves in a pandemic can be mathematically justified by the following theorem.

Assume that the disease-free equilibrium, , admits a non-trivial basin of attraction if , and that the endemic equilibrium admits a non-trivial basin of attraction if . Let and denote two sets of parameters of system (1) such that and . Let and consider the sequence of Cauchy problems for , where is such that , for i odd, and , for i even. Then the solution , of the latter sequence of Cauchy problems, exhibits pseudo-oscillations between a neighborhood of and a neighborhood of in Ω.

The initial condition of the sequence of Cauchy problems (16) has been chosen in , which is assumed to be the basin of attraction of . Since , then the solution is attracted to and thus reaches a neighborhood of after a finite time . Now it is assumed that belongs to the basin of attraction of . Since , then the solution is now attracted to and thus reaches a neighborhood of after a finite time . Next, we assume that belongs to the basin of attraction of . Since , then the solution is now attracted to and thus reaches again the neighborhood of after a finite time . Repeating those arguments a finite number of times leads to the desired conclusion, that is, the solution exhibits pseudo-oscillations between the neighborhoods and of and , respectively. □

Note that the assumption is directly satisfied for large enough, if belongs to the interior of in Ω. Similarly, the assumption is directly satisfied for large enough, if belongs to the interior of in Ω. By virtue of Theorem 3, the latter property is directly satisfied since . A numerical simulation of a solution of the sequence of Cauchy problems (16), exhibiting pseudo-oscillations between a neighborhood of and a neighborhood of in Ω, is depicted in Fig. 2 .

Fig. 2

Solution of the sequence of Cauchy problems (16) exhibiting pseudo-oscillations between a neighborhood of Σ0 and a neighborhood of Σ+ in Ω. Here I denotes the number of active infected individuals and p the fraction, 0 < p < 1, of susceptible individuals S that is transferred to the protected class P.

Dynamics of a complex network of non-identical SAIRP models

It is now widely admitted that mobilities play an important role on the dynamics of epidemics, at numerous stages of their development. In this section, we propose to study the propagation of the COVID-19 outbreak in Portugal by modeling this country by a complex network in which the six regions studied previously for the calibration of the model (1) are considered.

In a first stage, we show how to construct a complex network of SAIRP models. Let us consider the six regions of Portugal, as depicted in Fig. 3 . Those six regions are connected by a finite number of links that define a graph made of a set of 6 vertices, which correspond to the six regions (Norte, Centro, Lisboa e Vale do Tejo, Alentejo, Algarve, Pinhal Litoral), and of a set of edges, which model the main connections between those 6 regions. In order to describe the state of each region, we couple each vertex of the graph with one instance of the model (15). Since each region has its own specificity, we consider that the multiple instances of the model are non-identical, which means that the values of the parameters can differ from one region to another.

Fig. 3

Six regions of Portugal and some of their connections: (1) Norte; (2) Centro; (3) Lisboa e Vale do Tejo; (4) Alentejo; (5) Algarve; (6) Pinhal Litoral.

We introduce the following notations: where H is the matrix of coupling strengths defined by with non negative coefficients , , , and .

Next we define a matrix L of connectivity as follows. For each edge , , we set . If , , we set . The diagonal coefficients satisfy thus L is a matrix whose sum of coefficients in each column is null. For instance, the connectivity matrix of the graph corresponding to Fig. 3 is given by with In this complex network model, we consider that an edge , , models a connection between two regions k and j, which corresponds to human displacements from region k towards region j. Moreover, the parameter models the rate of susceptible individuals in region k which migrate towards vertex j. The parameters , , and are defined analogously. This implies that our model can take into account the situation where a part of the population is not concerned with the migrations. For instance, it is relevant to consider , while and . The set of edges and the coupling strengths stored in the matrix H define what is usually called the topology of the complex network.

Next, we make explicit the equations that describe the state of region : where the time dependence is omitted, in order to lighten the notations. In particular, the parameters can be determined by piecewise constant functions as in model (15).

In the next theorem, we establish the existence and uniqueness of global solutions to the complex network problem (17) as for models (1) and (15). Following [7], we introduce the minimum mortality rate defined by the positive coefficient defined by and the compact region

For any , the Cauchy problem given by (17) and admits a unique solution denoted by , defined on , whose components are non-negative. Furthermore, the region Θ defined by (18) is positively invariant.

The existence and uniqueness of local in time solutions is immediate. The non-negativity property is guaranteed by the quasi-positivity of the non-linear operator determined by the right hand side of system (17). The total population in the complex network, defined by satisfies since the matrix of connectivity L is a zero column sum matrix. It follows that which leads to the intended conclusion. □

In the next section, we apply previous complex network to the real data from COVID-19 in Portugal, since the first confirmed active COVID-19 case, in March 2, 2020, until September 17, 2020.

Portugal case study: complex network with 6 regions

In this section, we calibrate the number of active infected individuals I given by the model with piecewise constant parameters (15), to six distinct regions from Portugal. We show that the model allows to fit well the Portuguese real data available in [10]. Taking into account these parameter values, we perform numerical simulations of the complex network problem (17), where the main goal is to investigate the effect of the topology on the dynamics of the epidemics, and determine a topology that minimizes the level of active infected individuals.

Model with piecewise constant parameters: fitting to real data from 6 Portuguese regions

We first consider the model with piecewise constant parameters (15) and show that the class I fits well the active confirmed cases of infected individuals provided by the Portuguese National Authorities [10] for six regions of Portugal mainland, as depicted in Fig. 3, where (1) represents the region Norte, (2) Centro, (3) Lisboa e Vale do Tejo, (4) Alentejo, (5) Algarve and (6) Pinhal Litoral.

The initial conditions , , , and at time , where corresponds to the day of the first confirmed cases of COVID-19 disease in each region, which differ from one region to another, see Table 1 . Namely, the first cases occurred in each region at: Norte, March 2, 2020; Centro, March 3, 2020; Lisboa e Vale do Tejo, March 3, 2020; Alentejo, March 18, 2020; Algarve, March 8, 2020.

Table 1

Initial conditions for each of the six regions in Portugal.

Region	S0	A0	I0	R0	P0
Norte	3125804	20.15	2	0	0
Centro	1480664	10.15	1	0	0
Lisboa e Vale do Tejo	3659871	10.15	1	0	0
Alentejo	166726	20.15	2	0	0
Algarve	451006	10.15	1	0	0
Pinhal Litoral	271078	30.15	3	0	0

Some of the parameters take constant values for all the period of time under consideration and are the same for the six regions, see Table 2 . Others, even taking constant values for all time window, differ from one region to another, see Table 3 .

Table 2

Constant parameter values that take the same values for the six regions, in Portugal.

Parameter	Description	Value	Reference
Λ	Recruitment rate	0.19%×N_0365	[21]
μ	Natural death rate	181×365	[21]
θ	Modification parameter	1	[23]
v	Transfer rate from A to I	1	[23]
q	Fraction of A individuals that are confirmed to be infected	0.15	[23]

Table 3

Constant parameter values that differ from one region to another, in Portugal.

Region	ϕ (transfer rate from S to P)	δ (transfer rate from I to R)	w (transfer rate from P to S)
Norte	ϕ=112	δ=127	w=145
Centro	ϕ=111	δ=127	w=145
Lisboa e Vale do Tejo	ϕ=111	δ=127	w=145
Alentejo	ϕ = 1	δ=121	w=141
Algarve	ϕ=16	δ=121	w=145
Pinhal Litoral	ϕ = 1	δ=121	w=143

The parameters β, m and p are piecewise continuous, taking different values depending on the considered time interval and the region, see Table 4, Table 5 .

Table 4

Piecewise constant parameter values for the regions Norte, Centro and Lisboa e Vale do Tejo.

Region	β (transmission rate)	p (transfer fraction from S to P)	m (transfer fraction from P to S)
Norte
Time interval: [1;75]	β1 = 1.40	m1 = 0.09	p1 = 0.675
Time interval: [75;122]	β2 = 0.15	m2 = 0.15	p2 = 0.60
Time interval: [122;170]	β3 = 1.28	m3 = 0.14	p3 = 0.56
Time interval: [170;200]	β4 = 1.46	m4 = 0.17	p4 = 0.52
Centro
Time interval: [1;69]	β1 = 1.351	m1 = 0.10	p1 = 0.675
Time interval: [69;94]	β2 = 0.45	m2 = 0.09	p2 = 0.60
Time interval: [94;164]	β3 = 1.10	m3 = 0.16	p3 = 0.56
Time interval: [164;199]	β4 = 1.33	m4 = 0.17	p4 = 0.54
Lisboa e Vale do Tejo
Time interval: [1;36]	β1 = 1.45	m1 = 0.16	p1 = 0.675
Time interval: [36;54]	β2 = 1.13	m2 = 0.13	p2 = 0.69
Time interval: [54;129]	β3 = 1.41	m3 = 0.13	p3 = 0.56
Time interval: [129;160]	β4 = 1.10	m4 = 0.11	p4 = 0.57
Time interval: [160;199]	β5 = 1.43	m5 = 0.17	p5 = 0.54

Table 5

Piecewise constant parameter values for the regions Alentejo, Algarve and Pinhal Litoral.

Region	β (transmission rate)	p (transfer fraction from S to P)	m (transfer fraction from P to S)
Alentejo
Time interval: [1;11]	β1 = 3.50	m1 = 0.40	p1 = 0.675
Time interval: [11;29]	β2 = 7.80	m2 = 0.45	p2 = 0.45
Time interval: [29;84]	β3 = 4.40	m3 = 0.38	p3 = 0.55
Time interval: [84;104]	β4 = 6.90	m4 = 0.46	p4 = 0.48
Time interval: [104;124]	β5 = 2.70	m5 = 0.30	p5 = 0.56
Time interval: [124;184]	β6 = 4.10	m6 = 0.36	p6 = 0.44
Algarve
Time interval: [1;36]	β1 = 1.11	m1 = 0.10	p1 = 0.60
Time interval: [36;69]	β2 = 1.37	m2 = 0.20	p2 = 0.52
Time interval: [69;91]	β3 = 0.80	m3 = 0.20	p3 = 0.55
Time interval: [91;109]	β4 = 1.90	m4 = 0.32	p4 = 0.45
Time interval: [109;149]	β5 = 1.35	m5 = 0.22	p5 = 0.58
Time interval: [149;194]	β6 = 1.35	m6 = 0.30	p6 = 0.55
Pinhal Litoral
Time interval: [1;7]	β1 = 3.00	m1 = 0.09	p1 = 0.675
Time interval: [7;23]	β2 = 8.15	m2 = 0.25	p2 = 0.44
Time interval: [23;77]	β3 = 5.00	m3 = 0.16	p3 = 0.55
Time interval: [77;115]	β4 = 7.80	m4 = 0.25	p4 = 0.42
Time interval: [115;145]	β5 = 4.15	m5 = 0.15	p5 = 0.52
Time interval: [145;186]	β6 = 7.75	m6 = 0.27	p6 = 0.42

Fig. 4, Fig. 5, Fig. 6 show that the model (15) with piecewise constant parameters describes well the active infected cases in the six Portuguese regions under consideration.

Fig. 4

Active infected cases of COVID-19. The black dotted lines represent the real data and the continuous lines represent the active infected individuals, I, output of the SAIRP model with piecewise constant parameters from Table 4. Left: region Norte. Right: region Centro. (For interpretation of the colors in the figures, the reader is referred to the web version of this article.)

Fig. 5

Active infected cases of COVID-19. The black dotted lines represent the real data and the continuous lines represent the active infected individuals, I, output of the SAIRP model with piecewise constant parameters from Table 4, Table 5. Left: region Lisboa e Vale do Tejo. Right: region Alentejo.

Fig. 6

Active infected cases of COVID-19. The black dotted lines represent the real data and the continuous lines represent the active infected individuals, I, output of the SAIRP model with piecewise constant parameters from Table 5. Left: region Algarve. Right: region Pinhal Litoral.

Moreover, the active infected cases of all Portugal can also be modeled, see Fig. 7 .

Fig. 7

Active infected cases of COVID-19 in Portugal. The black dotted lines represent the real data and the continuous lines represent the active infected individuals, I, output of the SAIRP model with piecewise constant parameters. Initial conditions: S0 = 10295894, I0 = 2, , R0 = 0, Q0 = 0. The constant parameters take the values ϕ = 1, , and the remaining constant parameters take the values in Table 2. The first infected case with COVID-19 in Portugal, t = 1, occurred on March 2, 2020. Piecewise constant parameters: t ∈ [1;73], β1 = 1.505;m1 = 0.09;p1 = 0.675; t ∈ [73;90], β2 = 0.50;m2 = 0.09;p2 = 0.65; t ∈ [90;130], β3 = 1.15;m3 = 0.18;p3 = 0.58; t ∈ [130;163], β4 = 0.96;m4 = 0.16;p4 = 0.61; t ∈ [163;200], β5 = 1.50;m5 = 0.17;p5 = 0.58.

Complex network model: numerical simulations for COVID-19 in Portugal

Considering the six regions from Fig. 3, and the parameter values from Section 6.1, we perform numerical simulations of the complex network problem (17). The main goal is to investigate the effect of the topology on the dynamics of the epidemics. In particular, we analyze the existence of a topology which minimizes the average number of active infected individuals during a fixed time interval. Moreover, we analyze if other topologies are likely to worsen the level of infection.

In order to model the mobilities of susceptible and asymptomatic individuals, we set , whereas we fix . We test a sample of 1000 randomly generated topologies among 216 topologies (id est sets of edges), for . For each randomly generated topology, we perform a numerical integration of the complex network problem (17), with the same parameters as considered in Section 6.1. The computation has been performed with Python 3.5 language, in a Debian/Gnu-Linux environment. For each integration, we have computed the number of average number of active infected individuals, that is, individuals of the class I, per day. The results are depicted in Fig. 8 , where the average number of active infected individuals, per day, is depicted for the randomly generated topologies and , and in Fig. 9 , where the average number of active infected individuals is depicted with respect to the coupling strengths for a couple of remarkable topologies. The numerical results reveal the existence of a certain number of topologies that decrease the level of infection, compared to the empty topology, which corresponds to the situation where individuals do not migrate from one region to another. Among those topologies, one minimizes the level of cumulated infected individuals. This optimal topology is depicted in Fig. 10 (b). It is worth noting that this optimal topology connects the regions Norte, Centro, Lisboa e Vale do Tejo and Pinhal Litoral, but it does not connect those four regions to Alentejo and Algarve. In parallel, the numerical simulations exhibit numerous topologies that increase the level of infection. We have presented in Fig. 10 (c) and (d) two examples of such topologies. The corresponding levels of infection are presented in Fig. 9. Topology (c) leads to a level of infection that overcomes the level of the empty topology for only a weak coupling strength, whereas topology (d) seems to permanently overcome the level of the empty topology. We remark that the topologies that increase the level of infection seem to connect the 6 regions of Portugal to the region Alentejo. On the other hand, it can be helpful to authorize parsimonious migrations between regions that have a satisfying control of the epidemic.

Fig. 8

Average number of active infected individuals, per day, for a sample of 1000 randomly generated topologies, for σ = σ = 0.01. The black line shows the level of infection for the empty topology. The green circle shows the optimum topology which minimizes the level of infection, whereas the red circle shows the topology which leads to the highest level of infection.

Fig. 9

Average number of active infected individuals, for each topology, with respect to the coupling strengths σ = σ. Left: optimum topology which minimizes the level of infection of the epidemics. Right: two examples of topologies that can increase the level of infection, compared to the empty topology, which corresponds to the situation where individuals do not migrate from one region to another. Topology (c) leads to a level of infection that overcomes the level of the empty topology for only a weak coupling strength, whereas topology (d) seems to permanently overcome the level of the empty topology.

Fig. 10

Four remarkable topologies. (a) Empty topology, which corresponds to the situation where individuals do not migrate from one region to another. (b) Topology that minimizes the level of infection. (c) Topology that leads to a level of infection greater than the level of the empty topology for only a weak coupling strength. (d) Topology that permanently overcomes the level of the empty topology.

Discussion

One of the factors associated with the spread of the pandemic is the mobility of populations. It was based on this evidence that, in the first phase of the pandemic response, a generalized containment of the population was implemented in a big number of countries. As a result of this restriction in the mobility, it was possible to obtain a decrease in the number of infected. In the numerical simulations considered in this work, for the complex network problem (17), it is assumed that only susceptible (S) and asymptomatic infected (A) are circulating among the 6 identified regions of mainland Portugal. Starting from the number of active infected people (I) for the empty topology, in which it is assumed that there are no flows between regions, the results of the simulations indicate that there are more topologies that increase the number of active infected people than those that cause this number to decrease. Hence, this result reinforces the need to reduce circulation between regions as one of the tools to reduce the spread of the virus and decrease the number of active infected people. Still, it is interesting to note that there are other topologies (around 1 out of 3) in which the number of active infected individuals decreases with the movement of people between regions. This finding can be seen in Fig. 8, where the average number of active infected individuals for the 1000 tested topologies is shown. According to the topology that minimizes the number of active infected individuals, whose representation is shown in Fig. 10 b), it is possible to verify that the strategy of preventing the flow to the regions of Alentejo (4) and Algarve (5) contributes to reduce the number of active infected. This scenario was to be expected as those regions have a much lower number of cases than the rest of Portugal. Bidirectional circulation between Norte (1) and Centro (2) regions leads to a reduction in the number of active infected people (if one compares with the empty topology), but the movement between Pinhal Litoral region (6) and Lisboa e Vale do Tejo region (3) should only be done in the direction to prevent circulation from the region where there are less infected (6) to the region (3) that has a high number of active cases. In Fig. 10 c) it is possible to verify that the existence of this connection in the direction (6) to (3) increases the number of active cases. The topologies involving the regions of Alentejo (4) and Algarve (5) increase the number of cases. One extreme example of such topology can be found in Fig. 10 d). In this case, the direction of the flows presents a topology that converges to the Alentejo (4) and it is expected that this will translate into an increase in the number of infected (even in the case where the connection strength is small, as shown in Fig. 10 c) and Fig. 9 at right), as it is the region where the number of cases is lesser. The way the epidemic evolves is different between the different regions of a country, but also between countries. These specific characteristics of context (cultural, demographic, economic, …) seen in the number of active infected individuals (which is not uniform across regions), can be modeled by pseudo-periodic piecewise functions, as proved in this work. Accordingly, the application of such methodology to identify flows between regions or countries is a tool with enormous potential in the current pandemic context, and can be applied in the management of outbreaks (in regional terms) but also to manage the opening/closing of borders.

Conclusion and future work

The way the epidemic evolves is different between the different regions of a country, but also between countries. These specific characteristics of context (cultural, demographic, economic, …) seen in the number of active infected individuals, which is not uniform across regions, can be modeled by pseudo-periodic piecewise functions, as proved in this work. Accordingly, the application of such methodology to identify flows between regions or countries is a tool with enormous potential in the current pandemic context, and can be applied in the management of outbreaks (in regional terms) but also to manage the opening/closing of borders. In the implemented simulations, the intensities of the flows ( and ) are assumed to be equal. Still, and as future work, it would be useful, in terms of epidemiological management, to be able to model the intensity of flows by class (Asymptomatic, Susceptible, Active Infected, Recovered, and Protected) between regions.