Modeling nonlocal behavior in epidemics via a reaction-diffusion system incorporating population movement along a network.

The outbreak of COVID-19, beginning in 2019 and continuing through the time of writing, has led to renewed interest in the mathematical modeling of infectious disease. Recent works have focused on partial differential equation (PDE) models, particularly reaction-diffusion models, able to describe the progression of an epidemic in both space and time. These studies have shown generally promising results in describing and predicting COVID-19 progression. However, people often travel long distances in short periods of time, leading to nonlocal transmission of the disease. Such contagion dynamics are not well-represented by diffusion alone. In contrast, ordinary differential equation (ODE) models may easily account for this behavior by considering disparate regions as nodes in a network, with the edges defining nonlocal transmission. In this work, we attempt to combine these modeling paradigms via the introduction of a network structure within a reaction-diffusion PDE system. This is achieved through the definition of a population-transfer operator, which couples disjoint and potentially distant geographic regions, facilitating nonlocal population movement between them. We provide analytical results demonstrating that this operator does not disrupt the physical consistency or mathematical well-posedness of the system, and verify these results through numerical experiments. We then use this technique to simulate the COVID-19 epidemic in the Brazilian region of Rio de Janeiro, showcasing its ability to capture important nonlocal behaviors, while maintaining the advantages of a reaction-diffusion model for describing local dynamics.

Introduction

The outbreak of COVID-19, which started in 2019 and is still continuing, has caused unprecedented disruption in terms of both human cost and economic damage. In order to better understand the dynamics of the disease spread, with the hope of ultimately improving policy and public health outcomes, there has been an explosion in the study of mathematical modeling of infectious disease. These models have taken many forms, and a comprehensive survey of the literature is beyond the scope of the current work. However, we note that many different approaches have been used to model the epidemic, including machine-learning and data-driven approaches [1], [2], [3], [4], [5], models using a classical compartmental approach, together with parameter estimation techniques [6], [7], delay differential equations [8], partial differential equations [9], [10], [11], [12], network-based methods [9], [13], [14], [15], as well as agent-based [16], [17] and multiscale models [18]. We note that the cited works represent just a small sample of the total literature, and further that the various approaches discussed are not mutually exclusive. For a recent survey, see [14].

The basis for most mathematical models of infectious disease is rooted in the compartmental theory of Kermack and McKendrick [19] (see [20] for a more modern treatment). This approach divides a given population into different states corresponding to disease status, with the underlying equations governing the evolution of the population composition in terms of disease-state. This mechanism forms the backbone behind the classic S–I–R (susceptible–infected–removed) model and its variants: for instance, the S–I–S (susceptible–infected–susceptible) and S–E–I–R (susceptible–exposed–infected–removed), among others. While the initial model presented by Kermack and McKendrick is quite general, and presented in the form of an integral renewal equation, under certain assumptions, namely a mass-action transmission mechanism and exponentially-distributed sojourn times, we obtain the familiar systems of ordinary differential equations (ODE) widely used today. These systems should be recognized, however, as special cases within a more general framework [19], [20]. We note that many models, including agent-based and evolving network models, even if not expressed as integral or differential equations, still generally utilize some sort of compartmental structure [16], [21].

Such compartment models have historically been applied with great success in the modeling of infectious disease, and offer a clear mechanistic description of epidemic progression. However, in their most widely-utilized ODE formulations, relevant factors may not be taken into account. Among the most important such factors, and the subject of this work, is the introduction of spatial information into the epidemic model. Many approaches have been proposed to describe epidemic evolution in space as well as time. A common approach is to further stratify the compartmental structure by geographic location, such that different compartments also correspond to different geographic areas, often delimitated by political jurisdiction. For COVID-19, this approach has been utilized in e.g. [5], [6], [13], [15] and offers the advantage of computational efficiency and relative ease of implementation. Additionally, the spatial evolution is not inherently local, and nonlocal processes (ie, transmission between individuals in distant regions) can be easily incorporated. However, a continuous description of spatial dynamics with this approach is impossible, and as the spatial resolution increases, the number of necessary compartments can quickly become intractable.

