Mathematical modeling of spatio-temporal population dynamics and application to epidemic spreading.

Agent based models (ABMs) are a useful tool for modeling spatio-temporal population dynamics, where many details can be included in the model description. Their computational cost though is very high and for stochastic ABMs a lot of individual simulations are required to sample quantities of interest. Especially, large numbers of agents render the sampling infeasible. Model reduction to a metapopulation model leads to a significant gain in computational efficiency, while preserving important dynamical properties. Based on a precise mathematical description of spatio-temporal ABMs, we present two different metapopulation approaches (stochastic and piecewise deterministic) and discuss the approximation steps between the different models within this framework. Especially, we show how the stochastic metapopulation model results from a Galerkin projection of the underlying ABM onto a finite-dimensional ansatz space. Finally, we utilize our modeling framework to provide a conceptual model for the spreading of COVID-19 that can be scaled to real-world scenarios.

Introduction

Recently, epidemiological models have received lots of attention due to the COVID-19 pandemics. State of the art models for analyzing the epidemics spreading and doing forecasts include agent-based models (ABMs) [1], [2], [3], [4], metapopulation models [5], [6], [7] and deterministic models based on ordinary differential equations (ODE) [8], [9], [10], [11], [12], [13], [14]. While ABMs have a spatial resolution on a microscale, capturing the spatial movement and interactions of each individual agent, metapopulation models distinguish between subpopulations as groups of agents which are homogeneously mixed. Both of these two models can take into account stochastic effects in the interaction dynamics and the spatial exchange. In contrast, ODE models are purely deterministic and assume homogeneous mixing in the overall space of motion, not involving any mobility dynamics. Although ABMs often reflect the reality better than metapopulation and ODE-based models, they are very costly to simulate and to re-calibrate to new parameter sets, especially for large population numbers that are considered when modeling a pandemic. This motivates the consideration of model reduction techniques which allow to keep the characteristic properties of the detailed ABM, but significantly reduce the computational effort for simulation and analysis. In this work we present a step-wise model reduction which decreases both spatial resolution and stochasticity of an ABM to a level appropriate for fast and accurate simulations of real-world systems.

In the ABM, which we will introduce as a first modeling approach, each individual agent owns a position in some given spatial environment, as well as a status which can, e.g., refer to the agent’s state of health, opinion or level of knowledge [15]. In the course of time, the agents move within the environment and interact with each other, thereby possibly adopting the status of other agents in their spatial vicinity. The propensity for such adoption events to occur is defined by a adoption rate function. Simulations follow the movement and interaction of each individual agent and produce a high computational burden.

In order to approximate these detailed dynamics by the less complex stochastic metapopulation model (SMM), we assume that the spatial environment splits up into areas between which mixing is rare, e.g., areas of metastability where the transition rate to other areas is small. Groups of agents that are located in the same area are termed subpopulations. In application to epidemic dynamics, one can think of these subpopulations as inhabitants of a city or country, but also as smaller groups of people, e.g. pupils of the same school, age classes, working groups or households. Within each of the subpopulations, individuals of the same status are considered as indistinguishable. An agent can anytime interact with all other agents of the same subpopulation, while an interaction with agents of a different subpopulation requires a preceding spatial transition to the associated spatial area. Both the interactions within and the transitions between the subpopulations are treated as stochastic events inducing a jump in the population state, which is defined by the number of agents of each status in the different subpopulations. Analytically, this stochastic metapopulation model results from a Galerkin projection of the underlying ABM based on a spatial partition into metastable subsets [16], [17], [18]. In this work, we will define the corresponding projection operator and derive explicit formulas for the matrix representations of the projected diffusion and interaction operators for both first- and second-order status-changes. We will analyze the approximation quality for a simple example. Simulations of the Markov jump process described by the SMM can be done by means of the event-based stochastic simulation algorithm [19]. This is much more efficient than the agent-based simulations, but the effort still scales linearly with the population size.

Given that the number of agents in each subpopulation is large, a further approximation of the dynamics is possible, replacing the purely stochastic jump dynamics by a piecewise-deterministic Markov processes. This approach, discussed for, e.g., chemical reaction networks [20], has well-established mathematical foundations, going back to Kurtz [21], [22], who studied under which conditions the stochastic dynamics of a large population are well approximated by deterministic processes. In the resulting hybrid modeling approach, the interactions within each subpopulation are described by continuous, deterministic dynamics in form of ODEs, while the comparatively rare exchange between the subpopulations are still treated as stochastic, discrete events. These piecewise-deterministic dynamics can be efficiently simulated by a combination of the stochastic simulation algorithm and an ODE solver [23], [24], and the effort is independent of the population size.

There exist several approaches to describe epidemic (or other types of) population dynamics by semi-deterministic processes [25], [26]. A piecewise-deterministic metapopulation model (PDMM), where deterministic dynamics within the subpopulations are combined with stochastic transitions between the subpopulations, is formulated in [27] for modeling a cattle trade network. We will here use a similar hybrid approach in order to model and analyze epidemic spreading kinetics. This efficient PDMM-approach will allow to easily simulate and compare several scenarios induced by various measures that are taken to control the spreading of the COVID-19 within a population. We will investigate the interrelation of the kinetics within the different subpopulations and analyze the critical transition time of the virus spreading between separate subpopulations depending on several choices of containment measures.

Outline.

