A Spatial Kinetic Model of Crowd Evacuation Dynamics with Infectious Disease Contagion

Abstract

This paper proposes a kinetic theory approach coupling together the modeling of crowd evacuation from a bounded domain with exit doors and infectious disease contagion. The spatial movement of individuals in the crowd is modeled by a proper description of the interactions with people in the crowd and the environment, including walls and exits. At the same time, interactions among healthy and infectious individuals may generate disease spreading if exposure time is long enough. Immunization of the population and individual awareness to contagion is considered as well. Interactions are modeled by tools of game theory, that let us propose the so-called tables of games that are introduced in the general kinetic equations. The proposed model is qualitatively studied and, through a series of case studies, we explore different scenarios related to crowding and gathering formation within indoor venues under the spread of a respiratory infectious disease, obtaining insights on specific policies to reduce contagion that may be implemented.

1 Introduction

During the recent COVID-19 crisis, it has become evident to the whole society that controlling and eradicating an epidemic from a population is a highly challenging task. Mathematical modeling became imperative to understand the epidemiological dynamics of disease spreading and to evaluate the impact of control measures implemented over time, with complex dynamics coming into play when trying to predict disease propagation. Used as a tool to guide public health authorities, epidemiological models were able to give insights on how to best use available strategies to contain and mitigate the ongoing outbreaks. However, there are still many aspects to explore regarding the development of mathematical models considering spatial features, with some recent contributions in Aguiar et al. (2020, 2021); Bellomo et al. (2020, 2022). In particular, a variety of social distancing measures has been observed all around the globe (see Qian and Jiang 2022; Thu et al. 2020; Sarcinschi 2020 among others). Most indoor venues have been closed during large periods of times and in most cases many restrictions regarding limitations in the number of people were implemented. It is worth stressing, however, that this concern is not new and not limited to COVID-19 contagion. Indeed, deal of attention has been paid to aerosol and airborne transmitted diseases in general, like the severe acute respiratory syndrome (SARS) outbreaks of 2003, the human avian influenza A (H5N1) infections and the more recent pandemic influenza A(H1N1/2009) cases in 2009 (Levy 2006; Tang et al. 2011; Teunis et al. 2010; Modchang et al. 2012), among others.

Most epidemiological models are based on the average of the behaviors of a large population over a given period of time. In particular, compartmental models (starting from the celebrated model by Kermack and McKendrick (1927)) use mean-field approximations which are attractive because of their simplicity. However, these models involve complex parameters that depend on many factors, which makes it difficult to predict how a change in a single environmental, demographic or epidemiological condition will affect the whole population (Brauer et al. 2019). Moreover, these models are not valid if the population size is small-to-medium, as happens in some spatial domains (neighborhoods, stations, schools, etc.) that are very relevant in the dynamics around the development of an epidemic (Aguiar et al. 2021). Consequently, when it comes to assess particular responses to specific outbreaks, studying how the proximity of people plays a role in the diffusion of a disease and what can be done in crowded areas and mass gatherings could inform more targeted responses. In this present paper, a SEIR (Susceptible-Exposed-Infectious-Recovered) approach is considered to model disease spreading in a general situation. However, it will be shown that, due to its versatility, the epidemiological model can be tuned to consider specific dynamics, e.g., saturation of the infection force (Kolokolnikov and Iron 2021), immunization of the population (Qiu and Feng 2010) or isolation (Zuzek et al. 2015).

On the other hand, crowd dynamics has so far been studied to tackle safety problems (Bellomo et al. 2019; Degond et al. 2017), for instance by supporting crisis management in critical situations such as sudden evacuations induced in complex venues by incidents. Modeling of crowd dynamics requires a multiscale approach, since it can be developed at three scales, namely microscopic (individual-based), macroscopic (hydrodynamic), and mesoscopic (kinetic). The mesoscopic representation, which is typical of the mathematical kinetic theory, is delivered by a probability distribution function over the state of the individual entities. As observed in the recent book (Aylaj et al. 2020) (see also Liao and Zhou 2022 for the movement of aggregated groups and Bellomo et al. (2022) for the modeling of dynamical emotional states), the selection of the scale is a key problem as each one presents advantages and withdraws.

For instance, the description of the dynamics at the micro-scale is delivered by a system of ordinary differential equations corresponding to a Newtonian type dynamics. The key problem consists in modeling the acceleration terms as well as of the overall crowd distribution in space. Hence, individual-based models focus on interactions between each individual with the others. One of the main contributions is the so-called social force model proposed by Helbing and coworkers (Helbing and Molnar 1995), which has been further developed to model evacuation dynamics (Helbing et al. 2002) and panic situations (Helbing et al. 2000). In principle, this approach can be generalized to the case of groups moving according to different behaviors, as an example, by including the presence of leaders, as well as multiple interactions. As stated in Aylaj et al. (2020), accounting for the heterogeneity of walkers, which is an important feature of crowd dynamics, still needs additional work to be exhaustively treated. On the other hand, the macroscopic scale requires the continuity assumption of the matter which is reasonable in the case of real fluids, but not in the case of crowds, where distances between individuals cannot be neglected. In addition, local averages hide the heterogeneous behaviour of individuals. Also, models derived by the kinetic theory approach require the continuity assumption of the probability distribution which is not realistic in the case of flows with a limited number of individuals (Cercignani et al. 1994). However, the main problem to take into account in the modeling at all scales is the heterogeneous behavioral feature of people in a crowd and its specific influence on interactions of walkers to be interpreted as active, rather than classical, particles (Helbing and Johansson 2009; Haghani and Sarvi 2017). Individual behaviors and heterogeneity, in the case of crisis circumstances, can lead to significant deviations with respect to the usual dynamics in rational flow conditions.

Crowd and epidemiological modeling have been so far treated as separated fields of research, with very few attempts to link them together (Kim and Quaini 2020; Terna et al. 2020). In this paper, we study the propagation of an epidemic through the development of mathematical models of crowd dynamics based on the kinetic theory of active particles. Our departure point is the model proposed in Agnelli et al. (2015), where the evacuation of a crowd from a bounded domain with exit doors is developed. The model considers mechanical transitions only, accounting for changes in the velocity of pedestrians (speed and direction) according to a decision-based strategy: pedestrians may decide to move towards exit areas and avoid the walls (dynamics induced by the geometry of the domain) while following the mainstream or trying to avoid large concentrations (dynamics induced by interactions among pedestrians). In this present paper, we assume that pedestrians are also carriers of a new micro-state variable describing their state related to an infectious disease: namely, they can belong at a given time t to one compartment (e.g., susceptible, exposed, infectious, recovered). Thus, we will couple the kinetic evacuation model with epidemiological dynamics occurring within the same domain. The complete model will be introduced and then used to analyze the interplay between evacuation and contagion, in order to understand how much proximity, residence time before exiting, and even vaccination influence on the overall spread of the disease.