In contrast, partial differential equations (PDE) incorporating the familiar compartmental structure have been offered as an alternative. Most commonly, models of this type incorporate evolution in space via a reaction–diffusion equation [22], [23]. Works incorporating PDE can be seen in [3], [9], [10], [11], [24], [25], [26]. Most of these used a reaction–diffusion process, with [9] introducing a multiscale approach utilizing convection for nonlocal processes and diffusion for local processes. In [27] convection is also employed to account for population mobility.

PDE models offer the large advantage of a true, spatially continuous description; however, they too have downsides. From the practical point of view, they are generally more difficult to implement and require greater computational resources than ODE models. From the modeling point of view, the ability of a linear diffusion model to accurately describe the spatial evolution of an epidemic in human populations, given the nature of human movement, is questionable. While many of the COVID-19 related applications of these models mitigate this problem somewhat through, for instance, nonlinear diffusion and depensation effects [3], [10], [11], [24], [25], [26], the model still behaves locally. Although [9] does make an excellent attempt at a nonlocal description, ultimately, despite the presence of multiple scales of spatial evolution and convection, the underlying process is still local, as contagion is still spread along with the intermediate points between origin and destination. It is well-known that bilaplacian [22] or fractional diffusion [28] operators can describe nonlocal dynamics. Some works have indeed applied such techniques to COVID-19 (see [29], [30]); however, the complexity of these models in terms of computational effort and stability, as well as their still-developing theory, makes their widespread adoption, at least at the current time, difficult. Further, the nonlocal movement in human populations is far from random; indeed, it often follows a predictable pattern [31].

In contrast to these approaches, we propose herein a combination of the reaction–diffusion model first proposed in [11] with network-based methods, akin to those proposed in [9], [13], [14], [15]. Similarly to the modeling paradigm introduced in [9], we employ a PDE model and postulate that diffusion occurs at the local (i.e., within-city) level. However, rather than defining nonlocal dynamics via convection, we instead model such dynamics as sources and sinks within a network. We embed a spatial network structure within the underlying PDE system, and define a network via transport operators between different areas in the spatial domain (i.e., inter-city mobility). These operators then allow portions of the population to, for instance, move from one city to another without diffusing the disease along the way, as is typically the case with vehicular transport. We believe that such a method offers a reasonable, easily implemented, stable approach that maintains the advantage of a PDE-based model for local dynamics while also incorporating important nonlocal effects.

We outline the article as follows: we first introduce and provide a brief explanation of the underlying equations. We then discuss the idea behind the population transfer network and its mathematical formulation, as well as important details regarding its numerical implementation. We then provide numerical examples on an idealized problem, as well as a realistic problem in the Brazilian state of Rio de Janeiro, to provide evidence of the model’s ability to describe the relevant local and nonlocal dynamics found in human epidemics.

Governing equations

The governing equations are based on the spatio-temporal SEIRD model, presented in [11], [24], [32]. We postulate that only the susceptible, exposed and recovered compartments are able to move along the network. The system of equations reads: where , , , , and denote the densities of the susceptible, exposed, infected, recovered, and deceased populations, respectively. and denote the contact rates between infected and susceptible individuals and exposed and susceptible individuals, respectively (units days−1), denotes the incubation period (units days−1), corresponds to the exposed recovery rate (units days−1), the infected recovery rate (units days−1), represents the mortality rate (units days−1), and , , , are the diffusion parameters of the different population groups as denoted by the sub-scripted letters (units km persons−1 days−1). We additionally denote the sum of the living population as ; i.e., . Finally, , and are population transfer operators, defining the movement of susceptible, exposed and recovered sub-populations within the population transfer network as sources and sinks in the reaction–diffusion system. The precise definition of these operators will be discussed in detail in the following section. The system of Eqs. (1)–(5) must be supplemented with proper boundary and initial conditions; we assume hereafter no-flux boundary conditions.

The movement occurs between paired regions, ie, nodes that share an edge in the underlying network. Movement away from one region results in movement into the other region. In this sense, in a given time, when persons move away from a region, the population transfer operator acts as a sink term; likewise, when persons move into a region, the population transfer operator acts as a source term.

For the numerical solution of (1)–(5), we discretize in space using a Galerkin finite element formulation. We apply the second-order Backward Differentiation Formula (BDF2) in all cases, which offers second-order accuracy while remaining unconditionally stable. All simulations have been performed using the libMesh, a C＋＋ FEM open-source software library for parallel adaptive finite element applications [33].

Defining the transfer operator