We first introduce the mathematical formulations for each of the three modeling approaches and describe the details of the coarse-graining steps in Section 2, see Fig. 1 for an overview. This is illustrated by a guiding example for which we also analyze the approximation quality and highlight the gain in computational efficiency. Moreover, we derive explicit results for the relation between the microscopic adoption rate functions, defining the propensity for an agent to change status in the ABM, and the macroscopic adoption rate functions as an analogue in the SMM, see Section 2.2. Finally, we apply the PDMM framework in Section 3 to model the spreading of the COVID-19 virus.

Fig. 1

Illustration of hierarchy of modeling approaches considered herein and the related steps of the coarse-graining procedure. Dots refer to agents that have a status (blue or red), boxes on the mesoscale represent subpopulations. Straight lines between dots indicate connections between agents that enable interaction. Gray lines on the microscale show movement trajectories of the agents, arrows and lines between boxes on the mesoscale represent possible transitions between subpopulations. For the piecewise-deterministic metapopulation model, the population state is depicted by the proportion of blue and red.

Hierarchy of modeling approaches

In this section, we present three basic modeling concepts for interacting agent dynamics and show how they are related to each other. We start in Section 2.1 with the most detailed agent-based model, which tracks the movement and interactions of each individual agent present in the system. Assuming that within certain regions of the environment the spatial mixing is fast compared to the interaction dynamics, the dynamics may be approximated by a stochastic metapopulation model which is introduced in Section 2.2. A partial approximation of the stochastic dynamics by deterministic dynamics, which becomes reasonable in cases of large populations, leads to a piecewise-deterministic metapopulation process which is considered in Section 2.3.

Agent-based model

We consider a set of agents that are moving in a compact domain as independent realizations1 of a diffusion process on , where is a continuous time interval, e.g. . In addition to the position , each agent also has a status . Let denote the vector of all agents positions, while is the vector of all agents’ status and is the system state. We define the space of all possible system states as .

The dynamics in the agents’ status are modeled by a continuous-time Markov jump process. Each agent can change its status by transitions of the form for , called status-change/adoption events. We will distinguish between (i) first-order adoptions given by status transitions that require no interactions between agents and (ii) second-order adoptions based on pairwise interactions with close neighbors. The rate for agent to undergo such a transition may in general depend on the whole system state (i.e. the positions and status of all other agents),2 and is given by the adoption rate function Given that the agents are positioned according to and have a status according to , it is the probability per unit of time for agent to switch from status to status . Note that the adoption rate function does not explicitly depend on , i.e., it is not the case that every agent has its own propensity function. Dependence on is only indirect, with the propensity depending on the status and position of the agent. More concretely, as for first-order adoptions, we consider rate functions of the form where denotes the indicator function of a status, i.e.  if and otherwise, and gives the rate for a status transition from to depending on the spatial location of the acting agent. In particular, we set for all positions .

Equivalently, for second-order adoptions, we set where is a function defining the rate for a status adoption from to depending on the positions of two interacting agents. In a more special setting, this rate depends only on the distance between the interacting agents, e.g., they need to be closer than some interaction radius to interact, as in the Doi-model [28] in the context of chemical reaction systems. The underlying idea is that interactions of agents (which may induce adoptions of the status from other agents) require proximity of the agents in physical space. For this case, we set for a constant , where for is the distance indicator function: For there are no status transitions and we set .

The coupling of the diffusion process and the jump dynamics given by status-changes leads to a Markov process on the system state space . Let be the probability mass function for the process to be in the system state at time , where the marginal with respect to is a continuous density function.

We define for each status an operator that describes the change of the probability mass function through the motion of a single agent under the condition that he is currently in status . Then we can write down an operator for the movement of the agents, where is defined as for acting on a function with respect to the component of (see Example 1 for details). Note that, consequently, acts only on the position part of the probability mass function in (5).

As for the adoption dynamics, we define where denotes the th unit vector of .3 The first term on the right-hand side refers to the outflow from the current state through adoption events and the second term to the inflow through adoption events that would lead to the current state. The change of including movement and status transitions is then given by the set of differential equations

The strength of this general agent-based approach is that different types of restrictions regarding the interaction dynamics can be included and quite complicated dynamics can be formulated. On the other hand, the coupled differential equations which describe the system are usually not analytically solvable. Instead, Monte Carlo (MC) simulations of the dynamics are required to sample the quantities we are interested in, and the simulation of such complex systems can be numerically very costly, especially if local neighborhoods have to be computed in every time step. We realize the simulations with an algorithm that is a combination of the stochastic simulation algorithm, where we draw the time for the next adoption event, and a numerical scheme for stochastic differential equations (e.g. Euler–Maruyama) to update the positions of agents between the adoption events [24]. The computational cost can be very high, especially for small time steps and large agent numbers. As soon as there are second-order interactions, the effort does not scale linearly with the number of agents. This motivates us to consider approximate modeling approaches which reduce the numerical effort.

To illustrate the model reduction steps of the next section we introduce now the following guiding Example 1 for our ABM dynamics.

Two-status Dynamics in a Double-Well Potential

We consider the continuous space for the agent movements to be and the discrete status space to be . The movement of a single agent is defined as a diffusion process and is given by the following stochastic differential equation with (see Fig. 2), diffusion constant and a standard Brownian motion process in . The movement process of all agents is denoted by . As for the status dynamics, we consider second-order adoptions with rate functions defined by (2), (3), assuming that only the changes into status are possible, namely with rate constant , while transitions back to status are excluded by setting . In the initial state, all agents are in status 1 except for one agent in the left well, given by the subset , which has status 2. The critical transition event that we are interested in is the first time that an agent with status 2 makes the transition from the left to the right well . In Fig. 2(b) we see that this transition happens only after almost all agents in the left well have adopted status 2. This indicates that adoptions within one well happen faster than between two wells, which is due to the metastability of the diffusion process. In the following sections we will focus on this type of dynamics and derive reduced models that preserve its main properties. In particular, we will consider the critical transition event in order to compare the relevant statistics and quantify the approximation error of derived reduced models. ◊

