key: cord-0486601-vyc9xxnc
authors: Kim, Daewa; Quaini, Annalisa
title: A 2D kinetic model for crowd dynamics with disease contagion
date: 2021-07-23
journal: nan
DOI: nan
sha: 9932c6a0bbd97195a88041f9a954a9684e6c5a0b
doc_id: 486601
cord_uid: vyc9xxnc

We focus on the modeling and simulation of an infectious disease spreading in a medium size population occupying a confined environment, such as an airport terminal, for short periods of time. Because of the size of the crowd and venue, we opt for a kinetic type model. The paper is divided into two parts. In the first part, we adopt the simplifying assumption that people's walking speed and direction are given. The resulting kinetic model features a variable that denotes the level of exposure to people spreading the disease, a parameter describing the contagion interaction strength, and a kernel function that is a decreasing function of the distance between a person and a spreading individual. Such model is tested on problems involving a small crowd in a square walkable domain. In the second part, ideas from the simplified model are used to incorporate disease spreading in a kinetic theory approach for crowd dynamics, i.e. the walking speed and direction result from interaction with other people and the venue.

The goal of this paper is to extend to 2D a 1D model presented in [30] to study the early stage of an infectious disease spreading in an intermediate size population occupying a confined environment, e.g. a ER wait room, for a short period of time (minutes or hours). In this context, classical epidemiological models fail as they rely on averaged large population behaviors over a long time span (weeks or months). We focus on a kinetic theory approach [6, 9] for a disease that it spreads with close proximity of individuals, like, e.g. measles, influenza or COVID-19. In the first part of the paper, we present a model for the disease spreading based on the simplifying assumption that people's walking speed and direction are given. In the second part, ideas from the simplified model are used to incorporate disease spreading in a kinetic theory approach for crowd dynamics.

The COVID-19 pandemic has motivated a large body of literature in mathematical epidemiology. Several possible models have been investigated, including individual-based Markov models (see, e.g., [20] ), Susceptible, Exposed, Infectious or Recovered (SEIR) models (see, e.g., [25] ), and networks of nodes where each mode dynamics has a SEIR structure (see, e.g., [15, 21] ). Some works (e.g., [31] ) model the evolution of the COVID-19 epidemics in combination with other phenomena, such as fear. Other works focus a multiscale modeling approach. A notable example that features an interdisciplinary framework including applied mathematics, immunology, economics and virology is [8] . Similarly, modeling and simulation of disease spreading in pedestrian crowds has recently become a topic of increasing relevance. In [19] , contact evolution is combined with a stochastic infection spread model to simulate disease spreading in queuing pedestrians. Agent-based numerical simulations of pedestrian dynamics are used in [23] to assess the behavior of pedestrians in public places in the context of contact transmission of infectious diseases and gather insights about exposure times and the overall effectiveness of distancing measures. Finally, we mention the combination of a microscopic force-based model with a contact tracking method in [32] to simulate the initial spreading of a highly infectious airborne disease in a confined environment.

The model that we propose for disease spreading in a walking crowd takes inspiration from the work on emotional contagion (i.e., spreading of fear or panic) in Ref. [13, 36] . In order to focus on the spreading mechanism, we first introduce a simplifying assumption: the walking speed and direction are given. The key features of the resulting simplified model are:

-a variable that denotes the level of exposure to people spreading the disease, with the underlying idea that the more a person is exposed the more likely they are to get infected; -a model parameter that describes the contagion interaction strength; -a kernel function that is a decreasing function of the distance between a person and a spreading individual. We show preliminary results for a problem involving a small crowd in a square walkable domain with no obstacles or walls.