The paper is structured as follows: Section 2 resumes the kinetic model of crowd evacuation and couples it with the contagion dynamics, by considering a microscopic variable accounting for the disease-related state of particles. Section 3 presents numerical results in which the model is tested through three exploratory case studies. In particular, we perform a parameter sensitivity analysis and explore the role of the infection rate, immunization, and population awareness on the overall dynamics. Finally, in Section 4 we conclude the paper and present some possible further extensions.

2 A Mathematical Model

2.1 Functional Subsystems and Representation

Let us consider a large system of pedestrians moving in a bounded environment $\Omega \subset \mathbb {R}^2$. We assume that the boundary $\partial \Omega$ includes the exit zone $E\subset \partial \Omega$, while the remaining part of the boundary constitutes the wall $W\subset \partial \Omega$. It is worth noticing that E could be the finite union of disjoint sets, i.e., the domain may have several exit zones. We consider, for simplicity, a convex domain, as shown in Fig. 1. Even though the presence of internal obstacles is not included in the following treatment, it can be done through a straightforward technical generalization of the model, see (Kim and Quaini 2019). The domain $\Omega$ is also characterized by a parameter $\alpha \in [0,1]$ accounting for the quality of the environment, where $\alpha = 0$ corresponds to the worst quality which forces pedestrians to stop, while the value $\alpha = 1$ corresponds to the best quality, that contributes keeping high speeds. We can think about $\alpha$ as an indicator of some features like signaling, lighting, presence of obstacles or fire, etc.

Fig. 1

Geometry of the bounded domain $\Omega$ with boundary $\partial \Omega = W \cup E_1 \cup E_2$

Each pedestrian is seen as an active particle carrier of a microscopic state. Following the ideas presented in Agnelli et al. (2015), we consider a representation with continuous-discrete hybrid features, where the micro-state of each pedestrian is defined by position and velocity direction and disease-related state at time t. More in detail:

• The position $\textbf{x}= (x, y)$ is supposed to be a continuous variable defined over $\Omega$.

• For the velocity, denoted by $\textbf{v}= v(\cos \theta ,\sin \theta )$ in polar coordinates, it is assumed that the speed v is a continuous deterministic variable which evolves in time and space according to macroscopic effects determined by the overall dynamics, while the velocity direction ${\theta }$ is a discrete variable -heterogeneously distributed among pedestrians- attaining values in the set

  $$\begin{aligned} I_{{\theta }}= \left\{ {\theta }_i =\frac{i-1}{n} 2\pi :\,i=1, \ldots , N_{\theta }\right\} , \end{aligned}$$

  with cardinal $N_{\theta }$.

• A disease-related state given by a categorical variable in the set $I_\lambda$ with cardinal $N_\lambda$.

Remark 1

The set $I_\lambda$ will be adjusted according to the case under study. For instance, if we consider four compartments for susceptibles (S), exposed (E), infected (I), and recovered (R), then $I_\lambda =\{S,E,I,R\}$ and $N_\lambda =4$. If a vaccinated class (V) is considered, then $I_\lambda =\{S,E,I,R,V\}$ and $N_\lambda =5$. This way of introducing heterogeneity gives flexibility to choose the most suitable epidemiological model.

Remark 2

The term functional subsystem (FS) will be used to identify groups of particles sharing a common micro-state. For instance, the i-FS refers to those active particles moving with direction ${\theta }_i$ or the E-FS refers to those active particles belonging to the exposed class E.

The overall state of particles is given by the distribution function $f_i^j(t,\textbf{x})=f(t,\textbf{x},{\theta }_i, \lambda _j)$, which is interpreted as the number of particles that, at time t, are located in $\textbf{x}$, move with direction ${\theta }_i$ and are carriers of the disease-related state $\lambda _j$. The knowledge of the distribution function lets us obtain several macroscopic quantities of interest through the computation of its moments. For instance, the local density in $\textbf{x}$ is given by

$$\begin{aligned} \nonumber \rho (t,\textbf{x}) = \sum _{j\in I_\lambda }\sum _{i=1}^{N_{\theta }} f_{i}^{j}(t,\textbf{x}), \end{aligned}$$

the local density of individuals belonging to the $\lambda _j$-class is

$$\begin{aligned} \nonumber \rho ^j(t,\textbf{x}) = \sum _{i=1}^{N_{\theta }} f_{i}^{j}(t,\textbf{x}), \end{aligned}$$

and the total population within the domain is given by

$$\begin{aligned} \nonumber N(t) = \int _\Omega \rho (t,\textbf{x}^*)d\textbf{x}^*. \end{aligned}$$

In this way, the number of individuals belonging to each infectious compartment within the domain $\Omega$ can be computed as follows:

$$\begin{aligned} \begin{aligned} S(t)&= \int _\Omega \rho ^S(t,\textbf{x}^*)d\textbf{x}^*, \quad E(t) = \int _\Omega \rho ^E(t,\textbf{x}^*)d\textbf{x}^*, \\ I(t)&= \int _\Omega \rho ^I(t,\textbf{x}^*)d\textbf{x}^*, \quad R(t) = \int _\Omega \rho ^R(t,\textbf{x}^*)d\textbf{x}^*, \end{aligned} \end{aligned}$$

where S(t), E(t), I(t) and R(t) are the susceptible, exposed, infectious and recovered populations within the domain, respectively. Notice that these quantities may change over time due to the effect of the two coupled dynamics under study, namely contagion dynamics and exit (resp. entrance) of individuals from (resp. to) the domain.

Remark 3

In general, physical dimensions are removed through nondimensionalization, in such a way that the maximum reachable local density is equal to 1 under normal conditions. We refer to Agnelli et al. (2015) and Section 3 for more details.

2.2 Mathematical Structure