Fig. 2

Two-status dynamics in a double-well potential, see Example 1. The adoption rate constant is chosen as and the interaction radius as (a) Snapshot of an ABM trajectory for a simulation with agents with contour lines of the potential from Example 1 with diffusion constant . The blue dots refer to spatial positions of agents in status 1, the red dots represent agents in status 2. Gray lines show movement trajectories of the agents. (b) Projected trajectory of the same example, showing the temporal evolution of the total number of agents in each of the two wells (given by and ) depending on the status. The dashed black line refers to the time of the snapshot.

Stochastic metapopulation model

The diffusion process from Example 1 exhibits metastable behavior: an agent remains for a comparatively long period of time within one well of the potential before it eventually jumps to the other well. For such metastable dynamics we consider a coarse-graining such that on the longest time-scale the dynamical properties of the original process are well approximated. More precisely, we assume that there exists a partition of the spatial domain into disjoint metastable sets , i.e., we assume rare transitions between the sets and comparatively fast mixing inside the sets. Due to the fast mixing assumption we do not need to distinguish between the agents positions further than being in one of the metastable sets. Following these ideas, we will introduce here a stochastic metapopulation process, where the considered metastable sets (of the original dynamics) will represent the subpopulations and the rare transitions between the subpopulations will be reduced to jump dynamics. The dynamics within the subpopulations are affected by these transition events and by the adoption dynamics within subpopulations themselves. For more details, see [18] where such a coarsening has been proposed for chemical reaction–diffusion systems. Unlike in [18] where we focused on reactive particles, we will here consider the whole metapopulation state including several different statuses. In the following, we will first describe the coarse-graining approach and then introduce details of the resulting stochastic metapopulation model.

We denote the system state of the stochastic metapopulation model at time by an matrix , where refers to the number of members (agents) of subpopulation in status at time . The set of all possible system states, given the total number of agents, is denoted by , A member of status in subpopulation transitioning to subpopulation at time leads to the immediate change , where is an matrix with all entries zero except for the entry at index being one. Similarly, an adoption event from status to status in subpopulation leads to jump in the system state of the form .

Let denote the probability to find the system in state at time . In analogy to Eq. (7) of the ABM we consider the equation for operators given by where for is the rate to go from to by a spatial transition event between the subsets, while . Analogously, for is the rate to go from to by an adoption event, while . In the following we specify the shape of and and show their connection to the ABM operators and by means of Galerkin projection methods.

Going from ABM to stochastic metapopulation

Given the assumptions of well-mixed behavior, we will make use of the approach given in [18] and construct a Markov state model with respect to the movements of the agents by applying a Galerkin projection (derived from the metastable partition) to the generator defined in (7), see also [17]. In particular, we will study how this projection acts on the transition rate function and on the rate functions for both first-order and second-order adoption events. For deriving the analytical results, we will focus on the case of a full partition of the state space, i.e., it holds for disjoint sets .

For any we define the indicator ansatz functions with where denotes Kronecker delta and set based indicator functions, as well. That is, has the value whenever there are, for each , exactly agents with position and status , otherwise it is zero. These ansatz functions are thus non-negative and fulfill for all , thus they form a partition of unity [29].

Next, we define the inner product of two functions as where denotes the Lebesgue measure. For the indicator ansatz functions defined above we observe that it holds for , while for we have , where denotes the constant -function on .

By means of this inner product we can consider the full-partition projection to the ansatz space given by [29] Given any linear operator , a Galerkin projection with yields the projected operator . The goal is now to find the matrix representations and (see Eq. (10)) of the projected operators and for the operators and defined in (5), (6), respectively.

At first, we consider the spatial dynamics. Define where refers to the standard scalar product for functions in and denotes the constant -function on  [18].

The matrix representation of the projected generator is given by with

The proof can be found in Appendix Appendix A.2

Now, we derive the matrix representations of the projected generator for the interaction dynamics. Here, we consider the two fundamental cases of first- and second-order adoptions separately.

For first-order interactions with adoption rate functions given in (1), we define the conditional expectation of given that : Then, we obtain the following result.

For first-order adoptions with an ABM rate function of the form (1) , the projected generator has the matrix representation with where

The proof can be found in Appendix Appendix A.3.

For the second-order status-changes, we consider the propensity function defined in (2). Assuming the interaction distance to be small (compared to the size of the sets ), second-order adoptions will mainly take place between agents of the same subpopulation as these agents are located relatively close to each other in space. Nevertheless, near the boundaries of the spatial subsets also agents of different subpopulations can be close enough to interact with each other, even if the interaction distance is chosen to be small. However, in a metastable system, the sojourn of agents near the boundaries becomes unlikely and the probability for cross-over interactions (between different subpopulations) approaches zero. More precisely, let for given in (4) denote the conditional probability for two agents to be close enough to each other to interact, given that they are located in the sets and , respectively. Given the number state of the SMM system, we define as the equilibrium propensity for an adoption event to take place in subpopulation by means of a cross-over interaction with agents from a different subpopulation , and as the macroscopic rate constant for adoptions within subpopulation . Using these definitions, we obtain the following result.

For second-order adoptions with an ABM rate function given by (2), (3) , the projected generator has the matrix representation with where