In this section we define the population transfer operators within a network. For simplicity, we consider a single compartment and a network consisting of two nodes; extending the definitions to larger networks and other compartments follows immediately.

Model problem

We consider Eq. (1), defined in a domain , with two distinct subsets , such that .

Transfer network definition

We define the directed, weighted, and without self-loops transfer network as , where and represent, respectively, the set of nodes and edges between nodes and . For the time being, let consist of only two nodes, i.e., . We also assume that is a time-varying graph (TVG) [34], where the edges depend on time according to the following assumptions (note this dependence on time is assumed if not explicitly denoted):

Every units of time, a percentage of the population in move to ;

Every units of time, a percentage of the population in move to .

This is equivalent to stating that the probability that an individual moves from to (or vice versa) over time units is .

For example, if one assumes that this probability is distributed exponentially with parameter , we find that: which after some straightforward algebra gives:

Transfer operator definition

We now define the population transfer operator along the network . The dependence of on the network is assumed and not denoted explicitly hereafter. Such an operator is not unique, and in this section we outline the properties that an operator of this type must satisfy. We split the operator into two components, a donor operator which transfers a population from to at a rate and the corresponding receiver operator which receives a population from into at a rate . is then given as the sum of these components: . (and its individual components) is a linear operator in , although it may have additional dependencies on , any of the compartments, and/or in general, and need not be linear in these variables.

must satisfy the following conditions:

Non-negativity of the donor operator: for all ;

Population (mass)-conservation: for all .

Both conditions are physically motivated; the first ensures that the quantity of persons transferred away from given node must not exceed the total persons present in . Additionally, this condition is necessary to guarantee the coercivity (and thereby well-posedness) of the variational problem corresponding to (1)–(5) (see e.g. [35]).

The population conservation condition is necessary to guarantee that the transfer operator only models movement among the different regions, and does not change the overall quantity of persons (or mass) in the system.

We now discuss the donor and receiver operators separately.

Donor operator

In the movement of populations from to (and vice versa) at the rate (and ), the donor operator defines the removal of persons from the origin node.

We propose the following, simple definition of : where are characteristic functions: Assuming the , it follows immediately from that , ensuring that the non-negativity condition is satisfied.

Receiver operator

The receiver operator receives the populations taken from the donor region by the donor operator and distributes them in the receiver region, and has the following form: The kernel functions are defined such that: and

The transfer function defined by (9) , (11) satisfies the mass-conservation property.

One may verify this with straightforward computation:

Some examples of transfer operators

Transfer operators of the type introduced above are non-unique, and any operator satisfying the necessary properties may be applied in practice. Besides the clearly problem-dependent , the definition of the kernel functions , which corresponds to the distribution of the transferred population within the receiver region, is the most important component in constructing .

As shown in the preceding sections, as long as the conditions (12)–(13) are satisfied, then one may define an acceptable transfer function. However, these conditions are quite non-restrictive, and many such functions will satisfy them. Defining is therefore far from straightforward, and problem-specific considerations, computational constraints, and available data are likely to inform such choices in practice. Below we list two possible, obvious definitions.

Uniform distribution. A possible definition is given by: This defines the simplest possible transfer function, which distributes the population uniformly among the receiver domain. It is computationally simple and, therefore, the approach applied in the simulations of this work. However, it is questionable how well this definition captures reality in most instances.

Receiver-proportional distribution. Another possible, natural choice is: which distributes the population in proportion to the existing population density in . This is likely more realistic than the uniform distribution; however, it requires the additional evaluation of an integral, must be updated in time, and introduces a nonlinearity into the problem.

Depending on how one chooses to discretize the term (16) in time, it may become necessary to apply some nonlinear iteration scheme at each time step. We note, however, that such schemes are generally necessary when solving (1)–(5) anyway, and this may not necessarily represent a large additional computing cost.

Results

Simple test problem

In order to verify the consistency of the proposed approach, we first consider a rectangular domain divided in two regions. The dimensions of each region are for Region 1 and for Region 2. We start with 1000 people in each region, i.e., Region 1 has 1000 and Region 2 has 500 (Fig. 1).

Fig. 1

Simple test problem schematic.

The objective is to seasonally send people from Region 1 to Region 2 and people from Region 2 to Region 1. We define a weekly pattern in which people leave Region 1 on Fridays and go back home on Sundays. For instance, assuming day 1 is Monday, people go from Region 2 to Region 1 on day 5 (Friday), and on day 7 (Sunday), people go from Region 1 to Region 2. We do not consider movement on day 0.