The modeling approach takes into account interactions of pedestrians with all other pedestrians and with the environment. The derivation of the mathematical structure can be obtained by a suitable balance of particles in the elementary volume of the space of microscopic states, considering the net flow into such volume due to transport and interactions

$$\begin{aligned} \partial _{t}f_i^j(t,\textbf{x}) + \textrm{div}_\textbf{x}(\textbf{v}_i[\rho ](t,\textbf{x})f_i^j(t,\textbf{x}) ) = \mathcal{J}_i^j [f](t,\textbf{x}), \;\;j\in I_\lambda , \;i=1,\dots ,N_{\theta }, \end{aligned}$$

(1)

where $\textbf{v}_i[\rho ] = v[\rho ](\cos {\theta }_i,\sin {\theta }_i)$, the left-hand term models the transport of particles, while $\mathcal{J}_i^j[f]$ represents the net balance for those particles in the $\lambda _j$-FS that move with direction ${\theta }_i$ due to interactions.

As in Agnelli et al. (2015), to model interaction dynamics, we suppose that pedestrians modify their walking direction by a decision strategy which takes into account the following trends:

• 1. Moving towards the exit;

• 2. Avoiding collision with walls;

• 3. Moving towards less congested areas;

• 4. Attraction to follow the main stream.

The first two types of dynamics are related to purely geometric aspects of the domain, while the last two take into consideration that pedestrians’ behavior is strongly affected by the presence of other people. In this sense, items 3 and 4 are related to another interaction which modifies the internal state of pedestrians, namely: In this sense, items 3 and 4 are related to another kind of interaction which modifies the internal state of pedestrians that is now considered in this work, namely:

• 5. Interaction between susceptible and infectious pedestrians may produce contagion and spread of the infectious disease.

Hence, we split the interaction term into two terms

$$\begin{aligned} \mathcal{J}_{i}^{j}[f] = \mathcal{J}^G_i[f] + \mathcal{J}_{i}^{j P}[f], \end{aligned}$$

where $\mathcal{J}^G_i$ is the difference between the gain and the loss of particles moving with direction ${\theta }_i$, due to geometrical effects (notice that it is independent of the disease-related state) and $\mathcal{J}_{i}^{j P}$ accounts for the balance due to interactions among particles belonging to the $\lambda _j$ and i-FSs with the others subsystems.

Taking inspiration from Agnelli et al. (2015), we model changes in velocity direction by using probabilistic rules, through the so-called transition probabilities $\mathcal{A}_h(i)$ for $h,i=1,\dots ,N_{\theta }$ and $\mathcal{B}^1_{hk}(i)$ for $h,k,i=1,\dots ,N_{\theta }$.

The quantity $\mathcal{A}_h(i)$ is the probability that an h-particle, i.e., a pedestrian moving with direction ${\theta }_h$, adjusts his/her direction into ${\theta }_i$ as a consequence of the domain geometry (e.g., walls or exit doors), while $\mathcal{B}^1_{hk}(i)$ is the probability that an individual moving with direction ${\theta }_h$ changes his/her direction into ${\theta }_i$ after an interaction with an individual walking with direction ${\theta }_k$.

In addition, in order to account for disease-related transitions, we introduce the transition probability density $\mathcal{B}^2_{sr}(j)$ giving the probability that a $\lambda _s$-individual may undergo a transition into the state $\lambda _j$ as a consequence of an interaction with a $\lambda _r$-individual.

Remark 4

Notice that the symbol $\mathcal{A}$ is used to describe transitions involving only the geometry of the domain, while $\mathcal{B}$ is used for inter-pedestrian interactions.

The compact form of the transition probability $\mathcal{B}$ is given by the product

$$\begin{aligned} \mathcal{B}_{hk}^{sr}(i,j)[\rho ] = \mathcal{B}^1_{hk}(i)[\rho ] \times \mathcal{B}^2_{sr}(j), \end{aligned}$$

where $\mathcal{B}_{hk}^{sr}(i,j)[\rho ]$ is the probability that a pedestrian with infectious state $\lambda _s$ moving with direction $\theta _h$ undergoes a transition into the infectious state $\lambda _j$ and direction $\theta _i$ after an interaction with a pedestrian with state $\lambda _r$ moving with direction $\theta _k$.

Remark 5

The transition probabilities $\mathcal{A}$ and $\mathcal{B}$ satisfy

$$\begin{aligned} \sum _{i=1}^{N_{\theta }}\mathcal{A}_h(i)=1,{} & {} \forall \, h = 1, \ldots , N_{\theta }, \nonumber \\ \sum _{j=1}^{N_\lambda } \sum _{i=1}^{N_{\theta }}\mathcal{B}_{hk}^{sr}(i,j)=1,{} & {} \, \forall \, h,\!k, = 1, \ldots , N_{\theta }, \,\, \forall \, s,\!r \in I_\lambda . \end{aligned}$$

(2)

By considering the above-defined terms, Eq. (1) can be written as

$$\begin{aligned}{} & {} \underbrace{\partial _{t} f_i^j(t,\textbf{x}) + \textrm{div}_\textbf{x}(\textbf{v}_i[\rho ](t,\textbf{x})f_i^j(t,\textbf{x}))}_{\text{ transport }} = \underbrace{ \mathcal{J}_i^j [f](t,\textbf{x})}_{\text{ interactions }} \nonumber \\{} & {} \quad = \mu [\rho (t,\textbf{x})] \left( \sum _{h=1}^{N_{\theta }} \mathcal{A}_{h}(i) f_h^j(t,\textbf{x}) - f_i^j(t,\textbf{x}) \right) \nonumber \\{} & {} \qquad +\eta [\rho (t,\textbf{x})]\left( \sum _{s,r=1}^{N_\lambda } \sum _{h,k=1}^{N_{\theta }} \mathcal{B}_{hk}^{sr}(i,j)[\rho ] f_h^s(t,\textbf{x}) f_k^r(t,\textbf{x}) - f_i^j(t,\textbf{x}) \rho (t,\textbf{x})\right) \qquad \end{aligned}$$

(3)

for $i=1,\dots ,N_{\theta }$, $j\in I_\lambda$, and where $\textbf{v}_i[\rho ] = v[\rho ](\cos {\theta }_i,\sin {\theta }_i)$. The terms $\mu$ and $\eta$ are the interaction rates modeling the frequency of interactions with the geometry and with other pedestrians, respectively. Moreover, conditions (2) guarantee the conservation in the number of total pedestrians.