The proof is given in Appendix Appendix A.4.

Theorem 3 shows how the projected rate functions of second-order adoptions can be decomposed into: (1) a part coming from adoptions that take place between agents of the same subpopulation with propensities given by ; and (2) a part coming from adoptions that take place between agents of different subpopulations. As discussed above, in a metastable system and for a good choice of sets , the probability of cross-over interactions is small and thus the value of will be negligibly small; in fact, as we will see below in Example 1, it often is orders of magnitude smaller than the discretization error of the Galerkin projection.

In the SMM we only take the first type (1) of interactions, i.e., we assume that there are no (cross-over) adoptions taking place between agents of different subpopulations. Using the results from Theorem 1, Theorem 2, Theorem 3, the SMM equation given by (9), (10) can be written as the following spatio-temporal master equation: where the first two terms on the right-hand side refer to the change caused by the exchange between subpopulations (given by the operator ), while the other two lines are describing the change through status adoptions inside the subpopulations (referring to operator ).4

Using the definitions of the functions given in Theorem 2, Theorem 3, the interaction propensities agree with the standard law of mass-action from the chemical context [20]. This is due to the fact that we assume the agents to interact independently of each other (and of the overall system state) — which we do by choosing the ABM adoption rate functions according to Eqs. (1), (2). As for the spatial dynamics, state-of-the art metapopulation models [30] usually assume the commuting flow between two subpopulations and to be of the form for exponents which tune the dependence with respect to each subpopulation size. In our setting, we assume that the spatial movement of each agent is independent of the population sizes, which corresponds to setting and .

Extension to core set approach.

The Galerkin projection does not have to be constructed on a full partition of the state space using indicator function as ansatz functions, but can also use ansatz functions that form a partition of unity. Appropriate ansatz functions of this form are defined via so-called core sets, i.e., sets that do not form a partition of state space but cover only the core areas of the metastable sets [16], [31], [32]. For example, for a diffusive process in a potential energy landscape, the core sets are given by the vicinities around the energy’s local minima (i.e., the valleys or wells of the landscape), while the transition regions around the energy’s local maxima are not explicitly assigned to any core set. The advantage of considering a core set approach is a reduced approximation error in the estimated transition rates compared to the full partition case [16], [32]. The ansatz functions associated with the core set is its committor function [33] , with denoting the probability that the agent visited the core set last, conditional on being in position , see [29]. Thus, for agents we have that and . Agents that are in the transition region can be affiliated to several core sets with different probabilities, but such that it holds .

For the core set approach, we have to redefine the conditional probabilities given in (16) according to Since in the core set approach we choose the sets to be only the core areas of metastability, the probability of cross-over interactions of agents from different subpopulations is extremely small and thus, the value of is even smaller than in the case of a full partition.

Approximation Quality

The step from the full-scale ABM to the SMM (20) involves two approximations: the discretization error originating from the Galerkin projection, and the error resulting from neglecting the cross-over interactions between different spatial domains. While latter error can easily be monitored by estimating the neglected cross-over rates, controlling the discretization error is more difficult. For the cases where the spatial movement exhibits metastable sets, estimation of the discretization error is possible, see [31] for the core set approach, but requires sufficient ABM simulation data.

In order to illustrate how the stochastic metapopulation model can approximate the ABM, we now return to our guiding example.

Continued

Given the double-well potential shown in Fig. 2 and using the Markov state model approach, we partition the space of movement into the two core sets and and the transition region . The jump rates and between the two subpopulations are the transition rates between and . We analyze the quality of the SMM approximation for two different values of the diffusion constant and , where the first case is more metastable than the other. We compare the SMM process to the projected ABM dynamics regarding the temporal distribution of a critical transition event given by the first agent with status switching from one to the other subpopulation (meaning for the projected ABM that an agent of status who visited core set last reaches core set for the first time), see Fig. 3. As this transition has a major impact on the overall dynamics, we consider the total approximation error to be small if the difference in the distributions of this critical transition event time is small. We observe that for smaller the approximation is better, which is due to an increase in metastability of the dynamics.

Fig. 3

Distribution of the critical transition time for the ABM (blue) and the SMM (orange) given in Example 1, sampled over MC-simulations. The overlap of the two distributions is colored brown. For small the distribution is very well matched, while for larger the critical transition happens faster in the SMM.

This becomes even more clear when comparing the temporal evolution of the average number of agents of status 2 in the two subpopulations, see Fig. 4. For the smaller value , the first-order moments agree very well (Fig. 4(a)). In contrast, there is a significant difference between the model outcomes regarding these first-order moments for the larger diffusion constant (Fig. 4(b)). This is due to the approximation quality of the Markov state model being worse because the diffusion process is less metastable and thus the first spatial transitions in the SMM happen too fast on average. The deviation of the spatial transition dynamics is also the main contribution to the approximation error of the critical transition event time. ◊

Fig. 4

Comparison of the time-dependent mean population size of agents with status over MC-simulations for the projected ABM and the corresponding stochastic metapopulation model. The projection is based on a core set approach with and given in Example 1 and defining two subpopulations denoted by SP1 and SP2. In scenario (a) the solid blue and dashed blue line are indistinguishable due to the good approximation quality.

Stochastic simulation.

For the stochastic metapopulation model we realize the sampling with the stochastic simulation algorithm , which produces statistically exact realizations of the process, without any numerical approximation error [19]. The effort is much lower than for ABM-simulations but still scales at least linearly with the number of agents. This means that for large population numbers the sampling of quantities of interest can still be infeasible due to high computational costs. This motivates us to consider a further model approximation as discussed in the following subsection.