We consider a time-step of 1/4 days, and define the population movement as 5% of the population in a given region per time-step. Hence: We define as in (15).

We use four different meshes to study population (mass) conservation, each discretized with triangles. Three of the meshes are structured with different element mesh sizes (25100, and 100 × 200). The fourth mesh is unstructured, with 18,231 elements of sizes varying between the largest and smallest elements of the structured meshes. In Fig. 2 we show the coarsest structured and the unstructured meshes.

Fig. 2

Simple test problem meshes. (A) 25 × 50 elements mesh. (B) unstructured mesh.

As we are not considering diffusion in this example, the only change in population distribution is due to . Therefore, we may promptly obtain an analytical solution.

In Fig. 3, Fig. 4, we plot the total populations in each time-step for Region 1 and 2, for each mesh. Fig. 5, Fig. 6 show, for each mesh, the relative error in time between the simulated and analytically-computed number of people in Regions 1 and 2. Fig. 7 shows the total population for each mesh. Note that the total population here considers all compartments, including the deceased compartment. Since there is no movement across the boundaries, we do not consider new births, from the mass-conservative property of the transfer operator, the total population (including deceased) should remain the same at all time steps, up to numerical error. We stress that the total population in this context should not be confused with the living population in Eqs. (1)–(4). Finally, Fig. 8 shows the total population (mass) relative error between the simulations and analytical solution.

Fig. 3

Population at Region 1. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4

Population at Region 2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 5

Population relative error at Region 1. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 6

Population relative error at Region 2.

Fig. 7

Total population.

Fig. 8

Total population relative error for each grid. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

We briefly note that, as we keep all of the percentages fixed at 5% per time-step, we observe population drift over time, with different numbers of people moved at each instance. This occurs due to the fact that, after sending a given percentage of the population in Region 1, when we send the same percentage of the population in Region 2 back to Region 1, this population now includes the originally transferred persons. Noting that 5% of 1000 is 50, we may avoid this problem, and ensure a constant number of persons transferred, by defining the transfer rates instead as: While defining in such a manner requires integral evaluations, we note that these are already computed in order to define (11).

Depending on the mesh refinement, we start the simulation with a different amount of total population, close to 2000. These discrepancies occur due to errors in the numerical representation of the initial condition in Region 1 (Fig. 5). This occurs due because the nodes along interface between the two regions must be defined as belonging to only one region. Therefore, when setting the initial condition, we incur small numerical error near the interface, as it can only appear in one region. As our numerical experiments show, refining the mesh reduces this problem. We also note, however, that this problem is more pronounced for the unstructured mesh; this suggests that using a mesh that is conformal to the interface boundaries may also help to alleviate this effect.

Rio de Janeiro state - Brazil

In [26] the spread of COVID-19 in the state of Rio de Janeiro, Brazil was simulated. However, difficulties in reproducing the virus spread dynamics in some counties, particularly in Cabo Frio (Region 2) and Campos dos Goytacazes (Region 3) (see Fig. 9), were encountered. In that simulation, there were no initial infected or exposed populations in these counties (Fig. 10), making diffusion the only possible way for the virus to reach those areas. Given the large distance between these areas and the hotspots near the city center of Rio de Janeiro metropolitan area (state capital, Region 1, respectively 157/279 km distance from Regions 2 and 3), as well as the large sparsely-populated areas between them (see Fig. 11), population-weighted diffusion was not able to properly represent these dynamics. Clearly, nonlocal effects played a significant role in the spread of COVID-19 in these areas, motivating the application of the introduced nonlocal population-transfer operator.

Fig. 9

Sketch of the Rio de Janeiro state, Brazil, highlighting the Rio de Janeiro, Cabo Frio and Campos de Goytacazes counties.

Fig. 10

Initial exposed/infected population at Rio de Janeiro state. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 11

Initial susceptible population at Rio de Janeiro state. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Note that, for all simulations considered in the following, parameters unrelated to the population transfer network are exactly as considered in [26], which are: ,  ,  ,  . These values are based on available data from the literature regarding the mortality, incubation period, and recovery time for infected and exposed patients [24]. The contact rate and the diffusion parameters are time-dependent parameters, with values informed based on estimates of the effective reproduction number and the perception of average weekly social isolation based on the movement of cell phones in the state of Rio de Janeiro. More details on how such measures informed parameter values can be found in [26].