This mathematical structure offers a general framework suitable to derive specific models. It is more general than the model proposed in Agnelli et al. (2015) since it includes the presence of different groups featured by different walking abilities and strategies to develop the dynamics.

2.3 Modeling Interactions

Let us now focus on the modeling of interactions that results in the specification of the right-hand side of Eq. (3). Interactions involve three types of particles (Bellomo et al. 2021): with micro-state $(\textbf{x},\theta _i,\lambda _j)$ which are representative of the whole system; with micro-state $(\textbf{x},\theta _k,\lambda _r)$, whose presence triggers the interactions of the candidate particles; and with micro-state $(\textbf{x},\theta _h,\lambda _s)$, which can reach in probability the state of the test particles after individual-based interactions with field particles or with the environment. In what follows we refer to an i-particle to mean a pedestrian moving with direction ${\theta }_i$.

Two types of interactions are considered: (i) those between candidate and field particles and (ii) those between candidate particles and the environment where the dynamics occurs, that is with its geometrical and qualitative properties.

2.3.1 Towards the Selection of the Desired Direction

We introduce the dynamics that generates changes in the direction of movement of active particles and describe how the decision process, which consists in selecting the desired direction for each pedestrian, is modeled.

Fig. 2

a We denote the distance to the exit of a particle located in $\textbf{x}$ by $d_E(\textbf{x})$ and the vector pointing from $\textbf{x}$ to the exit by ${\vec {\nu }}(\textbf{x})$. b A particle in $\textbf{x}$ moving with direction ${\theta }_h$ is expected to collide the wall in $\textbf{x}_W$, then it computes the tangent direction to the wall that would take it toward the exit

We assume that not only binary interactions induce a change in the state of a particle, but pedestrians make a decision according to a combination of different causes, as detailed below.

• Trend to move toward the exit. During an evacuation, pedestrians may try to reach the exit by moving through the shortest path from their current location. Given a candidate particle at the point $\textbf{x}$, we define its distance to the exit as

  $$\begin{aligned} d_E(\textbf{x}) = \min _{\textbf{y}\in E} \Vert \textbf{x}- \textbf{y}\Vert , \end{aligned}$$

  (4)

  where $\Vert \cdot \Vert$ denotes the Euclidian norm in $\mathbb {R}^2$, and we consider the unitary vector ${\vec {\nu }}(\textbf{x})$, pointing from $\textbf{x}$ to the exit, see Fig. 2(a).

• Trend to avoid the collision with walls. If a particle at $\textbf{x}$ moves with direction ${\theta }_h$, and if this direction does not point it towards the exit, then it is expected to collide (unless it changes direction) with the wall at a point $\textbf{x}_W(\textbf{x},{\theta }_h)$, which is located at a distance $d_W(\textbf{x},{\theta }_h)$ from the particle, as shown in Fig. 2(b). Then, the particle should choose a suitable direction capable to prevent it hitting the wall after some time. Following the model in Agnelli et al. (2015), the unitary tangent vector ${\vec {\tau }}(\textbf{x}, {\theta }_h)$ to $\partial \Omega$ at $\textbf{x}_W$ pointing in the direction that would let a pedestrian get closer to the exit is chosen.

• Tendency to move towards less congested areas. In order to facilitate its movement, a candidate particle at $\textbf{x}$, moving with direction ${\theta }_h$, may decide to change direction by moving towards less congested areas. This can be achieved by choosing the direction that gives the minimal directional derivative of the density at the point $\textbf{x}$. This direction is denoted by the unitary vector ${\vec {\gamma }}({\theta }_h,\rho )$.

• Tendency to follow the stream. Binary interactions, at each time t and position $\textbf{x}$, involve test, candidate, and field particles. A candidate particle modifies its state, in probability, into that of the test particle, due to interactions with field particles, while the test one loses its state as a result of these interactions. This dynamics is inserted in the model in order to take into account the fact that a (candidate) pedestrian interacting with a (field) pedestrian may decide to follow him/her. We define the unitary vector ${\vec {\sigma }}_k = (\cos {\theta }_k, \sin {\theta }_k)$ to describe the movement of the field k-particle.

Let us observe that the first two effects are related to purely geometric aspects of the domain, meaning that candidate particles take into account the presence of doors or walls but without caring about other people’s behavior. Conversely, the last two effects take into consideration that people’s behavior is strongly affected by that of the others. In fact, on the one hand a candidate particle is capable to scan its surroundings in order to choose, at each moment and position, a proper direction that will prevent it to move into congested areas, while on the other hand the interaction with a field particle will try to bring its direction closer to that of the latter.

The modeling approach is based on the following assumptions:

• (A1) The trend to the exit increases as particles get closer to it.

• (A2) Particles are subject to a stronger influence to avoid the wall as they get closer to it.

• (A3) The tendency to look for less congested areas depends on the local density.

• (A4) The tendency to follow the stream depends on the local density.

Notice that the effects related to assumptions (A3)-(A4) compete with each other. In other words, higher densities will induce a higher tendency to look for less congested areas but at the same time to follow the stream. We introduce a parameter $\varepsilon \in [0,1]$ that reinforces one effect or the other according to the particular situation to be modeled. The value $\varepsilon =0$ corresponds to the situation in which only the research of less congested areas is considered, while $\varepsilon =1$ corresponds to the situation in which only the tendency to follow the stream is taken into account.

2.3.2 Interaction Terms and Selection of the Desired Direction

Pedestrians change their velocity direction according to a not purely deterministic decision which is taken by considering all the above introduced effects. This feature can be efficiently modeled in a probabilistic manner, since different pedestrians are not expected to react in the same way when facing a certain particular situation. More precisely, at each interaction, each pedestrian is assumed to play a game whose payoff is the updating of his/her direction accordingly to his/her strategy.

Dynamics induced by the shape of the environment.

This type of dynamics is modeled by means of the following two interaction terms.

• The interaction rate $\mu [\rho ]$ models the frequency of interactions between candidate particles and the boundary of the domain. We suppose that $\mu$ decreases with local density, since the lower this quantity is, the easier is for pedestrians to realize about the presence of walls and doors. Following this idea, we assume $\mu [\rho ]\sim 1-\rho$.

• The transition probability $\mathcal{A}_{h}(i)$ is the probability that an h-candidate particle adjusts its direction into that of the test particle ${\theta }_i$, induced by the presence of walls and exit areas.