In the second part of the paper, the simplifying assumption is removed: the walking speed and direction result from a kinetic model of crowd dynamics, which incorporates a term to account for disease spreading. Kinetic models are derived from observing the system at the mesoscale (see, e.g., [1, 5, 6, 7, 10, 11, 12] ), i.e. the intermediate scale between the macroscopic one (see, e.g., [17, 26, 34] ) and the microscopic one (see, e.g., [3, 16, 18, 24] ). The framework for mesoscale models comes from the kinetic theory of gases, with the main difference that interactions of "active" particles are irreversible, non-conservative and, in some cases, non-local and nonlinearly additive [2] . In general, one derives the kinetic model from interactions at the microscopic scale and the same modeling principles can lead to macroscopic (or hydrodynamic) models. See, e.g., Ref. [4] . We remark that the approach to model crowd dynamics with disease contagion presented here is different from the one in Ref. [29] , where two models (one for pedestrian dynamics and one for disease spreading) are coupled.

The particular kinetic model for crowd dynamics that we extend to account for disease spreading was first presented in [28] . Based on earlier works [1, 10] , the model's main features are -discrete walking directions; -interactions modeled using tools of stochastic games; -heuristic, deterministic modeling of the walking speed corroborated by experimental evidence [33] . In [28] , the model has been shown to compare favorably against experimental data related to egression from a facility of a medium-sized group of people (40 to 138 pedestrians) [35] . In addition, we demonstrated that realistic scenarios, such as passengers moving through one terminal of Hobby Airport in Houston (USA), can be reproduced [30] .

The paper is organized as follows. We introduce our simplified contagion model in Sec. 2 and its numerical discretization in Sec. 3. Numerical results are shown in Sec. 4. Sec. 5 presents the kinetic theory approach for crowd dynamics extended to incorporate disease spreading. Conclusions are drawn in Sec. 6.

We start from an agent-based model. At the microscopic level we consider a group of people divided between healthy or not spreading the disease yet ( ℎ ) and actively spreading ( ), with = ℎ + . If person belongs to the group of healthy (or not spreading) people, we denote with ∈ [0, 1) their level of exposure to people spreading the disease, with the underlying idea that the more a person is exposed the more likely they are to get infected. If person is actively spreading the disease, then = 1 and this value stays constant throughout the entire simulation time. Let [ 0 , ] be a time interval of interest. Let ( ) = ( ( ), ( )) and ( ) denote the position and speed of person for ∈ [ 0 , ], with initial position ( 0 ) and initial speed ( 0 ) being given. The microscopic model reads for = 1, 2, 3, . . . , :

where the walking speed and walking direction are assumed to be given. In model (1), * corresponds to a weighted average "level of sickness" surrounding person , with , serving as the weight in the average. We define , as follow

Notice that if person is spreading the disease, the interaction kernel is a decreasing function of mutual distance between two people and is parametrized by an interaction distance . The value of is set so that the value of , is "small" at about 6 ft or 2 m. More details are given in Sec. 4 . The meaning of parameter in (1) is contagion strength: for = 0 there is no contagion, while for ≠ 0 the contagion occurs and it is faster for larger the values of . Note that obviously the level of exposure can only increase over time. In addition, the second equation in (1) ensures that the people spreading the disease will constantly have = 1 in time. In order to derive a kinetic model from the agent-based model (1), we introduce the empirical distribution:

where is the Dirac delta measure. We assume that all people remain in a fixed compact domain Ω for the entire time interval under consideration: ( ( ), ( )) ∈ Ω ⊂ R 3 , for = 1, 2, 3, . . . , and ∈ [ 0 , ]. Let P (R 3 ) be the space of probability measures on R 3 . The sequence {ℎ } is relatively compact in the weak * sense (see, e.g., [14] ). Therefore, there exists a subsequence {ℎ } such that ℎ converges to ℎ with weak * -convergence in P (R 3 ) and pointwise convergence in time as → ∞.

where · , means integration against both and . Let us define

where · means integration only in . Then, we can rewrite eq. (3) as

Via integration by parts, eq. (4) leads to

where * is the local average sickness level weighted by (2):

Letting → ∞, we obtain the limiting kinetic model from the subsequence ℎ :

where ℎ( , , ) is the probability of finding at time and position a person with level of exposure if ∈ [0, 1) or a person spreading the disease if = 1.