We now describe the definition the population transfer network. In all instances, the number of people moving between each region was based on collected data from the Brazilian Institute of Geography and Statistics (IBGE).1

Cabo Frio (Region 2) is a coastal county, where there are seasonal visitors during the weekends. Therefore, we define as people going from Rio de Janeiro to Cabo Frio on Fridays and coming back on Sunday, as in the previous simple test problem. The main difference is the percentages of people that come and go from each city. Since Rio has 6.7 million inhabitants and Cabo Frio has about 200 thousand, movement of the same amount of persons will correspond to different percentages of the population. We assume that about 30,000 people have this weekly pattern movement, or approximately 15% of the population of Cabo Frio and 0.5% of the population of Rio de Janeiro. We apply the same rate to , and compartments. We note that, due to the previously-discussed population drifting effects, diffusion among the population, as well as the fact that these are applied to both and , they may not correspond to exactly the same number of persons moved at each time period. Given that is several orders of magnitude larger than other compartments, the simplified rates defined above do not create serious issues in general. Indeed, these differences in number and composition of the population transferred in time may in fact be more realistic, as we would not expect such numbers to remain constant. Nevertheless, one may also define the similarly to (18) in order to ensure that population movements remain consistent in time.

Another region that was underestimated in [26] for similar reasons to Cabo Frio is Campos dos Goytacazes (Region 3). Accordingly, we postulate movement between Rio de Janeiro and Campos dos Goytacazes, but consider a different weekly pattern, as Campos dos Goytacazes is not a vacation area. Rather, many people in the oil industry work there for one week at a time on the offshore oil platforms, returning to Rio de Janeiro the next week. We assume that there are about 25,000 people who move in this manner per week, corresponding to 0.375% the population of Rio de Janeiro, and 5% the population of Campos dos Goytacazes. Each week, 25,000 people are exchanged between Campos dos Goytacazes and Rio de Janeiro, corresponding to 50,000 persons alternating weekly shifts. Analogously to (19) we then define this pattern as: for , and .

Lastly, to complete the network, we also consider movement between Campos dos Goytacazes and Cabo Frio. We define this movement analogously to (19); however, rather than moving approximately 30,000 people between the regions, we move about 10,000, corresponding to about 2% of the population of Campos dos Goytacazes and 4.45% of the population of Cabo Frio:

In total, we simulate eight different cases, considering different networks as presented in Table 1 and Fig. 12. The mesh of the state of Rio de Janeiro territory contains 11,632 elements and 6082 nodes, and, although unstructured, the element size is very close to uniform, with a mean area of 1 km and 95% of elements between 0.5 and 1.5 km. Simulations are done considering piecewise linear finite elements.

Table 1

Definition of the movement between Rio de Janeiro city (RJ), Cabo Frio (CF) and Campos dos Goytacazes (CG) for each case.

		Movement between
Case	PDE	RJ and CF	RJ and CG	CF and CG
1	X
2	X	X
3	X		X
4	X			X
5	X	X	X
6	X	X		X
7	X		X	X
8	X	X	X	X

Fig. 12

Definition of the movement for each case. The region in blue is Rio de Janeiro city (RJ), red is Cabo Frio (CF) and green is Campos dos Goytacazes (CG). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

In Fig. 13, Fig. 14, Fig. 15 we show the comparison of cumulative deaths between the observed data and the simulated cases. Note that, differently from the toy model, diffusion is considered in all simulations; hence, we do not have an analytic solution as in the previous example. Nonetheless, the total population in the region should be conserved. In Fig. 16 we show the total population relative error with time, which is substantially under 1% throughout the time period.

Fig. 13

Cabo Frio cumulative deaths. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 14

Campos dos Goytacazes cumulative deaths. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 15

Rio de Janeiro city cumulative deaths. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 16

Total population relative error considering population movement. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