Piecewise-deterministic metapopulation model

For the case of a large population of interacting agents, the stochastic simulation of the stochastic metapopulation dynamics becomes computationally very expensive, because it tracks every single adoption event. Here, a further model reduction can be very useful. Given that the number of agents is large in each subpopulation, we can apply standard convergence results and approximate the jump process which describes the internal adoption dynamics by a deterministic evolution equation [21]. Such approximations, which are based on the law of large numbers, are well-known in the context of chemical reaction systems, where they are used to reduce the model complexity for systems with large molecular populations [20]. For the exchange process between the subpopulations, on the other hand, we assume that transitions are quite rare in time and occur sporadically at variable time points. In order to reduce the computational effort for simulations, while keeping the discrete, stochastic nature of the (inter-nodal) transition events, we approximate the overall dynamics by a piecewise-deterministic Markov process, see [27], [34], [35], [36] for details. This model will be called piecewise-deterministic metapopulation model (PDMM).

Going from stochastic metapopulation dynamics to PDMM

The stochastic process given by (20) can be rewritten in a pathwise notation of the form where and refer to independent, unit-rate Poisson processes [35], [37].

Assuming that the jumps induced by the Poisson processes , which refer to spatial transitions between metastable domains, occur with much less frequency than jumps induced by the Poisson processes referring to the adoption dynamics within the subpopulations, we can apply standard convergence results for Markov processes [38] in order to approximate the stochastic dynamics given by the second line of Eq. (22) by deterministic dynamics and obtain the PDMM process given by the equation This means that the population dynamics on the local scale are modeled by a system of ODEs (adoption dynamics given in the second line), while the rare interactions between subpopulations are modeled as a stochastic jump process (first line), just as in [27]. It is well-known [20], [39] that the relative error produced by this approximation (relative with respect to the population size) decreases with an increasing number of agents with a Monte Carlo like rate. We will here consider a finite sized population which is large enough for the approximation to be reasonable.

For our example of two-status dynamics in a double-well potential, the deterministic status-adoption dynamics from the second line of (23) are given by for , where denote the time points of two subsequent stochastic transition events induced by the first line of (23). Using the definition (19) of , we get the following ODE for the number of agents in subpopulation having status : for . Let denote the total number of agents in subpopulation at time . Between two transition events this number is constant, so we can substitute in Eq. (24) to arrive at The solution is given by the logistic function, that is, we obtain an analytical solution for . Treating the diffusive transitions between the subpopulations as stochastic events which induce jumps in the state of the PDMM process, we obtain trajectories as depicted in Fig. 5(b).

Fig. 5

Comparison between (a) SMM and (b) PDMM trajectories for . In (b) the rare jump events are marked by vertical dotted lines.

The main difference between the trajectories (see Fig. 5) is that in the stochastic metapopulation model we have many discontinuous jumps, while in the PDMM only a few remain. The error for the estimation of the critical event time caused by the PDMM approximation (see Fig. 6) is small compared to the error originating from the spatial discretization (see Fig. 3). Even though our population number of 100 individuals is not very large, the critical transition time distribution of the SMM is already approximated well by the PDMM for both choices of . ◊

Fig. 6

Distribution of the critical transition time for SMM (blue) and PDMM (orange) sampled over 10 000 MC-simulations. The overlap of the two distributions is colored brown. For both values of the distribution is very well matched.

In order to simulate a PDMM process given by (23) one has to simultaneously integrate the deterministic flow of the ODE part and the rate functions for the stochastic jumps, see [37]. Using the concept of the temporal Gillespie algorithm [23] one can thereby determine the time point of the next stochastic jump. After a jump event the initial conditions for the ODE are updated accordingly.5 The numerical effort does not scale significantly with the number of agents, and thus, from the approaches considered in this paper, the PDMM is the most efficient choice for simulating systems with large agent numbers, see Fig. 7 for a comparison of the computational effort for the three modeling approaches.

Fig. 7

Numerical effort for the simulations of the guiding example system of Example 1 for different choices of agent numbers , depending on the number of agents. Even for low agent numbers the approximate models are at least two orders of magnitude more efficient. For low agent numbers the PDMM effort is higher compared to the SMM as the critical transition event happens significantly later and thus more time steps are needed in the PDMM computation.

Modeling COVID-19 epidemic spreading

In this section, we discuss how the above presented approaches can be applied to the modeling of epidemic spreading. Agent-based models, as the most detailed models, have been extensively used for data-driven modeling with newly available large data-sets, e.g. on human mobility [1], [40]. However, most ABMs lack formalizations that would allow for more thorough mathematical analysis. Also, due to their large computational complexity, ABMs are usually focusing on epidemic spreading on smaller spatial scales, e.g. cities. Stochastic metapopulation models offer precise mathematical descriptions but, to the best of our knowledge, have not yet been investigated in the context of epidemic modeling. Piecewise-deterministic metapopulation models have been largely applied to model pathogen spreading [25], [26] due to (1) their smaller computational complexity compared to ABMs; (2) better spatial resolution as opposed to ODE models, which are typically used. Despite the enormous recent advancement in adapting these models to realistic epidemic scenarios, this topic is still at the infancy of its development. In particular, mathematical model formalization with respect to data-based descriptions on different levels needs to be explored in much more detail.