It is convenient to switch to non-dimensional variables. Let be the largest distance a pedestrian can cover in domain Ω and the largest speed a person can reach. Then, we can define characteristic time = / . Moreover, let be the maximum people density per square meter. The non-dimensional quantities are: positionˆ= / , walking speedˆ= / , timeˆ= / , and distribution function ℎ = ℎ/ . For ease of notation though, we will drop the hat with the understanding that all quantities are non-dimensional unless otherwise specified.

In this section, we present the numerical discretization of model (7) .

Let us start from the space discretization. Divide the spatial domain into a number of cells

] of length Δ and Δ , respectively. The discrete mesh points and are given by

Here, , and are the total number of points in −, − and − directions, respectively.

Let us denote ℎ , , = ℎ( , , , ). The average level of sickness is computed using a midpoint rule for the integrals in (6)

where if at ( , ) there is nonzero probability of finding a person that is spreading the disease:

, otherwise , = 0. We start with a first-order semi-discrete upwind scheme for problem (7) , which reads:

where , , = cos ℎ , , , , , = sin ℎ , , ,

, ,

Next, let us discretize in time. Let Δ ∈ R, = 0 + Δ , with = 0, ..., and = 0 + Δ , where is the end of the time interval under consideration. Moreover, we denote by the approximation of a generic quantity at the time . For the time discretization of problem (9), we use the forward Euler scheme:

To construct a scheme that is of second order, we add a flux limiter. Let be a slope limiter function. for example, one could choose the Van Leer function:

The space discretized eq. (7) now reads:

where

In eq. (17), the subscript k is k-1 if , , − 1 (15) and , , in (16) are defined in (10), while , , + 1 2 is defined in (11) . Thanks to the flux limiter, scheme (14) is of second-order scheme in velocity.

For the time discretization of problem (14), we adopt again the forward Euler scheme

The time step Δ for scheme (12) or (18) is set as

to satisfy the Courant-Friedrichs-Lewy (CFL) condition.

We assess the approach presented in Sec. In all the tests, we take the initial density to be constant in space and equal to . So, the group occupying the walkable domain consists of 700 people.

We set the interaction distance = 1 m since this choice makes the value of the kernel function relatively small at a distance of 2 m (or about 6 ft). See Fig. 1 . To discretize the walkable domain, we take Δ = Δ = 0.1 m, while we set Δ = 0.01. We will consider two values for the contagion strength: = 50 and = 10. The time step is set to Δ = 5 · 10 −5 s for = 50 and to Δ = 25 · 10 −5 s for = 10 according to (19) . First, we keep the group of people still (i.e., = 0) to make sure that the level of exposure evolves as expected. Then, in a second set of tests, we let people walk and observe how the motion affects the spreading of the disease. We run each simulation for ∈ (0, 10] s.

We adapt two initial conditions from [30] :

-IC1: people that are certainly spreading (i.e., = 1) are located outside the circle of radius 1 m centered at (5, 5) m, while within that circle we place people that have certainly not been exposed (i.e, = 0). See Fig. 2 (left) . -IC2: people that are certainly spreading (i.e., = 1) are located outside the circle of radius 2 m centered at (5, 5) m, while the rest of the people located within that circle have certainly not been exposed (i.e, = 0). See Fig. 2 (right) . The difference between the two boundary conditions is that all the healthy people in IC1 are exposed to spreading people, while some healthy people in IC2 are not exposed to spreading people. Thus, we expect that the level of exposure for the people that are centrally located rises faster for IC1 than for IC2. 3 shows the evolution of the distribution density ℎ for initial conditions IC1 and IC2 with = 50, together with a view of the results on a section of the 3D domain. We observe that the level of exposure of the central group of healthy people in IC1 increases quickly. In particular, we see that the the shape from * quickly moves away from being close to a sharp discontinuity and gets closer to a paraboloid, since the level of exposure increases faster for the people closer to the circle separating healthy people from spreading people. As expected, the rise in the level of exposure is much slower for the simulation with initial condition IC2. This is evident when comparing the second and fourth row in Fig. 3 . Moreover, notice that the healthy people close to the center of the circle in IC2 get very little exposure to spreading people.