When considering only movement between Cabo Frio and Campos dos Goytacazes (Case 3) we see very little change in the spread of the virus in these counties. As both counties both have no infected or exposed at the beginning of the simulation, no nonlocal propagation of the virus occurs, despite the population transfer between regions. However, when incorporating movement from Rio de Janeiro, the primary hotspot of the epidemic (as in Cases 1 and 2), we see a marked improvement in the simulated agreement with measured data. In particular, we observe that the nonlocal movement defined by the population transfer network enables the epidemic to spread into Cabo Frio and Campos dos Goytacazes. We also note that the deceased population in Rio shows little change from the introduction of nonlocal population movement.

To compare the effects of the population transfer network, Fig. 17 shows the spatial distribution of the deceased population at days for Case 1 (no population transfer) and Case 8 (full network incorporating population transfer between Rio de Janeiro, Cabo Frio, and Campos dos Goytacazes), while Fig. 18 reports the corresponding relative error. In Rio de Janeiro, we see that the purely diffusive model with no nonlocal population transfer (Case 1) already works well, and the introduction of the population transfer network does cause a deterioration in accuracy (Case 3). In contrast, in Cabo Frio and Campos dos Goytacazes, the purely diffusive model fails to recreate the observed dynamics. However, upon the introduction of the population transfer network, in Cabo Frio, the results closely match the observed data. In Campos dos Goytacazes, while a significant discrepancy remains between the simulated and measured data, the results are nonetheless significantly improved, as fatalities are no longer close to zero. The discrepancy may be due the simplistic assumptions about the nature of the population movement; under more realistic assumptions, these results may improve. For instance, the current model does not consider population movement from the bordering state of Espírito Santo.

Fig. 17

Spatial distribution of the deceased population at days for Case 1 (no population movement) and Case 8 (population movement network between RJ, CF and CG). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 18

Relative error between the deceased population at days for Case 1 (no population movement) and Case 8 (population movement network between RJ, CF and CG). Left: negative values. Right: positive values. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Discussion and future directions

In the present work, we have introduced a new approach to modeling nonlocal dynamics in reaction–diffusion models of infectious disease transmission. The proposed method involves defining different subregions in the domain of interest as nodes in a network, with the movement between nodes defined via edges in a time-varying graph. While conceptually simple, defining an operator properly representing these dynamics within a reaction–diffusion PDE is nontrivial. We discussed the proper definition of such an operator that both respects sensible physical constraints and maintains the mathematical well-posedness of the underlying system. In particular, the operator must be non-negative and population (mass)-conservative. We then demonstrated conditions that ensure that the transfer operator satisfies these important properties, and proposed some natural candidates for its definition.

We then performed two numerical experiments with different goals. The first was designed to confirm that our numerical implementation of the network transfer operator demonstrated the necessary conservation and non-negativity properties as outlined in Section 3. We found that our implementation generally respects mass conservation, with numerical errors below one percent for all meshes considered, and tending to zero with mesh refinement. Unsurprisingly, we also found that structured meshes generally demonstrate superior mass conservation properties.

Our second numerical experiment was a proof-of-concept designed to demonstrate the potential of this approach on a realistic problem. We considered a three-node network of three distinct regions in the state of Rio de Janeiro in Brazil. Without considering population transfer along the network, the PDE model is able to accurately capture the dynamics in the city of Rio de Janeiro. However, it is unable to account for transmission in the Cabo Frio and Campos dos Goytacazes regions, due to the inherent lack of nonlocal dynamics in a purely diffusive model. We then showed that, under reasonable assumptions, introducing a network structure endowed with a population transfer operator is able to account for such nonlocal dynamics, improving performance of the model in the Campos de Goytacazes and Cabo Frio regions.

There are several important directions for future research. We considered only the uniform distribution of the population transfer operator in the present work; however, it is important to know the sensitivity of the model to this specific definition. Given the observed mesh dependence on population conservation, it may also be beneficial to explore adaptive meshing algorithms that take the geography of the network regions into account, as this may improve performance. Techniques for non-conformal immersed boundary problems, such as the Finite cell method, [36] may also be helpful for problems in which the boundaries of the subregions in the network are highly irregular. Perhaps most importantly, we stress that the simulations shown here were also preliminary and designed primarily to show the potential of this methodology to account for nonlocal dynamics, and the considered network was very simple. However, applying the described approach to more complex networks with more nodes is important. Similarly, we defined the edges in a rudimentary way. In future work, these edges must be defined more realistically, incorporating more accurate data. Such extensions represent an important step toward the development of a realistic model which may ultimately inform public health decision-making authorities.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.