The results presented above show us how ABMs, SMMs and PDMMs are coupled. However, integrating all three types of models in one application would go beyond the scope of this article. Therefore, we will focus here on the most coarse-grained model, namely the piecewise-deterministic metapopulation model, and demonstrate how it can be applied to conceptually analyze COVID-19 spreading dynamics.

In the case of COVID-19, metapopulation models provide a good approximation of the original epidemic dynamics since the metastability assumption can be observed in mobility data [41], e.g. in rare spatial transitions between different cities or countries caused by mobility restrictions. Also large population sizes are realistic, such that an approximation by piecewise deterministic dynamics is justified and allows to drastically reduce the model complexity compared to the other two approaches. Moreover, the PDMM can more easily be calibrated to real world data than the underlying ABM, as we have good estimates for rates on the population scale, but not for individual interactions.

For simplicity, we will consider two subpopulations that have frequent local interactions and rare transitions in-between. The model will be calibrated to the parameters estimated in recent studies, for details see below. Note that in this paper, the model is not applied to analyze a particular real-world dataset, but representative results are used in order to show how the model can be applied to possible real-world scenarios. Additionally, we will analyze the effectiveness of different containment measures taken within subpopulations and of global measures that impact the spreading between subpopulations. However, the main goal of this section is not to identify the optimal choice of containment measures, but to demonstrate the applicability of the PDMM on stochastic spreading processes and their efficiency on large real-world systems.

The PDMM dynamics

As a first approach to formulate the piecewise-deterministic spreading dynamics, we will consider the rate constants , which define the propensity of agents to change their status, to be independent of the evolution of the process. In Section 3.2 we will generalize the dynamics by letting these constants depend on time and on the process’ history.

The system state

On a local scale of each subpopulation, we use a compartmental Susceptible–Exposed–Infected–Recovered–Deceased (SEIRD) model [8], [42] to describe the COVID-19 dynamics. In this model, each individual can be in one of five possible statuses: susceptible (S), exposed (E), infected (I), recovered (R) and deceased (D). Susceptible individuals are the ones who have not yet been in contact with the virus and have no immunity against it. After being exposed to the virus, a susceptible individual is first in an asymptomatic status E that changes to a symptomatic status I after the end of an incubation period. Recovered and deceased individuals are considered to be immune for the time scale of our model and are no longer able to be infected and transmit the disease.

Given these statuses as well as a set of subpopulations, a possible state of the system has the form where denotes the number of individuals in status within subpopulation .

Deterministic local interaction dynamics

We assume that exposed individuals are already able to transmit the virus [43], so the status change from to can be caused by the second-order interactions of type either or , with respective transition rate constants and , see Fig. 8 for an illustration. In general, it is possible to distinguish between infection rates from a contact with an exposed individual and from a contact with an infected individual. However, for simplicity, we assume here that infectiousness is constant from the moment of exposure until recovery or death, i.e. .

Fig. 8

Visualization of the SEIRD model. Black arrows stand for possible status transitions, blue arrows indicate an impact by interaction.

Given the system state , the rate functions for the second-order status-change are given by and respectively, see (19).6 The remaining status transitions that we consider are given by first-order events of the form and with respective rate constants . For the transition we accordingly obtain and equivalently for and The resulting ODE-system describing the local interaction dynamics within a subpopulation is then given by

Stochastic dynamics for spatial exchange

While the local interaction dynamics are given by deterministic evolution equations (25), the spatial transitions between the subpopulations are described by stochastic jump events. At any point in time , an agent of subpopulation can switch to another subpopulation , which induces a discrete change in the system’s state of the form depending on the agent’s status . (For we naturally assume that jumps cannot take place.) The PDMM process combines these discrete, stochastic jump events between the subpopulations with the ODE dynamics (25) for local status-changes, which in total, leads to a stochastic process , described by an equation of the form (23). The terms in the first line of (23) thereby correspond to an integrated version of the ODE (25) for each , and the second line describes the spatial jumps for given rate constants between subpopulations and . In the following subsection, both the rate constants for the local interactions and the rate constants for the spatial transitions will depend on the evolution of the overall stochastic process .

modeling choices presented above are made for simplicity and to demonstrate how the PDMM can be used to conceptually analyze COVID-19 spreading. Using rich data-sets and extensive literature on COVID-19 pandemic, our model can be easily extended to account for more realistic scenarios. For example, considering additional compartments such as symptomatic, asymptomatic, quarantined individuals [44]; including more general infection rates with possible time dependency [45], adding demographics information [46], [47], introducing vaccination effects [46] are only some of the extensions that would make this model more realistic.

Adaptive regulation of rate constants

Until effective pharmaceutical treatment for COVID-19 is found, a lot of effort is taken to slow down the virus spreading by introducing measures for reducing social contacts. This is achieved, for example, by targeting the individual interactions (e.g. social distancing and wearing masks), by reducing the number of interactions (e.g. closing of schools, offices), but also by introducing measures on a global level, such as travel bans between countries and continents.

Having this in mind, the choice of transition rate constants that are independent of time and of the process’ evolution (as assumed in Section 3.1) appears to be unrealistic because containment measures are taken depending on the dynamics in order to influence the future evolution of the process. This is why we will in the following consider rate constants that are adapted in the course of time according to given rules.

The transition rates are depending on the local contacts within a population, which are changing over time with the installation of measures. Additionally, the dependence can also be on the interaction type, e.g. exposed individuals might be less contagious than infected ones or symptomatic cases could cause less infections than expected because they have already reduced their number of contacts. In order to include many types of possible dependencies, we define the transition rates in a quite general way as functions of the process’ history and time That is, these rates not only depend on the current state but also on the history of the process. By this we can define rules such as implementing a strict lock-down when case numbers are rising for the first time. In our model, the rate for developing symptoms is a constant, while the recovery and case fatality rates depend on the capacity of the health care system of a population, so they are defined as state-dependent rates and .