Next, we consider initial condition IC1 and compare the evolution of the distribution density ℎ for = 50 and = 10. See Fig. 4 . We see that the level of exposure of the central group of healthy people increases faster for larger values of , as expected. Parameter plays a key role in the evolution of the disease spreading according to model (7) . If one was interested in modeling a realistic scenario, such as the spreading of COVID-19 in a ER waiting room, would have to be carefully tuned. and corresponding results on section = 0.5 (second row), evolution of ℎ for = 10 (third row) and corresponding results on section = 0.5 (fourth row). In both cases, the initial condition is IC1. Note that in all the subfigures time is dimensional while space is non-dimensional. The white dashed line in the images on row two and four represents * . The legend for all the images is the same as in Fig. 2 .

This first set of tests was meant to verify our implementation of method described in Sec. 3 and to check that the disease spreading term in eq. (7) (i.e., the third term on the left-hand side) produced the expected outcomes. Next, we will set people in motion.

We assign to all people walking speed = 1 m/s and walking direction = /4, as if they were headed to the upper right corner of the domain in the -plane. Once the people in the spreading phase of the disease have left the domain, we assume they cannot spread to the people in the walkable domain anymore.

We consider again IC1 and IC2 with = 50. Fig. 5 shows the evolution of the distribution density ℎ for both initial conditions. As one would expect, the motion contributes to lowering the exposure level in the both cases, since some of the spreading people leave the domain as soon as the motion starts. This difference is obvious when one compares the second (resp., fourth) row in Fig. 3 with the second (resp., fourth) row in Fig. 5 .

Finally, we modify the initial conditions to show that our model can handle also more complex scenarios featuring uncertainties. IC1 and IC2 are changed to: -IC1-bis: people are positioned like in IC1 but the probabilities of finding people with = 1 and = 0 is reduced and another value of is assigned at a given ( , ), i.e. at every location in the walkable domain the initial distribution has two bumps in . See Fig. 6 (left). -IC2-bis: people are positioned like in IC2 but the probabilities of finding people with = 1 and = 0 is reduced and, similarly to IC1-bis, another value of is assigned at a given ( , ). See Fig. 6 (right). Fig. 7 shows the evolution of the distribution density ℎ for initial conditions IC1-bis and IC2-bis with = 50.

In this section, we remove the simplifying assumption used in Sec. 2, i.e. walking speed and direction are no longer assumed to be given. Instead, they result from the interaction with the surrounding people. Following [1, 28, 30] , we account for the granular feature of the system (i.e., the fact that the distance between pedestrians can range from small to large) by discretizing . Instead of being continuous, variable can only take values in the set: Fig. 5 : Tests with prescribed walking velocity: the first row shows the evolution of the distribution density ℎ for initial condition IC1 and highlighted in green is the section whose results are reported in the second row. The third row shows the evolution of ℎ for initial condition IC2 and highlighted in green is the section whose results are reported in the fourth row. Note that in all the subfigures time is dimensional while space is non-dimensional. For both cases, = 50. The white dashed line in the images on row two and four represents * . The legend for all the images is the same as in Fig. 2 . where is the maximum number of possible directions. As for the walking speed , we assume it depends on the local level of congestion as demonstrated by experimental studies. Thus, variable is treated as a deterministic variable. Let = ( , ) denote position. For the above mentioned assumptions on and , the distribution function can be written as

where ( , ) represents the people that, at time and position , move with direction , and denotes the Dirac delta function. Like for the model in Sec. 2, it is convenient to switch to non-dimensional variables: positionˆ= / , walking speedˆ= / , timeˆ= / , and distribution functionˆ= / . Once again, we will drop the hat for ease of notation though, with the understanding that all variables from now on are non-dimensional. Due to the normalization of and the functions, the dimensionless local density is obtained by summing the distribution functions over the set of directions:

The model in [1, 28, 30] features two key parameters: Fig. 7 : Tests with prescribed walking velocity: the first row shows the evolution of the distribution density ℎ for initial condition IC1-bis and highlighted in green the section whose results are reported in the second row. The third row shows the evolution of ℎ for initial condition IC2-bis and highlighted in green the section whose results are reported in the fourth row. Note that in all the subfigures time is dimensional while space is non-dimensional. In both cases, we set = 50. The white dashed line in the images on row two and four represents * . The legend for all the images is the same as in Fig. 6 .

-∈ [0, 1], which represents the quality of the walkable domain: if = 1 the domain is clear, while if = 0 an obstruction is present. Note that could be dependent on space (and possibly time) in a prescribed way.

-∈ [0, 1], which represents the overall level of stress: = 0 indicates no stress, while = 1 indicates a high stress stress situation. Parameter plays a role in determining the walking speed: where = 1 people can walk at the desired speed (i.e., ), while where = 0 people are forced to slow down or stop. See, e.g., [28] on how to define the walking speed = [ ] ( , , ), with the square brackets denoting that depends on in a functional way. Parameter plays a role in the selection of the walking direction when interacting with other pedestrians. Indeed, the choice of the waking direction among the possible directions in results from four competing factors:

1. The goal to reach a target, like, e.g., the gate of an airport terminal. 2. The desire to avoid collisions with the walls and other obstacles. 3. The tendency to look for less crowded areas. 4. The tendency to follow the stream or "herd". Parameter weights between 3 (prevailing when stress is low) and 4 (prevailing when stress is high).

In order to explain how the four factors are modeled, we need to introduce some terminology. Interactions involve three types of active particles (or people):

test particle with state ( , ): they are representative of the whole system; candidate particle with state ( , ℎ ): they can reach in probability the state of the test particle after individual-based interactions with the venue (walls, obstacles, or target) or with field particles; field particle with state ( , ): their interactions with candidate particles determines a possible change of state. Note that if a test particle changes their state (in probability) into that of the test of the test particle as a result of the interactions with field people, the test particle loses their state.

We call A ℎ ( ) the probability that a candidate particle with direction ℎ adjusts their direction to (the direction of the test particle) due to the interaction with the venue. The transition probability that the candidate particle changes their direction to (direction of the test particle) in the search for less crowded areas if their stress level is low or to (direction of the field particle) if their stress level is high is denoted by B ℎ ( ). The sets of all transition probabilities A = {A ℎ ( )} ℎ, =1,..., and B = {B ℎ ( )} ℎ, , =1,..., form the so-called tables of games that model the game played by active particle interacting with the venue and other particles, respectively. See, e.g., [28] for more details.

In order to state the mathematical model, we need two more ingredients: -the interaction rate with the venue [ ] : it models the frequency of interactions between candidate particles and the venue and it is assume to have a functional dependence on . In fact, it is easier for pedestrians to see the walls, possible obstacles, and their target if density is low.

-the interaction rate with other pedestrians [ ]: it defines the number of binary encounters per unit time. This rate depends on in a functional way too since if the local density increases the interaction rate also increases. For possible ways to set [ ] and [ ], see, e.g., [28] .

In order to derive the mathematical model, one has to take the balance of particles in an elementary volume of the space of microscopic states, considering the net flow into such volume due to transport and interactions. Let J [ ] be the net balance of particles that move with direction :

due to interactions with the venue (first term at the left-hand side of (22)) and with the surrounding particles (second term at the left-hand side of (22)). Then, the model is given by:

for = 1, 2, . . . , . Problem (22)-(23) is completed with equation (21) for the density and a suitable law = [ ] ( , , ) that relates the walking speed to the density. While modeling motion in 2 dimensions ( and ), eq. (23) is a 3D problem in variables , , and . Its numerical approximation requires a carefully designed scheme to contain the computational cost. Here, we are only going to present ideas on how to design such a scheme. Numerical results will be presented in a follow-up paper. Following what we have done in [27] for emotional contagion, we will split model (23) into subproblems that are easier to solve using operator splitting (see, e.g., [22] ). For example, with the Lie splitting algorithm one could break problem (23) into: (i) a pure advection problem that features the contagion term, i.e. the last term in eq. (23), and (ii) a problem involving the interaction with the venue and other pedestrians. Then, suitable finite difference scheme will be applied for the space and time discretization of problems (i) and (ii).