The modeling approach assumes that particles change direction, in probability, only to an adjacent (clockwise or anticlockwise) direction in the discrete set $I_{\theta }$. This means that a candidate h-particle may eventually end up into the states $h-1$, $h+1$ or remain in the state h. Notice that in the case $h=1$ we set $\theta _{h-1}={\theta }_{N_{\theta }}$, while in the case $h={N_{\theta }}$ we set ${\theta }_{h+1}={\theta }_1$.

The set of all transition probabilities $\mathcal{A}=\{\mathcal{A}_{h}(i)\}_{h,i=1,\dots ,N_{\theta }}$ forms the so-called table of games that models the game played by active particles interacting with the geometry of the environment.

According to assumptions (A1)-(A2), we define the vector

$$\begin{aligned} {\vec {\omega }}_G(\textbf{x},{\theta }_h) = \left[ 1 - d_E(\textbf{x}) \right] \,{\vec {\nu }}(\textbf{x}) + \left[ 1 - d_W(\textbf{x}, {\theta }_h) \right] {\vec {\tau }} \, (\textbf{x}, {\theta }_h), \end{aligned}$$

whose direction ${\theta }_G$ is the , meaning the ideal direction that a pedestrian should take in order to reach the exit and avoid the walls in an optimal way.

Since directions are discretized, an h-particle will update its direction by choosing among the three allowed outputs ${\theta }_{h-1}$, ${\theta }_h$ and ${\theta }_{h+1}$ the closest to ${\theta }_G$. The compact form of the table of games ${\mathcal {A}}$ is given by

$$\begin{aligned} {\mathcal {A}}_h(i) = \beta _h(\alpha ) \delta _{s,i} + \left[ 1 - \beta _h(\alpha ) \right] \delta _{h,i}, \quad i=h-1, h, h+1, \end{aligned}$$

where