One goal of the introduced measures is to control the number of infections such that the limits of the health care system are not reached. As part of the global measures, the transitions between the subpopulations will be reduced. Thus, the spatial transition rates between subpopulations are defined as functions of the entire process’ history and time for :

Concrete choice of rate constants for status changes

When modeling the implementation of measures for virus containment, we will assume that each measure is followed by a phase for which the infection rate remains constant. The transition between the phases can be triggered by deterministic as well as stochastic events, e.g. by the process crossing a threshold number of infections for the first time. In total, we will consider three different phases:

Initial phase: In the beginning of the pandemic, the infection rates and have values and , respectively, and the interaction dynamics start with an unmitigated spreading.

Strict measures phase: The first measures to reduce the infection rates are taken in subpopulation as soon as the number of infected individuals crosses a critical value for the first time. That is, the strict measures start at the random first hitting time with value in case of for all . In this phase, the infection rates and are reduced by a factor . These strict measures are kept until the number of infected individuals falls below the critical value , i.e. until the random time point

Moderate measures phase: After the strict measures are lifted, the interactions inside the population do not go back to normal, i.e. to the values from the initial phase. Instead, we introduce moderate measures where the infection rates are scaled with a factor , s.t. , allowing for more contacts than in the previous phase. These measures are maintained for the remaining time of the model even if the number of infected individuals crosses again the value .

Taken all together, this means that the infection rate function is defined by and equivalently for , possibly with different reduction factors . More generally, also the rate constants may depend on the subpopulation , but here we omit the corresponding indices for the purpose of simplicity.

In order to make our model more realistic, we include in each subpopulation a limited health care capacity given by a threshold . We assume that the case fatality rate increases from a given value to another value if the number of infected individuals exceeds this threshold , giving Vice versa, the recovery rate is reduced in case of an exhausted health care capacity, such that for constants . Additionally, within each subpopulation we consider to be constant, i.e. .

Finally, we assume that exposed individuals develop symptoms after an incubation period of average length and set for all .

Concrete choice of rate constants for spatial transitions

The global spatial transition rate functions between the subpopulations will depend on the local phases within each of the subpopulations. More precisely, we define to be the first time that one of the subpopulations initiates the lock-down phase and to be the first time that both subpopulations have ended the lock-down phase. Whenever in one of the subpopulations the strict measures are applied, the spatial transition rates are reduced by a factor . After the strict measures have ended in both populations, the spatial transition rates are scaled by a factor , where . Thus, the spatial transition rate functions are defined as for . We assume that people with symptoms do not travel, i.e. we set for all independently of time.

PDMM-based simulations of COVID-19 spreading

We simulate the dynamics for model scenarios with differences in the infection and spatial transition dynamics. In particular, we compare the following three scenarios:

Scenario 1: Choose constant infection and spatial transition rates, which can be interpreted that no measures are implemented.

Scenario 2: Let the infection rate depend on the process history as defined in (26), but assume constant spatial transition rates between subpopulations. This corresponds to introducing local measures to control the infection dynamics within subpopulations, but no additional travel restrictions in between.

Scenario 3: Combine the measures, i.e., let both infection and spatial transition rates change according to the epidemic dynamics following the rules defined in (26), (27).

Parameter choices

Recently, lots of research has been done on inferring the parameters of COVID-19 dynamics from available data. However, for most parameters there is a wide range of estimates and thus the choices for a conceptual model can seem arbitrary. Since parameter estimation is not at the core of this manuscript, but our modeling approach is, we will choose the parameters based on a few recent publications [48], [49], [50], [51], [52].

Parameters for status change : The average incubation period is estimated to be – days [50], [51], so we choose .

Parameters for status change and : The average time for transition from I to either R or D will be set to days [50], [51]. For the infection fatality rate of our model, we will use the estimate from [52], which leads to the choices of and for the recovery and case fatality rates of the populations. We choose and accordingly such that is fulfilled.

Parameters for status change : Here we consider two reactions that can lead to the status change , namely and . As discussed before, for simplicity, we assume that and thus . Estimates for the initial reproduction number vary depending on the region of choice [53] as well as on the estimation method [54]. This leads to a wide range of possible parameter choices that have the highest impact on the model outcome. For our model, we use which corresponds to the estimate for the New York City in [53], assuming our subpopulations to be well-mixed and an urban area like NYC to meet this assumption. This choice of leads to the infection rate . The remaining parameters will be subject to changes due to different containment measures. During a period of strict measures phase we will reduce the infection rate to 10% of the original value by setting a scaling factor for . In the phase of moderate measures we assume more interactions which lead to an increase of the infection rate. For illustration purposes, we consider different choices of moderate measures in each subpopulation, such that we set values of the infection rates to be 30% and 40% of the original , i.e.  and . For both subpopulations, the infection threshold for the first hitting time event is chosen to be 2% of the initial total population number in subpopulation , i.e.  and the threshold for the capacity of the health care system is reached when 10% of are infected, i.e. .

Parameters for spatial transitions: The spatial transition rates are chosen to be , which corresponds in our example to 3 out of 10 000 agents transitioning per time unit and fits to the assumption of a metastable setting with slow transitions between subpopulations compared to the infection dynamics within subpopulations. When at least one subpopulation is in the strict measures phase, we introduce travel restrictions by reducing the spatial transition rate to 5% of the original value, i.e. we set . When both subpopulations are in the moderate measures phase, we moderately relax the travel restrictions by increasing the spatial transition rate to 50% of the original value, i.e. we set .