We presented a way to design a kinetic type model to simulate the early stage of an infectious disease spreading in an intermediate size population occupying a confined environment for a short period of time.

In order to focus on the spreading mechanism, in the first part of the paper we adopted a strong simplifying assumption: people's walking speed and direction are given. The disease spreading is modeled using three main ingredients: an additional variable for the level of exposure to people spreading the disease, a parameter for the contagion interaction strength, and a kernel function that is a decreasing function of the distance between a person and a spreading individual. For such simplified model, we propose a first-order and a second-order finite difference scheme. We tested our numerical approach on simple 2D problems and verified that the level of exposure evolved as expected, even in scenarios featuring uncertainties.

In the second part of the paper, we used ideas from the simplified model to incorporate disease spreading in a kinetic theory approach for crowd dynamics. The appealing feature of kinetic (or mesoscopic) models for crowd dynamics is the flexibility in accounting for multiple interactions (hard to achieve in microscopic models) and heterogeneous behavior in people (hard to achieve in macroscopic models). Modeling crowd dynamics and disease contagion in two spatial dimensions requires the solution of a 3D problem, the additional variable being the contagion level. The design of an efficient numerical scheme for such problem will be object of a follow-up paper.

A kinetic theory approach to the dynamics of crowd evacuation from bounded domains

Biological systems as nonequilibrium structures described by kinetic methods

Microscopic pedestrian simulation model combined with a tactical model for route choice behaviour

A unified multiscale vision of behavioral crowds

On the modeling of crowd dynamics: Looking at the beautiful shapes of swarms

A quest towards a mathematical theory of living systems. Modeling and Simulation in Science, Engineering and Technology

From the microscale to collective crowd dynamics

A multiscale model of virus pandemic: Heterogeneous interactive entities in a globally connected world

What is life? a perspective of the mathematical kinetic theory of active particles

Toward a mathematical theory of behavioral-social dynamics for pedestrian crowds

Behavioral crowds: Modeling and Monte Carlo simulations toward validation

On the interplay between behavioral dynamics and social interactions in human crowds

Contagion shocks in one dimension

Convergence of probability measures

Modeling and simulating the spatial spread of an epidemic through multiscale kinetic transport equations

Modelling of Pedestrian and Evacuation Dynamics

Multiscale Modeling of Pedestrian Dynamics

Simulation of pedestrian counter flow through bottlenecks by using an agent-based model

Multiscale model for the optimal design of pedestrian queues to mitigate infectious disease spread

Individual-based markov model of virus diffusion: Comparison with covid-19 incubation period, serial interval and regional time series

Spread and dynamics of the covid-19 epidemic in italy: Effects of emergency containment measures

Finite element methods for incompressible viscous flow, in Handbook of numerical analysis

Agent-based simulation of pedestrian dynamics for exposure time estimation in epidemic risk assessment

Social force model for pedestrian dynamics

The mathematics of infectious diseases

The flow of human crowds

A kinetic theory approach for 2D crowd dynamics with emotional contagion

A kinetic theory approach to model pedestrian dynamics in bounded domains with obstacles

Coupling kinetic theory approaches for pedestrian dynamics and disease contagion in a confined environment

A kinetic theory approach to model crowd dynamics with disease contagion

Modeling the dynamics of coronavirus disease pandemic coupled with fear epidemics

A microscopic approach to study the onset of a highly infectious disease spreading

Empirical results for pedestrian dynamics and their implications for modeling

Optimization-based feedback control for pedestrian evacuation from an exit corridor

Understanding human queuing behaviour at exits: An empirical study

Efficient numerical methods for multiscale crowd dynamics with emotional contagion

Acknowledgements This work has been partially supported by NSF through grant DMS-1620384.