$$\begin{aligned} s= & {} {{\,\mathrm{arg\,min}\,}}_{j\in \{h-1,h+1\}} \{d({\theta }_G, {\theta }_{j})\}, \\ d({\theta }_*,{\theta }^*)= & {} \left\{ \begin{array}{lcl} |\theta _* - \theta ^* |, &{} \quad \hbox {if} &{} |{\theta }_*-{\theta }^* |\le \pi , \\ &{} &{} \\ 2\pi - |{\theta }_*-{\theta }^* |, &{}\quad \hbox {if} &{} |{\theta }_*-{\theta }^* |> \pi , \end{array} \right. \end{aligned}$$

$\delta _{j,i}$ denotes the Kronecker delta function, and the coefficient $\beta _h$, proportional to the parameter $\alpha$, is introduced to describe the fact that even in the case that the geometrical preferred direction ${\theta }_G$ is close to ${\theta }_h$, a transition may occur, and the more distant the two states are, the more probable is this transition:

$$\begin{aligned} \beta _h(\alpha ) = \left\{ \begin{array}{lcl} \alpha , &{} \quad \hbox {if} &{} d({\theta }_h,{\theta }_G)\ge \Delta {\theta }, \\ &{} &{} \\ \alpha \frac{d({\theta }_h,{\theta }_G)}{\Delta {\theta }}, &{}\quad \hbox {if} &{} d({\theta }_h,{\theta }_G) < \Delta {\theta }, \end{array} \right. \end{aligned}$$

where $\Delta {\theta }=2\pi /n$. Notice that if ${\theta }_G = {\theta }_h$, then $\beta _h=0$ and ${\mathcal {A}}_h(h)=1$, meaning that a pedestrian keeps the same direction, at least due to the geometrical effects, with probability 1.

Dynamics induced by interactions between pedestrians.

In this case, the interaction dynamics can be modeled as follows:

• The interaction rate $\eta [\rho ]$, defines the number of binary encounters per unit time. We assume that the interaction rate increases with increasing local density. As proposed in Bellomo et al. (2013), we assume $\eta [\rho ] = \eta _0 \rho ,$ where $\eta _0$ is the rate, to be measured by experimental data, corresponding to the spatially homogeneous case at low densities.

• The transition probability $\mathcal{B}^1_{hk}(i)[\rho ]$ is the probability that a candidate particle modifies its direction into that of the test particle ${\theta }_i$ induced by a decision that combines the research of less congested areas and the interaction with a field particle moving with direction ${\theta }_k$. The square brackets denote the dependence on the density $\rho$.

Concerning the so-called vacuum direction ${\vec {\gamma }}$, we have to consider how pedestrians react according to their perception of the density around them, and the game consists in the choice of the less congested direction among the three admissible ones. This direction can be computed for a candidate h-pedestrian in position $\textbf{x}$ by taking

$$\begin{aligned} \ell = \ell ({\theta }_h,\rho (t,\textbf{x}))= {{\,\mathrm{arg\,min}\,}}_{j\in \{h-1,h,h+1\}} \{\partial _j\rho (t,\textbf{x})\}, \end{aligned}$$

where $\partial _j\rho$ denotes the derivative of $\rho$ in the direction ${\theta }_j$. In this way, we have ${\vec {\gamma }}({\theta }_h,\rho )=(\cos \theta _\ell ,\sin \theta _\ell )$.

According to assumptions (A3)-(A4) and recalling that ${\vec {\sigma }_k}$ denotes the vector pointing in the direction of the field particle, we define the vector

$$\begin{aligned} {\vec {\omega }}_P({\theta }_h,{\theta }_k,\rho ) = \varepsilon \,{\vec {\sigma }_k} + (1 - \varepsilon ) {\vec {\gamma }} ({\theta }_h,\rho ), \end{aligned}$$

where the subscript P stands for pedestrians, and the direction ${\theta }_P$ of ${\vec {\omega }}_P$ is the , obtained as a weighted combination between the trend to follow the stream and the tendency to avoid crowded zones. Then, we propose the following table of games:

$$\begin{aligned} \mathcal{B}^1_{hk}(i)[\rho ] = \beta _{hk}(\alpha ) \rho \delta _{r,i} + \left[ 1 - \beta _{hk}(\alpha ) \rho \right] \delta _{h,i}, \quad i=h-1, h, h+1, \end{aligned}$$

where r and $\beta _{hk}$ are defined as in the previous case

$$\begin{aligned} r= & {} {{\,\mathrm{arg\,min}\,}}_{j\in \{h-1,h+1\}} \{d({\theta }_P,{\theta }_j)\}, \\ \beta _{hk}(\alpha )= & {} \left\{ \begin{array}{lcl} \alpha , &{} \quad \hbox {if} &{} d({\theta }_h,{\theta }_P)\ge \Delta {\theta }, \\ &{} &{} \\ \alpha \frac{d({\theta }_h,{\theta }_P)}{\Delta {\theta }}, &{}\quad \hbox {if} &{} d({\theta }_h,{\theta }_P) < \Delta {\theta }. \end{array} \right. \end{aligned}$$

2.4 Modeling the Contagion Dynamics

Each individual is carrier of an internal micro-state accounting for their condition related to an infectious disease. Let us first consider the simplest case in which the population is partitioned into four mutually exclusive compartments (that can be treated as well as functional subsystems in our theoretical framework) with $I_\lambda =\{S,E,I,R\}$. Here, S, E, I and R denote, respectively, susceptible, exposed, infectious and recovered hosts.

Since the time scale of interest for movement and evacuation from a room is too short (order of minutes) compared to the period from time of infection to time of being contagious or infectious (order of days), it is imperative to include in our model the disease latency period and, thus, the exposed compartment E. When a susceptible individual interacts with an infectious one, she may get infected but -during the latency period- will not transmit the disease, until a transition from E to I takes place. This may occur for sure after the evacuation and we are interested in monitoring the dynamics of contagion during evacuation.

The dynamics of an SEIR model can be described by the following reaction scheme:

$$\begin{aligned} \begin{array}{ccc} S + I &{} \xrightarrow {i_{r_S}} &{} E + I \\ E &{} \xrightarrow {p_{EI}} &{} I \\ I &{} \xrightarrow {p_{IR}} &{} R, \end{array} \end{aligned}$$

(5)

where permanent waning immunity is being assumed (thus, individuals belonging to the R compartment remain there).

Here, we assume that contagion occurs following the law of mass action, with infection rate $i_{r_S}$ at which a susceptible individual may become exposed after an interaction with an infected individual. The transitions from E to I and from I to R depend only on the compartment sizes, with $p_{EI}$ and $p_{IR}$ denoting the respective transition rates.

As mentioned above, since we are only interested in the evacuation time interval, our table of games shall only consider the first transition in Eq. (5). Thus, rates $p_{EI}$ and $p_{IR}$ can be neglected and the only non-trivial entry in the table of games for contagion dynamics takes the form $\mathcal{B}^2_{SI}(E)=i_{r_S}$, while all the other interactions do not undergo any kind of transition.

It is useful to include a vaccinated class, as will be considered in the numerical experiments afterwards. To do so, we simply add a compartment V including susceptible individuals that can get infected but with a lower infection rate. In the case of a perfect sterilizing vaccine, the infection rate would be 0 but in general we shall assume a non-negative value $i_{r_V}\le i_{r_S}$.

The reaction scheme for the so-called SEIRV model is:

$$\begin{aligned} \begin{array}{ccc} S + I &{} \xrightarrow {i_{r_S}} &{} E + I \\ V + I &{} \xrightarrow {i_{r_V}} &{} E + I \\ E &{} \xrightarrow {p_{EI}} &{} I \\ I &{} \xrightarrow {p_{IR}} &{} R, \end{array} \end{aligned}$$

(6)

and the transition probability densities shall be changed accordingly by adding the non-trivial entry $\mathcal{B}^2_{VI}(E)=i_{r_V}$.

Remark 6

It is worth noticing that the addition of a V-FS can model two possible scenarios:

$\bullet$ Firstly, the straightforward interpretation of having immunized particles that can contract the disease with a rate $i_{r_V}\le i_{r_S}$.

$\bullet$ On the other hand, V can also be thought as a class of individuals with a larger awareness level, namely following recommendations and taking cautions to avoid contagion (e.g., wearing masks). That is translated in a lower probability of getting infected, with $i_{r_V}\le i_{r_S}$.

2.5 Modeling the Velocity Modulus

The decay or increase of the velocity modulus depends on the interactions between pedestrians. It is assumed that people adjust their speed according to the level of congestion around them. We assume that the maximal reachable dimensionless velocity modulus $v_m=v_m(\alpha )$ depends linearly on the quality of the environment, in such a way that $v_m(0)=0$—any movement is hindered—and $v_m(1)=1$—the maximal speed can be reached. Here, we assume that the speed depends on the local perceived density. In particular, the highest reachable speed $v_m(\alpha )$ is kept under low density conditions (free flow regime), up to a certain critical density $\rho _c(\alpha )$ that can be experimentally measured. For values of $\rho$ greater than $\rho _c$, the velocity modulus decreases to zero (slowdown zone). In the slowdown zone, we choose a polynomial-type dependence of the velocity modulus on the local density, see Fig. 3. We refer to Sec. 2.3 of Agnelli et al. (2015) for more details.

Fig. 3

Dependence of the dimensionless velocity modulus v on the dimensionless density $\rho$ for different values of the parameter $\alpha$ representing the quality of the environment. In the free flow zone ($\rho \le \rho _c(\alpha )=\alpha /5$) pedestrians move with the maximal velocity modulus $v_m(\alpha )=\alpha$ allowed by the environment. In the slowdown zone ($\rho >\rho _c(\alpha )$) pedestrians have a velocity modulus which is here heuristically modeled by the third order polynomial joining the points $(\rho _c(\alpha ),v_m(\alpha ))$ and (1, 0) and having horizontal tangent in such points (Color figure online)

3 Numerical Results and Case Studies

This section presents some numerical simulations of the above introduced model, in order to capture the overall dynamics of disease spreading under a variety of scenarios.

As a starting point, we endow system (3) with initial conditions $f_i^j(0,\textbf{x})$, for $i=1,\ldots ,N_{\theta }$ and $j\in I_\lambda = \{S, E, I, R, V\}$. Following the reasonings in Agnelli et al. (2015), boundary conditions are not explicitly imposed, but are induced by the non-local action over the particles given by the term ${\mathcal {J}}_i^j$. Indeed, individuals feel the presence of walls at a distance and modify their dynamics in order to avoid them (see Sect. 2.3.1 for more details). The numerical solution of the system is obtained by using a splitting method (Holden et al. 2010), where the overall evolution operator is seen as the sum of evolution operators for each term in the model. Then, each term is solved by means of an appropriate scheme, and finally all the pieces are attached together. In particular, Eq. (1) is splitted into two subequations:

$$\begin{aligned} \partial _{t}f_i^j(t,\textbf{x}) + \textrm{div}_\textbf{x}(\textbf{v}_i[\rho ](t,\textbf{x})f_i^j(t,\textbf{x}) ) = 0, \end{aligned}$$

(7)

and

$$\begin{aligned} \partial _{t}f_i^j(t,\textbf{x}) = \mathcal{J}_i^j [f](t,\textbf{x}), \end{aligned}$$

(8)

for $j\in I_\lambda , \;i=1,\dots ,N_{\theta }$. Equation (7) is a 2D homogeneous transport equation that we solve using a finite volume scheme. This method guarantees conservation of the total number of particles and we refer to LeVeque (2002); Piccoli and Tosin (2011); Schäfer (2006) for more details on the implementation. On the other hand, Eq. (8) is solved by means of a first order Euler explicit method to go forward in time.

In the following we consider a square domain of side length 10 m with an exit door of width 2 m. This is the same domain used in Agnelli et al. (2015); Kim and Quaini (2019) and it is useful to analyze how evacuation and epidemiological dynamics are merged in the model. Details on how to compute the vectors pointing to the exits or letting pedestrian avoid the walls can be found in those references. The set $I_{\theta }$ is defined as a discrete size of $n=8$ directions:

$$\begin{aligned} I_{\theta }= \Big \{{\theta }_i=\frac{i-1}{8}2\pi , i = 1, \ldots , 8 \Big \}, \end{aligned}$$

while the velocity modulus is assumed to depend on the perceived density and on the quality of the environment, as described in Agnelli et al. (2015). We will deal with non-dimensional quantities, that are obtained by referring the spatial coordinates relative to the longest dimension of the domain $L = 10\sqrt{2}$ m, velocity modulus to $V_M = 2$ m/s and density to $\rho _M = 7$ $m^{-2}$ as in Buchmueller and Weidmann (2006). From these values, we get the reference time $T_M = 5\sqrt{2}$ s.

Figure 4 shows some snapshots of the evacuation process. We start with around 50 pedestrians initially moving with direction ${\theta }_1$. The proportion of individuals belonging to S, I and V compartments are, respectively, 60%, 25% and 15%, while none of them is initially exposed. In Fig 4 (d) we see how the proportion of exposed individuals evolves as the evacuation takes place. Notice that contagion is considered only within the domain as we do not monitor disease dynamics once individuals are outside the room. For this case, we took $\alpha = 1$.

Fig. 4

Figures a, b and c show three snapshots of the evacuation process and d the evolution of exposed population E(t) vs t. For this experiment, the population is initially distributed in the following proportions: 60% in S, 25 % in I, 15% in V and 0% in E (Color figure online)

In the following we present three case studies aiming to analyze the model performance under different scenarios. The first case study is devoted to the sensitivity to parameters $\alpha$ and $\varepsilon$, the second case study deals with the role of the infection rate, while the third one refers to the role of immunization and/or contagion awareness.

3.1 Case-Study 1: Sensitivity to $\alpha$ and $\varepsilon$

In Agnelli et al. (2015) a great deal of attention was given to understand the dependence of the evacuation time on the quality of the environment $\alpha$ and the parameter weighing “stream vs vacuum” effects $\varepsilon$. Now, since we are dealing with the coupled contagion model, we want to understand the influence of these two parameters on the spread of an infectious disease. Consequently, we study the evacuation and contagion dynamics as we move in the $(\alpha , \varepsilon )$ space, while keeping the infection rates $i_{r_S} = 0.1$ and $i_{r_V} = 0.01$ fixed. In particular, as shown in Fig. 5a it is clear that poor environment conditions lead to a tougher evacuation, and thus larger evacuation times. As observed in Fig. 5b for fixed $\alpha$ the variation with respect to $\varepsilon$ is not much significant.

Fig. 5

a Evacuation time as a function of $\alpha$ for different values of $\varepsilon$. b Evacuation time as a function of $\varepsilon$ for different values of $\alpha$ (Color figure online)

Figure 6 shows the final proportion of exposed pedestrians (when the evacuation is over) as a function of (a) $\alpha$ for different values of $\varepsilon$ and (b) $\varepsilon$ for different values of $\alpha$. We see that larger $\varepsilon$ values generate higher levels of contagion, since this is correlated with a stronger tendency to follow the stream and, consequently, larger gatherings: pedestrians remain closer and the transmission of a pathogen is facilitated.

On the other hand, we see that larger values of $\alpha$ are related, in general, to smaller fractions of exposed individuals, due to shorter exposure times (as observed above). However, the relationship is not always monotonic as for the evacuation time, since the role of $\varepsilon$ and crowding is acting as well.

Fig. 6

Final proportion of exposed pedestrians (when the evacuation is over) as a function of a $\alpha$ for different values of $\varepsilon$ and b $\varepsilon$ for different values of $\alpha$. One can observe that larger $\varepsilon$ values generate higher levels of contagion. On the other hand, larger values of $\alpha$ are related, in general, to smaller fractions of exposed individuals, due to shorter exposure times. However, the relationship is not always monotonic due since the role of $\varepsilon$ and crowding is acting as well (Color figure online)

Finally, Fig. 7 shows the combined sensitivity analysis on $\alpha$ and $\varepsilon$, where the above described phenomena is jointly observed.

Fig. 7

Final proportion of exposed pedestrians (when the evacuation is over). Each dot shows the evacuation time (horizontal axis) and the proportion of exposed individuals at exit time (vertical axis) for a given combination of $\alpha$ and $\varepsilon$. Notice that experiments with the same $\alpha$ value form clusters of points illustrated with different colors (moving to the right for decreasing $\alpha$’s). For each fixed $\alpha$, the evacuation time is roughly the same, and points move up in the plane (higher levels of contagion) as $\varepsilon$ increases (Color figure online)

These qualitative results are in agreement with the empirical evidence of a strong exposure time-risk relationships for airborne transmitted diseases. For instance, knowing the exposure time of healthy subjects in closed venues, is a useful tool to estimate the corresponding individual infection risk, as quantitatively discussed in Buonanno et al. (2020) for several scenarios in the case of SARS-CoV-2. In addition, in Costanzo and Flores (2022) authors present a contagion risk estimator model based on Wells-Riley probabilistic approach for interior spaces in which people share the same so-called “airborne shared space”. Starting from a single infectious individual in the simulated environments, results show a positive correlation between probability of infection and permanence time.

3.2 Case-study 2: On the Role of the Infection Rate

The infection rate measures the risk of infection upon interaction between infected and susceptible individuals. As a key parameter appearing in most epidemiological models, highly related to outbreak formation, appearance of endemic states and extinction of a disease (Brauer et al. 2019), it is important to understand its role within the context of our evacuation model. This case-study is devoted to a qualitative understanding of the role of the infection rate on the overall contagion dynamics.

To do so, let us now vary the infection rates $i_{r_S}$ and $i_{r_V}$ for non-vaccinated and vaccinated (or aware) pedestrians, respectively, while keeping fixed $\alpha =1$ and $\varepsilon =0.4$.

Initially, we consider that $75\%$ of the population is initially healthy, with $60\%$ belonging to the S compartment and $15\%$ the V compartment.

Figure 8 shows the exposed population after the evacuation is over for different combinations of the infection rates. It can be seen that increasing infection rates lead to larger proportion of exposed individuals, as expected. In particular, by comparing Figs. 8a and b we see that contagion is much more sensitive to $i_{r_S}$ than to $i_{r_V}$ in a predominantly susceptible population. This is especially relevant when designing social distancing and containment measures for the transmission of a respiratory disease within a closed domain. In order to keep the spread of the disease under certain levels, and particularly when large fractions of population are not immunized -as we will discuss in the next case-study- non-pharmaceutical interventions to keep $i_{r_S}$ and $i_{r_V}$ under control must be implemented.

Fig. 8

Case-study 2: On the role of the infection rate. Final proportion of exposed pedestrians (when the evacuation is over) as a function of a $i_{r_S}$ for different values of $i_{r_V}$ and b $i_{r_V}$ for different values of $i_{r_S}$. It can be observed that increasing infection rates lead to larger proportion of exposed individuals. Additionally, note that contagion is much more sensitive to $i_{r_S}$ than to $i_{r_V}$ (Color figure online)

3.3 Case-study 3: The Role of Immunization

In this last exploratory scenario, we study the role of immunization (or awareness) programs. In particular, the aim is to qualitatively understand if larger proportions of individuals with a lower probability of contracting the disease lead to a reduction in contagion.

In the following numerical experiments, $\alpha =1$ and $\varepsilon =0.4$ are still fixed, but now the infection rates $i_{r_S} = 0.1$ and $i_{r_V} = 0.01$ are kept fixed as well. The proportion of infected pedestrians I goes from $20\%$ to $70\%$ (thus $80\%$ to $30\%$ for initially healthy individuals). Among healthy individuals, we change the ratio between vaccinated or aware V(0) and non-vaccinated susceptible S(0) pedestrians.

In the following, we analyze the influence of the fraction of immunized (or aware) individuals on the overall dynamics.

Figure 9 shows the fraction of exposed individuals relative to initially healthy pedestrians after evacuation for different initial proportions of whether vaccinated individuals or individuals with larger levels of awareness. It is clear from this experiment that both immunizing the population and/or reaching larger proportion of individual awareness have a strong effect on reducing the spread of the disease. Two important conclusions arise from this experiment. First, notice that larger proportions of infected individuals lead to a higher level of contagion at the end of the evacuation. Second, as we move to the right in the figure (namely, increasing $V(0)/[V(0)+S(0)]$) the lines become closer. This supports the fact that reaching high levels of immunity or awareness helps to keep contagion under control independently of the force of infection in the range under study.

This result confirms the evidence that, for several diseases, application of appropriate measures (in this case, immunization or suitable awareness controls such as social distancing, mask wearing or hand hygiene) aimed at enlarging the proportion V related to $V+S$ is of major importance to reduce the likelihood of transmission and thereby protect the population (Morawska et al. 2020).

Fig. 9

Case-study 3: On the role of the immunization. Final proportion of exposed pedestrians when the evacuation is over relative to initially healthy population for different proportions of whether vaccinated individuals or individual awareness levels. Each line corresponds to a different fraction of infected pedestrians (Color figure online)

4 Conclusions

An extension of a kinetic model for crowd evacuation from bounded domains presented in Agnelli et al. (2015) to deal with the spread of an infectious respiratory disease within indoor venues has been presented in this work. Some contributions dealing with individual awareness to contagion were discussed and incorporated in Kim and Quaini (2020), while in Bellomo et al. (2020, 2022); Aguiar et al. (2021) treated specifically the epidemiological problem of contagion. The main contribution of this present paper is to merge both approaches, by considering that pedestrians moving within a domain can be described not only by their position and velocity but also by a microscopic variable accounting for their state regarding the disease (susceptible, exposed, infectious, recovered, vaccinated).

Even if a huge amount of epidemiological modeling peer reviewed articles appeared since the COVID-19 pandemic started (Else 2020), the model is suitable to explore any airborne transmitted infectious disease that needs proximity among individuals to be spread (Tang et al. 2011), e.g., influenza A (Teunis et al. 2010; Modchang et al. 2012), tuberculosis (Deol et al. 2022) or pneumococcal pneumonia (Torén et al. 2022). For all of them, inhalation of droplets, from either persons with the disease or healthy carriers, constitute a major transmission route.

As a result, we obtained a rich model that is highly versatile to explore several scenarios and situations. For example, the V functional subsystem might represent both the presence of an immunized class or of a group of individuals with larger awareness levels. Thus, through simple modifications of the transition probabilities the model is suitable to deal with a variety of different assumptions both from the epidemiological and the evacuation points of view.

The numerical experiments showed that the model is robust and let us explore several features of the dynamics while moving the parameter space or changing initial conditions. The proposed case-studies summarize situations that need to be tackled to propose specific policies to reduce contagion in closed venues. In particular, the considered exploratory scenarios show a strong agreement with empirical evidence related to exposure time and risk awareness. Very useful tools have been recently developed and confirm our results that, as time indoors increases, infections can occur in shared room air despite distancing (Peng et al. 2022; Peng and Jimenez 2021; Miller et al. 2021). Our results qualitatively agree with observed contagion dynamics in the case of SARS-CoV2, as discussed in Srikrishna (2020) in the case of closed room or a classroom with comparable number of individuals.