Simulation results

For the initial population sizes we choose the values . We start with one member of the subpopulation 1 (SP1) being in status E and all other members of SP1 and subpopulation 2 (SP2) being in status S. The critical transition event is the first time when an individual with status E jumps from SP1 to SP2. In Fig. 9, we see one outcome of Scenario 1, where no containment measures are taken. We observe one wave of infections in each subpopulation with the number of I cases quickly going above the threshold and staying there for a considerable amount of time. Almost all population members get infected. Due to the increased case fatality rate , in the end 3.8% of the total population are in status D. This is a scenario that in reality should be avoided, for example by flattening the curve through the implementation of containment measures.

Fig. 9

Trajectory of a PDMM simulation for Scenario 1 (no measures). The dotted magenta line marks the critical transition event, the horizontal red line marks the threshold of health care capacity.

As a result of the local measures that are present in Scenario 2, the curve of infections shows two smaller waves instead of one large wave, see Fig. 10. In SP1 the number of infections stays below the threshold during the whole simulation period, while in SP2 the number of cases crosses during the peak of the second wave. This is due to a higher infection rate within SP2 in the phase of moderate measures, which leads to a higher total number of infections and more fatal cases in SP2 at the end of the simulation. Nevertheless, the outcome in both subpopulations is a much smaller number of I and D individuals than in Scenario 1. The same is true for Scenario 3 which has the same local measures.

Fig. 10

Trajectory of a PDMM simulation for Scenario 2. The dotted magenta line marks the critical transition event, the horizontal red line marks the threshold .

Additionally, in Scenarios 2 and 3 the shape of the infection curves determined by the internal population dynamics is the same (see Fig. 11). However, the distribution of the critical transition time that starts the epidemic in SP2 is considerably different. Due to the introduction of travel restrictions in Scenario 3, we observe a later first infection in SP2 compared to the one from Scenario 2, see Fig. 11(b).

Fig. 11

Comparison between infection curves for single trajectories obtained in Scenarios 2 and 3. The vertical dotted magenta line marks the critical transition event, the horizontal red line marks the threshold . The dashed lines refer to the development in the subpopulation 1 and solid lines to subpopulation 2.

In order to compare the critical transition time distributions for different containment measures, we run MC-simulations for each of the three scenarios, see Fig. 12. For approximately one third of the simulations (no matter which scenario) the critical transition happens before the time when the measures are introduced in scenarios 2 and 3. After this time point we can observe the differences in the shapes of the critical transition time distributions. Namely, compared to Scenario 1, we observe for Scenario 2 a larger number of critical transitions happening later in time. This is due to the influence of the number of active cases in SP1, which is declining much faster in Scenario 2 compared to Scenario 1 due to the local measures. The mean time for the first infection in SP2 is 24.5 days for Scenario 1 and 43.9 days for Scenario 2. Due to the reduced spatial transitions in Scenario 3, the probability for a critical transition is much smaller during the measures. As a result, in (out of ) MC-simulations the critical transition event did not occur at all and the virus was successfully contained in SP1. Conditioned on the transition happening before the end of the simulation period, the mean critical transition time was 78.6 days, which shows the benefit of the travel restrictions on the spreading dynamics.

Fig. 12

Critical transition time distribution for different scenarios of the model. The number of MC-simulations for each scenario is .

In our conceptional model the outcome of the epidemic was improved by the implementation of measures. The spatial separation into multiple subpopulations has a high impact on the dynamics of the spreading process, especially when considering the scenario of combined measures where the spreading between populations could be delayed by a long time and sometimes even prevented.

Conclusion

In this work we introduced a hierarchy of modeling approaches for spatio-temporal dynamics of interacting agents. We showed how the stochastic metapopulation model can be derived from the underlying spatially continuous agent-based model by means of a Galerkin projection of the dynamics, considering both a full-partition approach and a core set approach. Especially, we specified the form of the projected generators and derived equations for the relation between the corresponding adoption rate functions for both first- and second-order status changes. In our guiding example, we analyze how the approximation error depends on the metastability of the dynamics.

Given the stochastic metapopulation model, we investigated a further approximation by piecewise-deterministic dynamics, where only the spatial transitions between metastable areas are modeled as stochastic jumps while the adoption dynamics within each subpopulation are approximated by deterministic dynamics. We analyzed the approximation quality as well as efficiency, finding that for large numbers of agents the PDMM delivers convincing results combined with enormous decrease in computational effort.

Based on this insight, we formulated a PDMM for the spatio-temporal spreading dynamics of COVID-19 and analyzed the impact of different measures. We compared three basic scenarios with respect to critical transition times of the virus spreading between two subpopulations, finding out that the spatial separation into subpopulations can have a high impact on the epidemic spreading process, e.g. by means of local measures and/or traveling restrictions.

Finally, in this paper we showed that, for large number of agents, stochastic metapopulation models and in particular PDMMs represent good approximations of ABMs which can be achieved with much less computational power, allowing for faster simulations of many different scenarios. Compared to other state of the art modeling approaches that were not discussed in detail in this manuscript, e.g. ODE-based models [9], [10], [11], [12], [13], [14], in real world applications, SMMs and PDMMs appear to be better suited as they include both a certain level of stochasticity and a spatial resolution on a mesoscale and still have a reasonable computational cost.

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.