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

Daewa Kim  Department of Mathematics, West Virginia University,
94 Beechurst Ave, Morgantown, WV 26506
<EMAIL_ADDRESS>Kaylie O’Connell  Department of Mathematics,
University of Houston,
3551 Cullen Blvd, Houston TX 77204
<EMAIL_ADDRESS>William Ott  Department of Mathematics, University of
Houston,
3551 Cullen Blvd, Houston TX 77204
<EMAIL_ADDRESS>Annalisa Quaini  Department of Mathematics, University of
Houston,
3551 Cullen Blvd, Houston TX 77204
<EMAIL_ADDRESS>
###### Abstract

In this paper, we present a computational modeling approach for the dynamics
of human crowds, where the spreading of an emotion (specifically fear) has an
influence on the pedestrians’ behavior. Our approach is based on the methods
of the kinetic theory of active particles. The model allows us to weight
between two competing behaviors depending on fear level: the search for less
congested areas and the tendency to follow the stream unconsciously (herding).
The fear level of each pedestrian influences her walking speed and is
influenced by the fear levels of her neighbors. Numerically, we solve our
pedestrian model with emotional contagion using an operator splitting scheme.
We simulate evacuation scenarios involving two groups of interacting
pedestrians to assess how domain geometry and the details of fear propagation
impact evacuation dynamics. Further, we reproduce the evacuation dynamics of
an experimental study involving distressed ants.

###### keywords:

Crowd dynamics; kinetic model; complex systems; fear propagation; evacuation

AMS Subject Classification: 35Q91, 65M06, 91C99, 92D99

##  1 Introduction

Due to the complexity of human behavior, the mathematical modeling of crowd
motion has attracted a lot of attention. While there exists a large body of
literature on this topic, the interest in designing more and more realistic
models is not declining. This is due to the fact that realistic modeling of
human crowd dynamics could bring great benefit to society. In fact, the
quantitative and computational analysis associated with a mathematical model
could support e.g. crisis managers in emergency situations. In this paper, we
focus on one important aspect of crowd behavior: the effect on crowd dynamics
of an emotional contagion (specifically fear) that spreads within the crowd
itself.

In broad terms, the models of crowd dynamics that have been proposed can be
grouped into three categories by scale. One category is the microscopic,
individual-based model inspired by Newtonian mechanics focusing on inter-
particle interaction. For a survey on the _microscopic_ approach, we refer to
Ref. RevModPhys.73.1067. The second category treats crowd dynamics as a
continuum flow, an approach that works well for large-scale, dense crowds. For
a survey on this _macroscopic_ approach, see e.g. Refs. Hughes2003 and
Cristiani2014_book.

The third category, called _mesoscopic_ or _kinetic_ , uses a scale of
observation in between the previous two. For introductory concepts, we refer
to Refs. Bellomo2013_new and BellomoBellouquid2015. Further developed in Refs.
Bellomo2017_book,Bellomo2015_new,Bellomo2019_new and kim_quaini, the mesoscale
approach derives a Boltzmann-type evolution equation for the statistical
distribution function of the position and velocity of the pedestrians, in a
framework close to that of the kinetic theory of gases. There is one key
difference though: the interactions in Boltzmann equations are conservative
and reversible, while the interactions in the kinetic theory of active
particles are irreversible, non-conservative and, in some cases, nonlocal and
nonlinearly additive. An important consequence is that often for active
particles the Maxwellian equilibrium does not exist. [2]

Multiscale approaches have also been investigated. See e.g. Ref. Bellomo2020
for a multiscale view of human crowds that is consistent at the three possible
modeling scales. Each modeling scale has its advantages and pitfalls. We
choose to work with a kinetic-type model because of its flexible ability to
account for multiple interactions and heterogeneous behavior among
pedestrians.

By and large, existing mathematical models assume that pedestrians behave
rationally. However, in emergency situations people may behave irrationally
when responding to fear. A majority of the current models therefore do not yet
effectively support crowd managers. Several works have shown that social
phenomena such as fear transmission are pervasive in crowds and modify the
interaction rules. [19, 20, 21, 27, 30, 35, 36] One of the first works on a
multiscale (from microscopic to macroscopic) approach to crowd dynamics with
emotional contagion is Ref. Bertozzi2015. Therein, fear is propagated by a
BGK-like model (see e.g. Ref. Cercignani1994_book) and results are limited to
one space dimension. Studies on the impact of social dynamics on individual
interactions and their influence at a higher scale are carried out in Refs.
Degond2013,Degond2017. More recently, a kinetic approach to modeling
pedestrian dynamics in the presence of social phenomena (e.g. propagation of
stress) is presented in Ref. Bellomo2019_new. The numerical results in Ref.
Bellomo2019_new show that stress propagation significantly affects crowd
density patterns and overall crowd dynamics.

In this paper, we extend the kinetic-type model in Ref. kim_quaini to account
for the propagation of stress conditions in time and space. To this end, we
draw inspiration from the 1D model in Ref. Bertozzi2015. The results produced
by the model in Ref. kim_quaini have been compared against data collected in a
recent empirical study reported in Ref. Kemloh. We have shown that for medium-
sized groups of people, the average population density and flow rate computed
with our kinetic-type model agree with the respective measured quantities in
the absence of fear. The model in Ref. Bertozzi2015 adapts the model in Refs.
Bertozzi2014ContagionSI,688c8c70f9c748e486b48806e7ce2fdb and has been designed
to track the level of fear within individuals under the assumption that fear
influences motion. The numerical results produced by the agent-based model in
Ref. 688c8c70f9c748e486b48806e7ce2fdb have compared more favorably to actual
crowd footage than other pedestrian models. See Ref. inproceedings.

We assess our kinetic model with emotional contagion through three test cases.
Our first two tests involve two groups of pedestrians with different initial
fear distributions. We have set up an evacuation scenario wherein the two
groups interact as they move toward the exit. We study how the evacuation
dynamics depend on domain geometry (specifically the exit size) and on the
spatial length scale over which fear can be directly transmitted. Our results
for these two tests confirm what observed in Ref. Bellomo2019_new: the
spreading of a stress condition significantly affects the overall evacuation
dynamics and leads to overcrowding in certain areas of the domain. Our third
test demonstrates that our model can recapitulate the results of experimental
ant panic studies.

The paper is organized as follows. In Sec. 2, we present a 2D extension of the
model in Ref. Bertozzi2015 under one simplifying assumption: pedestrians
change walking speed because of the spreading of fear, but they do not alter
walking direction. We discretize the simplified model and present numerical
results. In Sec. 3, we remove the simplifying assumption and present our full
2D model, where fear affects both walking speed and direction. Sec. 4
describes a splitting algorithm for the model in Sec. 3 and Sec. 5 presents
our numerical results for the three test cases mentioned above. We draw
conclusions in Sec. 6.

##  2 A simplified two-dimensional kinetic model

We extend the one-dimensional emotional contagion model in Ref. Bertozzi2015
to 2D with one important simplifying assumption: the pedestrians have a common
walking direction $\theta$ that does not change over time. This assumption
will be subsequently removed (see Section 3).

We start from an agent-based model at the microscopic level. Let particle $n$
represent a pedestrian, whose position and fear level are given by
${\bm{x}}_{n}(t)=(x_{n}(t),y_{n}(t))^{T}$ and $q_{n}(t)$, respectively. In
addition, let $q_{n}^{*}$ be an average fear level perceived by pedestrian
$n$. The microscopic model reads:

$\displaystyle\frac{d{\bm{x}}_{n}}{dt}=q_{n}(\cos\theta,\sin\theta)^{T},\quad\frac{dq_{n}}{dt}=\gamma(q_{n}^{*}-q_{n}),$
$\displaystyle
q_{n}^{*}=\frac{\sum_{m=1}^{N}\kappa_{n,m}q_{m}}{\sum_{m=1}^{N}\kappa_{n{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0},}m}},\quad
n=1,2,3,\dots,N,$ (1)

where $N$ is the total number of particles. This model assumes that: (i) at a
given time an agent walks faster if more scared (see the first equation in (1)
where the walking speed is equal to the fear level) and (ii) agents tend to
equilibrate their fear level with the average fear level they perceive. In
model (1), the interaction kernel $\kappa_{n,m}$ serves as the weights in the
average $q_{n}^{\ast}$. This kernel is a decreasing function of mutual
distance between two particles and is parametrized by an interaction distance
$R$:

$\displaystyle\kappa_{n,m}=\kappa(|{\bm{x}}_{n}-{\bm{x}}_{m}|)=\dfrac{R}{(|{\bm{x}}_{n}-{\bm{x}}_{m}|^{2}+R^{2})\pi}.$
(2)

Parameter $\gamma$ in (1) describes the contagion interaction strength and it
may vary from particle to particle for more general cases.

From the agent-based model (1), we derive a model at the kinetic level. Denote
the empirical distribution by

$h^{N}=\frac{1}{N}\sum_{n=1}^{N}\delta({\bm{x}}-{\bm{x}}_{n}(t))\delta(q-q_{n}(t)),\\\
$

where $\delta$ is the Dirac delta measure. We assume that the particles remain
in a fixed compact domain
$({\bm{x}}_{n}(t),q_{n}(t))\in\Omega\subset\mathbb{R}^{3}$ for all $n$ and for
the entire time interval under consideration. Prohorov’s theorem implies that
the sequence $\\{h^{N}\\}$ is relatively compact in the weak∗ sense.
Therefore, there exists a subsequence $\\{h^{N_{k}}\\}_{k}$ such that
$h^{N_{k}}$ converges to $h$ with weak∗-convergence in
$\mathcal{P}(\mathbb{R}^{3}$) and pointwise convergence in time as
$k\rightarrow\infty$. Here, $\mathcal{P}(\mathbb{R}^{3}$) denotes the space of
probability measures on $\mathbb{R}^{3}$.

Consider a test function $\psi\in C_{0}^{1}(\mathbb{R}^{3})$. We have

$\displaystyle\frac{d}{dt}\langle h^{N},\psi\rangle_{{\bm{x}},q}$
$\displaystyle=\frac{d}{dt}\biggl{<}\frac{1}{N}\sum_{n=1}^{N}\delta({\bm{x}}-{\bm{x}}_{n}(t))\delta(q-q_{n}(t)),\psi\biggr{>}_{{\bm{x}},q}$
$\displaystyle=\frac{d}{dt}\frac{1}{N}\sum_{n=1}^{N}\psi({\bm{x}}_{n}(t),q_{n}(t))$
$\displaystyle=\frac{1}{N}\sum_{n=1}^{N}\psi_{x}q_{n}\cos\theta+\psi_{y}q_{n}\sin\theta+\psi_{q}\gamma(q_{n}^{\ast}-q_{n})$
$\displaystyle=\langle\psi_{x}q\cos\theta,h^{N}\rangle_{{\bm{x}},q}+\langle\psi_{y}q\sin\theta,h^{N}\rangle_{{\bm{x}},q}$
$\displaystyle\quad+\frac{\gamma}{N}\sum_{n=1}^{N}\psi_{q}\biggl{(}\frac{\sum_{m=1}^{N}\kappa_{n,m}q_{n}}{\sum_{m=1}^{N}\kappa_{n,m}}-q_{n}\biggr{)},$
(3)

where $\langle\cdot\rangle_{{\bm{x}},q}$ means integration against both
${\bm{x}}$ and $q$ and the subindex of $\psi$ denotes the variable with
respect to which we take the partial derivative of $\psi$.

Let us define

$\rho_{h^{N}}({\bm{x}})=\frac{1}{N}\sum_{n=1}^{N}\delta({\bm{x}}-{\bm{x}}_{n})\\\
$

and

$m_{h^{N}}({\bm{x}})=\biggl{<}q,\frac{1}{N}\sum_{m=1}^{N}\delta({\bm{x}}-{\bm{x}}_{m})\delta(q-q_{m})\biggr{>}_{{\bm{x}},q}=\frac{1}{N}\sum_{m=1}^{N}\delta({\bm{x}}-{\bm{x}}_{m})q_{m},\\\
$

We have

$\frac{1}{N}\sum_{m=1}^{N}\kappa(|{\bm{x}}_{n}-{\bm{x}}_{m}|)=\biggl{<}\kappa(|{\bm{x}}_{n}-\tilde{{\bm{x}}}|),\frac{1}{N}\sum_{m=1}^{N}\delta(\tilde{{\bm{x}}}-{\bm{x}}_{m})\biggr{>}_{{\bm{x}}}=\kappa\ast\rho_{h^{N}}({\bm{x}}_{n}),\\\
$
$\frac{1}{N}\sum_{m=1}^{N}\kappa(|{\bm{x}}_{n}-{\bm{x}}_{m}|)q_{m}=\biggl{<}\kappa(|{\bm{x}}_{n}-\tilde{{\bm{x}}}|),\frac{1}{N}\sum_{m=1}^{N}\delta(\tilde{{\bm{x}}}-{\bm{x}}_{m})q_{m}\biggr{>}_{{\bm{x}}}=\kappa\ast
m_{h^{N}}({\bm{x}}_{n}),\\\ $

where $\langle\cdot\rangle_{\bm{x}}$ means integration only in ${\bm{x}}$.
Then, we can rewrite eq. (2) as

$\frac{d}{dt}\langle
h^{N},\psi\rangle_{{\bm{x}},q}=\langle\psi_{x}q\cos\theta,h^{N}\rangle_{{\bm{x}},q}+\langle\psi_{y}q\sin\theta,h^{N}\rangle_{{\bm{x}},q}+\gamma\biggl{<}h^{N},\frac{\kappa\ast
m_{h^{N}}}{\kappa\ast\rho_{h^{N}}}\psi_{q}-q\psi_{q}\biggr{>}_{{\bm{x}},q},\\\
$

which leads to

$h_{t}^{N}+\nabla\cdot(q(\cos\theta,\sin\theta)^{T}h^{N})=\gamma((q-q^{\ast})h^{N})_{q},$
(4)

via integration by parts.

Now letting $k\rightarrow\infty$, the subsequence $h^{N_{k}}$ formally leads
to the limiting kinetic equation

$h_{t}+\nabla\cdot(q(\cos\theta,\sin\theta)^{T}h)=\gamma((q-q^{\ast})h)_{q},$
(5)

where $h(t,{\bm{x}},q)$ is the probability of finding a person in an
infinitesimal neighborhood of $({\bm{x}},q)$. The quantity $q^{*}(t,{\bm{x}})$
is the local _average_ emotional contagion level weighted by the distance to
${\bm{x}}$:

$q^{\ast}(t,{\bm{x}})=\frac{\iint\kappa(|{\bm{x}}-\overline{{\bm{x}}}|)h(t,\overline{{\bm{x}}},q)qdqd\overline{{\bm{x}}}}{\iint\kappa(|{\bm{x}}-\overline{{\bm{x}}}|)h(t,\overline{{\bm{x}}},q)dqd\overline{{\bm{x}}}}.$
(6)

###  2.1 Full discretization

In this section, we present the numerical discretization of the kinetic
equation (5).

Let us start from the space discretization. Divide the spatial and velocity
domain into a number of cells
$[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]\times[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}]$
of the length $\Delta x$ and $\Delta y$. The discrete mesh points $x_{i}$ and
$y_{j}$ are given by

$\displaystyle x_{i}=i\Delta x,\quad x_{i+1/2}=x_{i}+\frac{\Delta
x}{2}=\Big{(}i+\frac{1}{2}\Big{)}\Delta x,$ $\displaystyle y_{j}=j\Delta
y,\quad y_{j+1/2}=y_{j}+\frac{\Delta y}{2}=\Big{(}j+\frac{1}{2}\Big{)}\Delta
y,$ (7)

for $i=0,1,\dots,N_{x}$ and $j=0,1,\dots,N_{y}$, The contagion level domain is
partitioned into subdomains $[q_{k-\frac{1}{2}},q_{k+\frac{1}{2}}]$ and the
length of $\Delta q$ with $k\in 1,2,\dots,N_{q}$.

$\displaystyle q_{k}=k\Delta q,\quad q_{k+1/2}=q_{k}+\frac{\Delta
q}{2}=\Big{(}k+\frac{1}{2}\Big{)}\Delta q.$ (8)

Here $N_{x}$, $N_{y}$ and $N_{q}$ are the total number of points in $x-$, $y-$
and $q-$ directions, respectively. For simplicity, we assume each cell is
centered at $x_{i}$, $y_{j}$ or $q_{k}$ with a uniform length $\Delta x$,
$\Delta y$ and $\Delta q$.

Let us denote $h_{i,j,k}=h(t,x_{i},y_{j},q_{k})$. The average fear level is
computed using a midpoint rule for the integrals in (6):

$\displaystyle q^{\ast}(t,x_{i},y_{j})\approx q^{\ast}_{i,j}$
$\displaystyle=\frac{\sum_{k}\sum_{\overline{j}}\sum_{\overline{i}}\kappa_{\overline{i},\overline{j}}h_{\overline{i},\overline{j},k}q_{k}\Delta
q\Delta x\Delta
y}{\sum_{k}\sum_{\overline{j}}\sum_{\overline{i}}\kappa_{\overline{i},\overline{j}}h_{\overline{i},\overline{j},k}\Delta
q\Delta x\Delta y},$ (9) $\displaystyle\kappa_{\overline{i},\overline{j}}$
$\displaystyle=\frac{R}{((x_{i}-x_{\overline{i}})^{2}+(y_{j}-y_{\overline{j}})^{2}+R^{2})\pi}.$

We consider the following scheme for eq. (5):

$\displaystyle\partial_{t}h_{i,j,k}$
$\displaystyle+\frac{\overline{\eta}_{i,j,k}-\overline{\eta}_{i-1,j,k}}{\Delta
x}+\frac{\overline{\phi}_{i,j,k}-\overline{\phi}_{i,j-1,k}}{\Delta y}$
$\displaystyle+\gamma\frac{\xi_{i,j,k+\frac{1}{2}}-\xi_{i,j,k-\frac{1}{2}}}{\Delta
q}+\gamma\frac{C_{i,j,k+\frac{1}{2}}-C_{i,j,k-\frac{1}{2}}}{\Delta q}=0,$ (10)

where

$\displaystyle\overline{\eta}_{i,j,k}$
$\displaystyle=\eta_{i,j,k}+\frac{\eta_{i+1,j,k}-\eta_{i,j,k}}{2}\varphi\left(\frac{\eta_{i,j,k}-\eta_{i-1,j,k}}{\eta_{i+1,j,k}-\eta_{i,j,k}}\right),~{}\eta_{i,j,k}=q_{k}\cos\theta\,h_{i,j,k},$
(11) $\displaystyle\overline{\phi}_{i,j,k}$
$\displaystyle=\phi_{i,j,k}+\frac{\phi_{i,j+1,k}-\phi_{i,j,k}}{2}\varphi\left(\frac{\phi_{i,j,k}-\phi_{i,j-1,k}}{\phi_{i,j+1,k}-\phi_{i,j,k}}\right),~{}\phi_{i,j,k}=q_{k}\sin\theta\,h_{i,j,k},$
(12) $\displaystyle\xi_{i,j,k+\frac{1}{2}}$
$\displaystyle=\dfrac{|q_{i,j}^{\ast}-q_{k+\frac{1}{2}}|+(q_{i,j}^{\ast}-q_{k+\frac{1}{2}})}{2}h_{i,j,k}+\dfrac{(q_{i,j}^{\ast}-q_{k+\frac{1}{2}})-|q_{i,j}^{\ast}-q_{k+\frac{1}{2}}|}{2}h_{i,j,k+1}$
$\displaystyle C_{i,j,k+\frac{1}{2}}$
$\displaystyle=\dfrac{1}{2}\left|s_{i,j,k+\frac{1}{2}}\right|\left(1-\dfrac{\Delta
t}{\Delta
q}\left|s_{i,j,k+\frac{1}{2}}\right|\right)W_{i,j,k-\frac{1}{2}}\varphi\left(\dfrac{W_{i,j,\textbf{k}-\frac{1}{2}}}{W_{i,j,k-\frac{1}{2}}}\right)$
(13)

with

$\displaystyle q_{k+\frac{1}{2}}=\frac{q_{k}+q_{k+1}}{2},\quad
s_{i,j,k+\frac{1}{2}}=q^{\ast}_{i,j}-q_{k+\frac{1}{2}},\quad
W_{i,j,k-\frac{1}{2}}=h_{i,j,k}-h_{i,j,k-\frac{1}{2}}.$ (14)

In eq. (11) and (12), $\varphi$ is a slope limiter function. We choose the Van
Leer function:

$\varphi(\theta)=\frac{|\theta|+\theta}{1+|\theta|}.$ (15)

In eq. (13), the subscript k is k-1 if $s_{i,j,k-\frac{1}{2}}>0$ and k+1 if
$s_{i,j,k-\frac{1}{2}}<0$. This scheme is second-order in velocity thanks to
the use of the flux limiter.

Next, let us discretize in time. Let $\Delta t\in\mathbb{R}$,
$t^{n}=t_{0}+n\Delta t$, with $n=0,...,N_{T}$ and $T=t_{0}+N_{T}\Delta t$,
where $T$ is the end of the time interval under consideration. Moreover, we
denote by $y^{n}$ the approximation of a generic quantity $y$ at the time
$t^{n}$. We adopt the forward Euler scheme for the time discretization of
problem (2.1)

$\displaystyle h_{i,j,k}^{l+1}=h_{i,j,k}^{l}-\Delta t\Big{(}$
$\displaystyle\dfrac{\overline{\eta}^{l}_{i,j,k}-\overline{\eta}^{l}_{i-1,j,k}}{\Delta
x}+\dfrac{\overline{\phi}^{l}_{i,j,k}-\overline{\phi}^{l}_{i,j-1,k}}{\Delta
y}$ $\displaystyle+$
$\displaystyle\gamma\dfrac{\xi^{l}_{i,j,k+\frac{1}{2}}-\xi^{l}_{i,j,k-\frac{1}{2}}}{\Delta
q}+\gamma\dfrac{C^{l}_{i,j,k+\frac{1}{2}}-C^{l}_{i,j,k-\frac{1}{2}}}{\Delta
q}\Big{)}$ (16)

The time step $\Delta t$ is chosen as

$\displaystyle\Delta t=\frac{1}{2}\min\Bigg{\\{}\frac{\Delta
x}{\max_{k}q_{k}},\frac{\Delta y}{\max_{k}q_{k}},\frac{\Delta
q}{2\gamma\max_{k}q_{k}}\Bigg{\\}}$ (17)

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

###  2.2 Numerical results

The computational domain in the $xy$-plane is $[-10,10]\times[-10,10]$ and the
fear level $q$ is in $[0,1]$. We choose $\theta=\pi/4$ as walking direction.
For the space discretization, we choose $\Delta x=\Delta y=0.05$ (we note that
these are dimensional values) and $\Delta q=0.005$. We set $R=0.0002$ and
$\gamma=100$. These relatively small interaction radii and rather large
interaction strength (i.e. quick interactions) are meant to model a dense
crowd setting. The time interval under consideration is $[0,1]$ and the time
step $\Delta t$ is set to $1.25\cdot 10^{-5}$.

We approximate the delta function as follows:

$\delta(q)\sim
E(q)=\cfrac{1}{\sqrt{\pi}R_{0}}e^{-\cfrac{q^{2}}{R_{0}^{2}}},\quad
R_{0}=0.04.$ (18)

Following Ref. Bertozzi2015, the initial condition for the distribution
density is set as:

$\displaystyle h(0,{\bm{x}},q)$
$\displaystyle=\dfrac{1}{2}\left(\sin\left(\cfrac{\pi
x}{10}\right)+2\right)\left(\cfrac{1}{4}\,\,E(q-q_{I}(x)-0.5)+\cfrac{3}{4}\,\,E(q-q_{I}(x)+0.3)\right)$
$\displaystyle+\dfrac{1}{2}\ \left(\sin\left(\cfrac{\pi
y}{10}\right)+2\right)\left(\cfrac{1}{4}\,\,E(q-q_{I}(y)-0.5)+\cfrac{3}{4}\,\,E(q-q_{I}(y)+0.3)\right),$

where

$q_{I}(0,x)=\cfrac{1}{2}\Big{(}3-\tanh x\Big{)}.$

Figure 1 shows the simulated evolution of $h(t,{\bm{x}},q)$. Initially, $h$
has two bumps in $q$ for every ${\bm{x}}$, as shown in the panels on the left.
As time evolves, $h$ starts to concentrate on $q^{*}(t,{\bm{x}})$, as
confirmed by the panels on the right.

\begin{overpic}[width=138.76157pt,grid=false]{./contagion_contagion_time_0.pdf}\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./contagion_contagion_time_1000.pdf}\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./contagion_contagion_time_4000.pdf}\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./contagion_contagion_time_0_transparency.pdf}
\put(47.4,35.6){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){33.7}}
\put(47.2,69.2){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){6.9}}
\put(40.3,8.3){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){33.5}}
\put(47.4,35.6){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){6.8}}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./contagion_contagion_time_1000_transparency.pdf}
\put(47.4,35.6){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){33.7}}
\put(47.2,69.2){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){6.9}}
\put(40.3,8.3){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){33.5}}
\put(47.4,35.6){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){6.8}}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./contagion_contagion_time_4000_transparency.pdf}
\put(47.4,35.6){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){33.7}}
\put(47.2,69.2){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){6.9}}
\put(40.3,8.3){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(0,1){33.5}}
\put(47.4,35.6){\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\line(-1,-4){6.8}}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./contagion_contagion_time_0_diagonal.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./contagion_contagion_time_1000_diagonal.pdf}\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./contagion_contagion_time_4000_diagonal.pdf}\end{overpic}

Figure 1: Computed distribution density $h$ at time $t=0$ (left), $t=0.0125$
(center), $t=0.5$ (right). The top panels show one view of the computed $h$,
the middle panels show the same view but with a reduced the level of opacity
to see through and highlighted in green the section whose results are reported
in the bottom panels. The white line in all the views represents the average
$q^{\ast}$.

##  3 Our full two-dimensional emotional contagion model

The model we develop is based on the model proposed in Ref. Agnelli2015.

Let $\Omega\subset\mathbb{R}^{2}$ denote a bounded domain. We assume that
pedestrians head to the exit $E$ within the domain. The case of multiple exits
(i.e., $E$ is the finite union of disjoint sets) can be handled too. See,
e.g., Ref. kim_quaini. At the boundary of the domain and within the domain, we
have walls $W$. Let ${\bm{x}}=(x,y)$ denote position and
${{\bm{v}}}=q(\cos\theta,\sin\theta)$ denote velocity, where the velocity
modulus $q$ corresponds to the fear level (like in Section 2), and $\theta$ is
the velocity direction. Unlike Section 2, here pedestrians can change their
walking direction as explained in what follows. For a system composed of a
large number of pedestrians distributed inside $\Omega$, the distribution
function is given by

$f=f(t,{\bm{x}},q,\theta)\quad\text{for all}\,\,\,t\geq
0,\,\,{\bm{x}}\in\Omega,\,\,q\in[0,1],\,\,\theta\in[0,2\pi).$

Under suitable integrability conditions,
$f(t,{\bm{x}},q,\theta)d{\bm{x}}dqd\theta$ represents the number of
individuals who, at time $t$, are located in the infinitesimal rectangle
$[x,x+dx]\times[y,y+dy]$ and have a velocity belonging to
$[q,q+dq]\times[\theta,\theta+d\theta]$.

Following Ref. Agnelli2015, we assume that the crowd size is not enough to
justify the continuity of the distribution function over the variable
$\theta$. Thus, variable $\theta$ is discrete and it can take values in the
set:

$I_{\theta}=\left\\{\theta_{i}=\frac{i-1}{N_{d}}2\pi:i=1,\dots,N_{d}\right\\},$

where $N_{d}$ is the maxim number of possible directions. Then, we can write
the kinetic type representation

$f(t,{\bm{x}},q,\theta)=\sum_{i=1}^{N_{d}}f^{i}(t,{\bm{x}},q)\delta(\theta-\theta_{i}),$
(19)

where $f^{i}(t,{\bm{x}},q)=f(t,{\bm{x}},q,\theta_{i})$ represents the active
particles that, at time $t$ and position ${\bm{x}}$, move with speed $q$ and
direction $\theta_{i}$. In equation (19), $\delta$ denotes the Dirac delta
measure. Notice that a discrete set of directions is also convenient from the
computational point of view. Our model could be extended to handle a continuum
of directions. See, e.g., Ref. Bellomo2019_new.

Let us introduce some reference quantities that will be use to make the
variables non-dimensional. We define:

  * -

$D$: the largest distance a pedestrian can cover in domain $\Omega$;

  * -

$V_{M}$: the highest velocity modulus a pedestrian can reach;

  * -

$T$: a reference time given by $D/V_{M}$;

  * -

$\rho_{M}$: the maximal admissible number of pedestrians per unit area.

The dimensionless variables are then: position $\hat{{\bm{x}}}={\bm{x}}/D$,
time $\hat{t}=t/T$, and distribution function $\hat{f}=f/\rho_{M}$. Notice
that the fear level $q\in[0,1]$ is already a non-dimensional variable.
Therefore, the velocity modulus is a non-dimensional variable and the
dimensional speed is $qV_{M}$. From now on, all the variables will be non-
dimensional and hats will be omitted to simplify notation.

Due to the normalization of $f$, and of each $f^{i}$, the non-dimensional
local density is given by:

$\displaystyle\rho(t,{\bm{x}})=\sum_{i=1}^{N_{d}}\int_{0}^{1}f^{i}(t,{\bm{x}},q)dq.$
(20)

###  3.1 Modeling interactions

Each pedestrian is modeled as a particle. Interactions involve three types of
particles:

  * -

test particles with state $({\bm{x}},\theta_{i})$: they are representative of
the whole system;

  * -

candidate particles with state $({\bm{x}},\theta_{h})$: they can reach in
probability the state of the test particles after individual-based
interactions with the environment or with field particles;

  * -

field particles with state $({\bm{x}},\theta_{k})$: their presence triggers
the interactions of the candidate particles.

The process through which a pedestrian decides the direction to take is the
results of several factors. We take into account four factors:

  1. (F1)

The goal to reach the exit.
Given a candidate particle at the point ${\bm{x}}$, we define its distance to
the exit as

$d_{E}({\bm{x}})=\min_{{\bm{x}}_{E}\in E}||{\bm{x}}-{\bm{x}}_{E}||,$

and we consider the unit vector ${\bm u}_{E}({\bm{x}})$, pointing from
${\bm{x}}$ to the exit. See Figure 2.

  2. (F2)

The need to avoid the collision with the walls.
Given a candidate particle at point ${\bm{x}}$ moving with direction
$\theta_{h}$, we define the distance $d_{W}({\bm{x}},\theta_{h})$ from the
particle to a wall at a point ${\bm{x}}_{W}({\bm{x}},\theta_{h})$ where the
particle is expected to collide with the wall. The unit tangent vector ${\bm
u}_{W}({\bm{x}},\theta_{h})$ to $\partial\Omega$ at ${\bm{x}}_{W}$ points to
the direction of the the exit. Vector ${\bm u}_{W}$ is used to avoid a
collision with the walls. See Figure 2.

  3. (F3)

The tendency to look for less congested area.
A candidate particle $({\bm{x}},\theta_{h})$ may decide to change direction in
order to avoid congested areas. This is achieved with the direction that gives
the minimal directional derivative of the density at the point ${\bm{x}}$. We
denote such direction by unit vector ${\bm u}_{C}(\theta_{h},\rho)$. In the
literature, this behavior is sometimes referred to as _individualistic_. [25]

  4. (F4)

The tendency to follow the stream.
A candidate particle $h$ interacting with a field particle $k$ may decide to
follow it and thus adopt its direction, denoted with unit vector ${\bm
u}_{F}=(\cos\theta_{k},\sin\theta_{k})$. In the literature, this behavior is
also referred to as _herding_. [25]

\begin{overpic}[width=130.08731pt,height=101.17645pt,grid=false,angle={-90},tics=10]{./domain.pdf}
\put(15.0,79.0){\text{$\Omega$}} \put(96.0,47.0){\text{$E$}}
\put(30.0,35.0){\scriptsize{\text{${\bm{x}}$}}}
\put(95.0,3.0){\scriptsize{\text{${\bm{x}}_{W}$}}}
\put(60.0,40.0){\scriptsize{\text{$d_{E}({\bm{x}})$}}}
\put(55.0,10.0){\scriptsize{\text{$d_{W}({\bm{x}})$}}}
\put(43.0,35.0){\scriptsize{\text{${\bm u}_{E}$}}}
\put(100.0,12.0){\scriptsize{\text{${\bm u}_{W}$}}}
\put(38.0,23.0){\scriptsize{\text{$\theta_{h}$}}} \end{overpic} Figure 2:
Sketch of computational domain $\Omega$ with exit $E$ and candidate particle
$h$ located at ${\bm{x}}$, moving with direction $\theta_{h}$. The pedestrian
should choose direction ${\bm u}_{E}$ to reach the exit, while direction ${\bm
u}_{W}$ is to avoid collision with the wall. The distances form the exit and
from the wall are $d_{E}$ and $d_{W}$, respectively.

Factors (F1) and (F2) are related to geometric aspects of the domain, while
factors (F3) and (F4) consider that people’s behavior is influenced by
surrounding crowd. Factors (F3) and (F4) compete with each other: (F4) is
dominant in a stress situation, while (F3) characterizes rational behavior. To
weight between (F3) and (F4), we use the fear level $q\in[0,1]$.

The interaction with the bounding walls is modeled with two terms:

  * -

$\mu[\rho]$: the interaction rate models the frequency of interactions between
candidate particles and the walls and/or an exit. If the local density is low,
it is easier for pedestrians to see the walls and doors. Thus, we set
$\mu[\rho]=1-\rho$.

  * -

$\mathcal{A}_{h}(i)$: the transition probability gives the probability that a
candidate particle $h$, i.e. with direction $\theta_{h}$, adjusts its
direction into that of the test particle $\theta_{i}$ due to the presence of
the walls and/or an exit. The following constraint for $\mathcal{A}_{h}(i)$
has to be satisfied:

$\sum_{i=1}^{N_{d}}\mathcal{A}_{h}(i)=1\quad\text{for
all}\,\,h\in\\{1,\dots,N_{d}\\}.$

We assume that particles change direction, in probability, only to an adjacent
direction in the discrete set $I_{\theta}$ for each possible change step. This
means a candidate particle $h$ may end up into the states $h-1,h+1$ or remain
in the state $h$. In the case $h=1$, we set $\theta_{h-1}=\theta_{N_{d}}$ and,
in the case $h=N_{d}$, we set $\theta_{h+1}=\theta_{1}$. The set of all
transition probabilities
$\mathcal{A}=\\{\mathcal{A}_{h}(i)\\}_{h,i=1,\dots,N_{d}}$ forms the so-called
table of games that models the game played by active particles interacting
with the walls.

To take into account factors (F1) and (F2), we define the vector

$\displaystyle{\bm u}_{G}({\bm{x}},\theta_{h})$
$\displaystyle=\frac{(1-d_{E}({\bm{x}})){\bm
u}_{E}({\bm{x}})+(1-d_{W}({\bm{x}},\theta_{h})){\bm
u}_{W}({\bm{x}},\theta_{h})}{||(1-d_{E}({\bm{x}})){\bm
u}_{E}({\bm{x}})+(1-d_{W}({\bm{x}},\theta_{h})){\bm
u}_{W}({\bm{x}},\theta_{h})||}$
$\displaystyle=(\cos\theta_{G},\sin\theta_{G}).$ (21)

Here $\theta_{G}$ is the preferred geometrical direction, which is the ideal
direction that a pedestrian should take in order to reach the exit and avoid
the walls in an optimal way. Notice that the closer a pedestrian is to an exit
(resp., a wall), the larger the weight for direction ${\bm u}_{E}$ (resp.,
${\bm u}_{W}$).

A candidate particle $h$ will update its direction by choosing the angle
closest to $\theta_{G}$ among the three allowed angles
$\theta_{h-1},\theta_{h}$ and $\theta_{h+1}$. The transition probability is
given by:

$\mathcal{A}_{h}(i)=\beta_{h}\delta_{s,i}+(1-\beta_{h})\delta_{h,i},\quad
i=h-1,h,h+1,$ (22)

where

$s=\operatorname*{arg\,min}_{j\in\\{h-1,h+1\\}}\\{d(\theta_{G},\theta_{j})\\},$

with

$d(\theta_{p},\theta_{q})=\begin{cases}|\theta_{p}-\theta_{q}|&\text{if}\,\,\,|\theta_{p}-\theta_{q}|\leq\pi,\\\
2\pi-|\theta_{p}-\theta_{q}|&\text{if}\,\,\,|\theta_{p}-\theta_{q}|>\pi.\end{cases}$
(23)

In (22), $\delta$ denotes the Kronecker delta function. Coefficient
$\beta_{h}$ is defined by:

$\beta_{h}=\begin{cases}1&\text{if}\,\,\,d(\theta_{h},\theta_{G})\geq\Delta\theta,\\\
\dfrac{d(\theta_{h},\theta_{G})}{\Delta\theta}&\text{if}\,\,\,d(\theta_{h},\theta_{G})<\Delta\theta,\end{cases}$

where $\Delta\theta=2\pi/{N_{d}}$. The role of $\beta_{h}$ is to allow for a
transition to $\theta_{h-1}$ or $\theta_{h+1}$ even in the case that the
geometrical preferred direction $\theta_{G}$ is closer to $\theta_{h}$. Such a
transition is more likely to occur the more distant $\theta_{h}$ and
$\theta_{G}$ are. 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
(in the absence of interactions other than with the environment) with
probability 1.

The interaction with other pedestrians is modeled with two terms:

  * -

$\eta[\rho]$: the interaction rate defines the number of binary encounters per
unit time. If the local density increases, then the interaction rate also
increases. For simplicity, we take $\eta[\rho]=\rho$.

  * -

$\mathcal{B}_{hk}(i)[\rho{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0},q}]$:
the transition probability gives the probability that a candidate particle $h$
modifies its direction $\theta_{h}$ into that of the test particle $i$, i.e.
$\theta_{i}$, due to the research of less congested areas (factor (F3)) and
the interaction with a field particle $k$ that moves with direction
$\theta_{k}$ (factor (F4)). The following constrain for
$\mathcal{B}_{hk}(i)[{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\rho,q]}$
has to be satisfied:

$\sum_{i=1}^{N_{d}}\mathcal{B}_{hk}(i)[\rho,q]=1\quad\text{for
all}\,\,h,k\in\\{1,\dots,N_{d}\\},$

where the square brackets denote the dependence on the density $\rho$ and fear
level $q$.

The game consists in choosing the less congested direction among the three
admissible ones. This direction can be computed for a candidate pedestrian $h$
situated at ${\bm{x}}$, by taking

$C=\operatorname*{arg\,min}_{j\in\\{h-1,h,h+1\\}}\\{\partial_{j}\rho(t,{\bm{x}})\\},$

where $\partial_{j}\rho$ denotes the directional derivative of $\rho$ in the
direction given by angle $\theta_{j}$. We have ${\bm
u}_{C}(\theta_{h},\rho)=(\cos\theta_{C},\sin\theta_{C})$. As for the tendency
to follow the crowd, we set ${\bm u}_{F}=(\cos\theta_{k},\sin\theta_{k})$.
This means that a candidate particle follows the direction of a field
particle.

To take into account factors (F3) and (F4), we define the vector

$\displaystyle{\bm u}_{P}(\theta_{h},\theta_{k},\rho,q)=\frac{q{\bm
u}_{F}+(1-q){\bm u}_{C}(\theta_{h},\rho)}{||q{\bm u}_{F}+(1-q){\bm
u}_{C}(\theta_{h},\rho)||}=(\cos\theta_{P},\sin\theta_{P}),$ (24)

where the subscript $P$ stands for pedestrians. Direction $\theta_{P}$ is the
preferred interaction-based direction, obtained as a weighted combination
between the tendency to follow the stream and the tendency to avoid crowded
areas.

The transition probability is given by:

$\mathcal{B}_{hk}(i)[\rho,q]=\beta_{hk}\rho\delta_{r,i}+(1-\beta_{hk}\rho)\delta_{h,i},\quad
i=h-1,h,h+1,$

where $r$ and $\beta_{hk}$ are defined by:

$r=\operatorname*{arg\,min}_{j\in\\{h-1,h+1\\}}\\{d(\theta_{P},\theta_{j})\\},$
$\beta_{hk}=\begin{cases}1&\text{if}\,\,\,d(\theta_{h},\theta_{P})\geq\Delta\theta\\\
\dfrac{d(\theta_{h},\theta_{P})}{\Delta\theta}&\text{if}\,\,\,d(\theta_{h},\theta_{P})<\Delta\theta.\end{cases}$

We recall that $d(\cdot,\cdot)$ is defined in (23).

######  Remark 3.1.

The kinetic approach to pedestrian dynamics in Ref. Agnelli2015 and
kim_quaini, which will be used for comparison in Section 5.1, models fear in a
simplified way. The first simplification is that the walking speed is not
affected by fear. The walking direction is adapted through a parameter
$\epsilon$ that takes the place of $q$ in (24). Unlike $q$, $\epsilon$ is
constant in space and time. So, the second simplification is that all the
people in the simulation have the same fear level that remains constant for
the entire time interval under consideration.
The kinetic type model in Ref. Bellomo2019_new uses a more realistic
description of stress behavior. In the future, it will be interesting to
compare the results produced by the model in this paper with the results
obtained from the model in Ref. Bellomo2019_new.

###  3.2 The mathematical model

The derivation of the mathematical model can be obtained by a suitable 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.
Taking into account factors (F1)-(F4), we obtain:

$\displaystyle\frac{\partial f^{i}}{\partial t}$
$\displaystyle+\nabla\cdot\left(q(\cos\theta_{i},\sin\theta_{i})^{T}f^{i}(t,{\bm{x}},q)\right)$
$\displaystyle=\mathcal{J}^{i}[f](t,{\bm{x}},q)+\gamma\frac{\partial((q-q^{\ast})f^{i})}{\partial
q}$
$\displaystyle=\eta[\rho]\left(\sum_{h,k=1}^{N_{d}}\int_{0}^{1}\int_{0}^{1}\mathcal{B}_{hk}(i)f^{h}(t,{\bm{x}},q)f^{k}(t,{\bm{x}},q)dqdq-f^{i}(t,{\bm{x}},q)\rho(t,{\bm{x}})\right)$
$\displaystyle\quad+\mu[\rho]\left(\sum_{h=1}^{N_{d}}\int_{0}^{1}\mathcal{A}_{h}(i)f^{h}(t,{\bm{x}},q)dq-f^{i}(t,{\bm{x}},q)\right)$
$\displaystyle\quad+\gamma\frac{\partial((q-q^{\ast})f^{i}(t,{\bm{x}},q))}{\partial
q}$ (25)

for $i=1,2,\dots,N_{d}$. Functional $\mathcal{J}^{i}[f]$ represents the net
balance of particles that move with direction $\theta_{i}$ due to interactions
with the environment and with the surrounding people. Notice that the last
term at the right-hand-side in (25) corresponds to the term at the right-hand-
side in (5), i.e. it accounts for the fact that people tend to equilibrate
their fear level with the average fear level they perceive.

Equation (25) is completed with equation (20) for the density. In the next
section, we will discuss a numerical method for the solution of problem (20),
(25).

##  4 Numerical method

Our approach is based on a splitting method. The idea is to split model (25)
into subproblems that are easier to solve and for which practical algorithms
are readily available. The numerical method is then completed by picking an
appropriate numerical scheme for each subproblem. Among the available
operator-splitting methods, we chose the Lie splitting scheme because it
provides a good compromise between accuracy and robustness, as shown in Ref.
glowinski2003finite.

###  4.1 The Lie operator-splitting scheme

Before applying it to problem (20),(25), we briefly present the Lie splitting
scheme for a generic method that is first-order system in time:

$\displaystyle\frac{\partial\phi}{\partial t}+A(\phi)$ $\displaystyle=$
$\displaystyle 0,\quad\textrm{in}\ (0,T),$ $\displaystyle\phi(0)$
$\displaystyle=$ $\displaystyle\phi_{0},$

where A is an operator from a Hilbert space into itself. Operator A is then
split, in a non-trivial decomposition, as

$A=\sum\limits_{i=1}^{I}A_{i}.$

The Lie scheme consists of the following. Let $\Delta t>0$ be a time
discretization step for the time interval $[0,T]$. Denote $t^{k}=k\Delta t$,
with $k=0,\dots,N_{t}$ and let $\phi^{k}$ be an approximation of
$\phi(t^{k}).$ Set $\phi^{0}=\phi_{0}.$ For $n\geq 0$, compute $\phi^{k+1}$ by
solving

$\displaystyle\frac{\partial\phi_{i}}{\partial t}+A_{i}(\phi_{i})$
$\displaystyle=$ $\displaystyle 0\quad\textrm{in}\;(t^{k},t^{k+1}),$ (26)
$\displaystyle\phi_{i}(t^{k})$ $\displaystyle=$
$\displaystyle\phi^{k+(i-1)/I},$ (27)

and then set $\phi^{k+i/I}=\phi_{i}(t^{k+1}),$ for $i=1,\dots.I.$

This method is first-order accurate in time. More precisely, if (4.1) is
defined on a finite-dimensional space, and if the operators $A_{i}$ are smooth
enough, then $\|\phi(t^{k})-\phi^{k}\|=O(\Delta t)$. [15]

We will apply Lie splitting to problem (20),(25), which will get split into
two subproblems:

  1. 1.

An advection problem that includes the term for emotional contagion at the
right-hand-side in eq. (25).

  2. 2.

A problem involving the interaction with the environment and other
pedestrians.

###  4.2 Lie scheme applied to problem (25)

Let us apply the Lie operator-splitting scheme to problem (20),(25). Given an
initial condition $f^{i,0}=f^{i}(0,{\bm{x}},q)$, for $i=1,\dots,N_{d}$, the
algorithm reads: For $n=0,1,2,\dots,N_{t}-1,$ perform the following steps:

  * -

Step 1: Find $f^{i}$, for $i=1,\dots,N_{d}$, such that

$\begin{cases}\dfrac{\partial f^{i}}{\partial
t}+\nabla\cdot\left(q(\cos\theta_{i},\sin\theta_{i})^{T}f^{i}(t,{\bm{x}},q)\right)=\gamma\frac{\partial((q-q^{\ast})f^{i}(t,{\bm{x}},q))}{\partial
q}\,\,\,\text{on }(t^{n},t^{n+1}),\\\
f^{i}(t^{n},{\bm{x}},q)=f^{i,n}.\end{cases}$ (28)

Set $f^{i,n+\frac{1}{2}}=f^{i}(t^{n+1},{\bm{x}})$.

  * -

Step 2: Find $f_{i}$, for $i=1,\dots,N_{d}$, such that

$\begin{cases}\dfrac{\partial f^{i}}{\partial
t}=\mathcal{J}^{i}[f](t,{\bm{x}},q)\,\,\,\text{on }(t^{n},t^{n+1}),\\\
f^{i}(t^{n},{\bm{x}},q)=f^{i,n+\frac{1}{2}}.\end{cases}$ (29)

Set $f^{i,n+1}=f^{i}(t^{n+1},{\bm{x}})$.

Notice that once $f^{i,n+1}$ is computed for $i=1,\dots,N_{d}$, we use
equation (20) to get the density $\rho^{n+1}$.

We associate with $M$ a time discretization step $\tau=(t_{f}-t_{0})/M$ and
set $t^{m}=t_{0}+m\tau$. For the discretization of problem (28), we use the
scheme reported in Section 2.1. As for the approximation of the problem (29),
we use the following explicit scheme:

$f_{j,k,l}^{i,m+1}=f_{j,k,l}^{i,m}+\tau\Big{(}\mathcal{J}^{m}[F^{m}]\Big{)},$

where $F^{m}$ is the approximation of distribution function (19) at time
$t^{m}$.

##  5 Numerical results

This section presents some sample simulations focused on the computation of
evacuation time from a room with one exit. We select three case studies in
order to (i) compare with the results obtained from a model that features no
emotional contagion, (ii) analyze the role of a geometric parameter (exit
size), model parameters (interaction radius and contagion interaction
strength), and initial fear distribution, and (iii) show that our model can
roughly reproduce experimental data.

For all the simulations, we consider eight different velocity directions
$N_{d}=8$ in the discrete set:

$I_{\theta}=\left\\{\theta_{i}=\frac{i-1}{8}2\pi:i=1,\dots,8\right\\}.$

###  5.1 Test 1: evacuation from a room with one door of variable width

This first test case is inspired from Ref. Agnelli2015. The computational
domain contains a square room with side $10$ m with one exit door located in
the middle of the right side. The exit size is $2.6$ m. The computational
domain is larger than the room itself to follow the motion of the pedestrian
also once they have left the room. We aim at simulating the evacuation of 46
people located inside the room and initially distributed into two equal-area
circular clusters. See Figure 3, first panel on the second row. The two groups
are initially moving against the each other with opposite initial directions
$\theta_{3}$ and $\theta_{7}$. For the purpose of comparing with the crowd
dynamics given by the model in Ref. kim_quaini, which features no emotional
contagion (see Remark 3.1), in this test we assign the initial fear level as
follows: the group with initial direction $\theta_{7}$ has an initial high
fear level ($q=0.8$), while the other group has an initial low fear level
($q=0.2$). We set contagion interaction strength $\gamma=1$ and interaction
distance $R=0.5$ m. This value of $R$ means that people adjust their fear
level to the average fear level of the people in their immediate vicinity.

In order to work with dimensionless quantities as described in Section 3, we
define the following reference quantities: $D=10\sqrt{2}$ m, $V_{M}=2$ m/s,
$T=D/V_{M}=5\sqrt{2}$ s, and $\rho_{M}=7$ people/m2. However, once the results
are computed we convert them back to dimensional quantities.

We consider two different meshes:

  * -

_coarse mesh_ with $\Delta x=\Delta y=0.5$ m and $\Delta q=0.05$;

  * -

_fine mesh_ with $\Delta x=\Delta y=0.25$ m and $\Delta q=0.05$.

Notice that the level of refinement changes only in the spatial directions. We
choose the time step to $\Delta t=0.0375$ s for both meshes. This time step is
such that $\tau=\Delta t/M$ satisfies the CLF condition (17) for $M=3$. Figure
3 shows the density computed with the fine mesh at times $t=0$, 1.5, 3, 6,
10.5, 13.5 s. At time $t=6$ s, we observe the formation of a region of high
density next to a region of low density, with a sharp divide between the two
that is aligned with one door end. There are several possible reasons for
this. First, the walking direction results from the interplay of different
factors and thus it may not be intuitive. We consider only 8 discrete
directions as stated at the beginning of Sec. 5. Maybe a higher number of
directions would lead to a less sharp divide, although it would not be
realistic for the small group of people under consideration (see Sec. 3).
Finally, the value of $R$ is small for this test, which means that it takes
longer for fear to propagate. As the group with high fear reaches the group
with low fear, it takes a while for the emotion to spread through the entire
low fear group.

\begin{overpic}[width=138.76157pt,grid=false]{./ML_combine_contagion_time_0.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./ML_combine_contagion_time_40.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./ML_combine_contagion_time_80.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./ML_combine_Density_time_0.pdf}
\put(24.0,54.2){\vector(0,-1){10.0}} \put(24.0,22.0){\vector(0,1){10.0}}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./ML_combine_Density_time_40.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./ML_combine_Density_time_80.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./ML_combine_contagion_time_160.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./ML_combine_contagion_time_280.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./ML_combine_contagion_time_360.pdf}
\end{overpic}

\begin{overpic}[width=138.76157pt,grid=false]{./ML_combine_Density_time_160.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./ML_combine_Density_time_280.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./ML_combine_Density_time_360.pdf}
\end{overpic}

Figure 3: Test 1: Evacuation process of 46 pedestrians grouped into two
clusters with opposite initial directions $\theta_{3}$ and $\theta_{7}$ using
the fine mesh: average fear level $q^{\ast}$ (top panel for each time) and
people density (bottom panel for each time) for $t=0,1.5,3,6,10.5,13.5$ s. The
density plots show also the computational mesh.

Figure 4 reports the number of pedestrians left in the room computed with the
coarse and fine mesh: we can see very good agreement for the two curves. In
addition, in Figure 4 we report the results obtained with no emotional
contagion, i.e. with the model in Ref. kim_quaini with either low fear
($\epsilon=0.2$) or high fear ($\epsilon=0.8$) in the entire domain (see
Remark 3.1) and for both meshes under consideration. Since the model with no
emotional contagion does not have a stringent condition like (17) on the time
step, we set $\Delta t=1.5$ s for the coarse mesh and $\Delta t=0.75$ s for
the fine mesh. See Ref. kim_quaini for more details. Fig. 4 shows that the
model with no emotional contagion produces very similar results for
$\epsilon=0.2$ (low fear) and $\epsilon=0.8$ (high fear), indicating that its
description of fear-induced behavior is too simplistic. The model with
emotional contagion predicts a faster evacuation dyanmics. Recall that in the
model with emotional contagion fear affects both the walking speed (which
increases with the fear level) and the walking direction, while in the model
with no emotional contagion it influences only the latter. Our results in
Figures 3 and 4 confirm what observed in Ref. Bellomo2019_new: the spreading
of a stress condition significantly affects the overall evacuation dynamics
and leads to overcrowding in certain areas of the domain.

\begin{overpic}[width=260.17464pt,grid=false]{./Total_Evacuation_Combine_1111.pdf}
\end{overpic} Figure 4: Test 1: number of pedestrians inside the room over
time computed with the model with and without emotional contagion, for both
coarse and fine mesh. In the model with no emotional contagion, we set either
low fear ($\epsilon=0.2$) or high fear ($\epsilon=0.8$) in the entire domain.

Next, we consider the same room but we vary the exit size. We locate the exit
symmetrically with respect to the room centerline and let the exit size vary
from 1.5 m to 4 m. We consider the coarse mesh mentioned before since Figure 4
show that is an appropriate choice for the problem under consideration. All
the other model and discretization parameters are unchanged. Figure 5 shows
the total evacuation time as a function of the exit size for the model with
and without emotional contagion. We see that for all exit sizes the model with
emotional contagion predicts a shorter evacuation time. The total evacuation
time obtained with both models decreases with the exit size, but that once the
exit is large enough for the crowd contained in the room the evacuation time
does not change significantly if the exit is further enlarged. This plateau is
reached at a smaller exit size for the model with emotional contagion.
Finally, we notice that our results for the model with no emotional contagion
are in very good agreement with the results reported in Ref. Agnelli2015.

\begin{overpic}[width=260.17464pt,grid=false]{./Line_Exit_compare_1103.pdf}
\end{overpic} Figure 5: Test 1: total evacuation time as a function of the
exit size for the model with and without emotional contagion.

The existing literature is not conclusive on how “herding” (i.e., the tendency
to follow the stream) influences the efficiency of collective crowd egress.
[25] Some works suggest that herding behavior is detriment to efficient
evacuation (e.g., Refs. ZHENG20114627,HAGHANI201949), some others conclude
that there may be scenarios where herding is beneficial to an escaping crowd
(see, e.g., Ref. KIRCHNER2002260), while according to other works the effect
of herding is still unestablished (see, e.g., Ref.
Pan_2007,Moussaid_2016,LOVREGLIO2016421). The phenomenon of mixed strategy
(i.e., a mixture of factors F3 and F4 in Sec. 3.1) has been investigated by
several numerical studies. Also on the efficacy of the mixed strategy, the
literature is inconclusive. Some works (e.g., Ref. Helbing2000) suggest that a
crowd can benefit from mixed strategies, while other (e.g., Ref.
HAGHANI2019101064) suggest that any percentage of herding strategy within the
crowd has a negative impact on the evacuation efficiency. Our results agree
with the former.

###  5.2 Test 2: evacuation from a room with one door and variable interaction
radius

The computational domain is the same as for the test in Sec. 5.1, with exit
size $2.6$ m. The number of pedestrian, their initial distribution and walking
direction are as in Sec. 5.1. The difference with respect to test in Sec. 5.1
is in the initial fear distribution. The group centered at $(-0.25,-0.25)$
(see Fig. 3, left-most panel on the second row), which occupies a circle of
radius 0.1414 and has initial direction $\theta=3$, has the following initial
probability function:

$f^{3}(0,{\bm{x}},q)=g_{1}(x,y)\,\delta_{0}(q),$

where

$g_{1}(x,y)=-25.8(x+0.25)^{2}-25.8(y+0.25)^{2}+0.52,$

and $\delta_{n}(q)$ is the Kronecker delta defined by $\delta_{n}(q)=1$ if
$q=n$ and $\delta_{n}(q)=0$ if $q\neq n$. As for the group centered at
$(-0.25,0.25)$, which also occupies a circle of radius 0.1414 but has initial
direction $\theta=7$, we have the following initial probability function:

$f^{7}(0,{\bm{x}},q)=g_{2}(x,y)\,\Big{[}(1/2)\,\delta_{1/4}(q)+(1/2)\,\delta_{3/4}(q)\Big{]},$

where

$g_{2}(x,y)=-25.8(x+0.25)^{2}-25.8(y-0.25)^{2}+0.52.$

See the $x=-0.25$ section of the initial distribution function
$f(0,{\bm{x}},q)$ and initial average fear level in Fig. 6.

\begin{overpic}[width=173.44534pt,grid=false]{./fine_slice_time_0.pdf}
\end{overpic}\begin{overpic}[width=173.44534pt,grid=false]{./fine_contagion_time_0.pdf}
\end{overpic}

Figure 6: Test 2: (left) $x=-0.25$ section of the initial distribution
function $f(0,{\bm{x}},q)$ and (right) initial average fear level.

We simulate evacuation from the room using three different choices of $R$:
$R=0.5$ (like in Sec. 5.1), $R=1$, and $R=10$, which corresponds to the case
where the fear level of each person is influenced by the fear level of
everybody in the room. Figure 7 reports the number of pedestrians inside the
room over time computed with the fine mesh for the three values of $R$. Figure
7 shows that there is a big difference in the evacuation dynamics when going
from $R=0.5$ to $R=1$, while the change is not substantial when going from
$R=1$ to $R=10$. See also Fig. 8, which shows the evolution of the average
fear level. When $R=0.5$, part of the group with an initially higher fear
level that does not get in contact with the lower fear group and leaves the
room quickly. This explains the initially fast evacuation dynamics for $R=0.5$
seen in Figure 7. A part of the group with an initially lower fear level gets
in contact with the higher fear group and thus picks up speed. Nonetheless,
the emotion spreads slowly and in fact the last 23 people in the room take
longer to leave the room. For $R=1,10$, the emotion spreads fast and people
evacuate much faster.

\begin{overpic}[width=260.17464pt,grid=false]{./Total_Evacuation_times_Rs_1122.pdf}
\end{overpic} Figure 7: Test 2: number of pedestrians inside the room over
time computed for $R=0.5,1,10$.

\begin{overpic}[width=108.405pt,grid=false]{./Test2_a_Rs_R_0_5_contagion_time_40_cut.pdf}
\put(-40.0,50.0){$R=0.5$}
\end{overpic}\begin{overpic}[width=108.405pt,grid=false]{./Test2_a_Rs_R_0_5_contagion_time_80_cut.pdf}
\end{overpic}\begin{overpic}[width=108.405pt,grid=false]{./Test2_a_Rs_R_0_5_contagion_time_160_cut.pdf}
\end{overpic}\begin{overpic}[width=108.405pt,grid=false]{./Test2_a_Rs_R_1_contagion_time_40_cut.pdf}
\put(-40.0,50.0){$R=1$}
\end{overpic}\begin{overpic}[width=108.405pt,grid=false]{./Test2_a_Rs_R_1_contagion_time_80_cut.pdf}
\end{overpic}\begin{overpic}[width=108.405pt,grid=false]{./Test2_a_Rs_R_1_contagion_time_160_cut.pdf}
\end{overpic}\begin{overpic}[width=108.405pt,grid=false]{./Test2_a_Rs_R_10_contagion_time_40_cut.pdf}
\put(-40.0,50.0){$R=10$}
\end{overpic}\begin{overpic}[width=108.405pt,grid=false]{./Test2_a_Rs_R_10_contagion_time_80_cut.pdf}
\end{overpic}\begin{overpic}[width=108.405pt,grid=false]{./Test2_a_Rs_R_10_contagion_time_160_cut.pdf}
\end{overpic}

Figure 8: Test 2: Evolution of the average fear level computed with the fine
mesh for $R=0.5$ (top row), $R=1$ (middle row), and $R=10$ (bottom row).

###  5.3 Test 3: experiment in ant colony

We aim at reproducing one of the experiments in Ref. SHIWAKOTI20111433. In the
experimental setting, 200 ants are placed in a square chamber (31 mm $\times$
31 mm) with exit at the corner of the walls that is 2.5 mm wide. At a certain
time, citronella (an insect repellent) is injected through a small hole
created at the center of the square chamber to create panic. We start our
simulations right after the citronella has been injected. In Ref.
SHIWAKOTI20111433, it is observed that when citronella is introduced, ants
rush to the exit in a manner reminiscent of humans in a crowd panic. The
evacuation time to the nearest second and the number of ants escaped reported
in Ref. SHIWAKOTI20111433 were measured by manual counting from playback of
digital video recordings. The measured mean (out of 10 trials) escape time is
11.2 s ($\pm$ 2.6 s standard deviation) for the first 50 ants. We aim at
showing that with a suitable parameter set-up and initial condition, our model
can match this evacuation time. For this study, the parameters were set using
our experience with the previous tests and by trial-and-error. More
sophisticated parameter estimation techniques will be adopted for a follow-up
paper.

In order to work with dimensionless quantities, we define the following
reference quantities: $D=31\sqrt{2}$ mm, $V_{M}=1$ mm/s,
$T=D/V_{M}=31\sqrt{2}$ s, and $\rho_{M}=4$ ants/mm2.

For the numerical simulation, we choose $\Delta x=\Delta y=1$ mm and $\Delta
q=0.05$. The time step is set to $\Delta t=0,0375$ s in order to satisfy
condition (17). We initially place the ants in a square with side 22 mm with
uniform density. See Fig. 9, left-most panel on the second row. We consider
two initial conditions. Both of them have high fear level ($q=1$) where the
citronella is injected, but they differ in the rest of the domain. In
condition 1 we prescribe low panic ($q=0.1$) everywhere else, while in
condition 2 we assign $q=0.65$ everywhere else. Fig. 9, top left panel, shows
the initial average fear level for condition 2. From Ref. SHIWAKOTI20111433,
we could not guess whether the repellent has some immediate effect on the
entire domain and that is why we opted for two initial conditions. Recall that
the fear level corresponds to the walking speed. We set $\theta_{2}$ as
initial walking direction. Based on the results from test 2, we set $R=1$ and
try different values of the contagion interaction strength $\gamma$.

Fig. 9 shows the evacuation process of the first 50 ants for initial condition
2 and $\gamma=0.1$. We report the computed average fear level and the density
at different times. Around $t=9$ s, we observe the formation of a high density
region close to the narrow exit, which gets worse as time passes.

\begin{overpic}[width=138.76157pt,grid=false]{./contagion_time_0.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./contagion_time_5.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./contagion_time_9.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./Density_time_0.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./Density_time_5.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./Density_time_9.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./contagion_time_17.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./contagion_time_25.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./contagion_time_29.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./Density_time_17.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./Density_time_25.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt,grid=false]{./Density_time_29.pdf}
\end{overpic}

Figure 9: Test 3: Evacuation process of 200 ants with initial direction
$\theta_{2}$ and $\gamma=0.1$: average fear level $q^{\ast}$ (top panel for
each time) and people density (bottom panel for each time) for
$t=0,1.5,3,6,9,10.875$ s. The density plots show also the computational mesh.

Figure 10 reports the number of ants inside the room over time computed with
initial condition 2 and $\gamma=0.1$ (i.e., same case as in Fig. 9), and
initial condition 1 and $\gamma=0.1,0.5,1$. With condition 2 and $\gamma=0.1$,
we obtain that the first 50 ants leave the room in 10.7 s. This compares
favorably with the experimental data in Ref. SHIWAKOTI20111433 (11.2 s $\pm$
2.6 s). With initial condition 1 instead, for no value of $\gamma$ we can
match the experimental data. In addition, notice the large difference in
evacuation dynamics for the different values of $\gamma$.

\begin{overpic}[width=260.17464pt,grid=false]{./Total_Evacuation_Time.pdf}
\end{overpic} Figure 10: Test 3: number of ants inside the room over time
computed with two initial conditions and different values of $\gamma$. The
graph also reports the average time for the first 50 ants to leave the room in
the experiments.

##  6 Conclusion

In this paper, we have utilized a kinetic-theoretic approach to model crowd
dynamics in the presence of a propagating emotional contagion such as fear. We
have applied a Lie splitting scheme to simulate our model. This scheme
separates the overall task into an advection problem that accounts for the
emotional contagion and a second problem that accounts for the physical
interactions of pedestrians with the environment and one another.

We have completed three tests involving evacuation from a room. One of our
tests highlights how crowd dynamics depend on whether an emotional contagion
is present or absent. Indeed, it is known that social phenomena can modify
interaction rules. Our numerical results confirm that the propagation of an
emotional contagion can significantly affect evacuation dynamics. We have
shown that evacuation dynamics depend sensitively on two key model parameters:
the interaction radius and the contagion interaction strength. Finally, we
have shown that our model can reproduce experimental data generated by
distressed ants provided that we appropriately tune initial data and model
parameters.

In the future, we plan to thoroughly validate our model using data generated
by studies of the collective behavior of panicked animals. Naturally there
exist fewer experimental studies of collective human dynamics under emergency
conditions.

## Acknowledgements

This work has been partially supported by the National Science Foundation
through grant DMS-1816315 (WO) and DMS-1620384 (AQ).

## References

  * [1] J. P. Agnelli, F. Colasuonno, and D. Knopoff,  A kinetic theory approach to the dynamics of crowd evacuation from bounded domains,  Math. Models Methods Appl. Sci. 25 (2015) 109–129.
  * [2] V. V. Aristov,  Biological systems as nonequilibrium structures described by kinetic methods,  Results Phys. 13 (2019) 102232.
  * [3] B. Aylaj, N. Bellomo, L. Gibelli, and A. Reali,  A unified multiscale vision of behavioral crowds,  Math. Models Methods Appl. Sci. 30 (2020) 1–22.
  * [4] N. Bellomo and A. Bellouquid,  On multiscale models of pedestrian crowds from mesoscopic to macroscopic,  Comm. Math. Sciences. 13 (2015) 1649–1664.
  * [5] N. Bellomo, A. Bellouquid, L. Gibelli, and N. Outada,  A Quest Towards a Mathematical Theory of Living Systems, In Modeling and Simulation in Science, Engineering and Technology (Birkhauser, 2017).
  * [6] N. Bellomo, A. Bellouquid, and D. Knopoff,  From the microscale to collective crowd dynamics,  SIAM Multiscale Model. Simul. 11 (2013) 943–963.
  * [7] N. Bellomo and L. Gibelli,  Toward a mathematical theory of behavioral-social dynamics for pedestrian crowds,  Math. Models Methods Appl. Sci. 25 (2015) 2417–2437.
  * [8] N. Bellomo, L. Gibelli, and N. Outada,  On the interplay between behavioral dynamics and social interactions in human crowds,  Kinet. Relat. Models 12 (2019) 397–409.
  * [9] A. L. Bertozzi, J. Rosado, M. B. Short, and L. Wang,  Contagion shocks in one dimension,  J. Stat. Phys. 158 (2015) 647–664.
  * [10] T. Bosse, R. Duell, Z. A. Memon, J. Treur, and C. Natalie van der Wal,  A multi-agent model for mutual absorption of emotions,  In J. Otamendi, A. Bargiela, J. L. Montes, and L. M. D. Pedrera, editors, Proc. 23rd European Conf. Modelling and Simulation (ECMS’09), (European Council on Modeling and Simulation, 2009), pp. 212–218.
  * [11] C. Cercignani, R. Illner, and M. Pulvirenti,  The Mathematical Theory of Dilute Gases,  (Springer-New York, 1994).
  * [12] E. Cristiani, B. Piccoli, and A. Tosin,  Multiscale Modeling of Pedestrian Dynamics,  (Springers-Cham, 2014).
  * [13] P. Degond, C. Appert-Rolland, M. Moussaïd, J. Pettré, and G. Theraulaz,  A hierarchy of heuristic-based models of crowd dynamics,  J. Stat. Phys. 152 (2013) 1033–1068.
  * [14] P. Degond, J. Liu, S. Merino-Aceituno, and T. Tardiveau,  Continuum dynamics of the intention field under weakly cohesive social interaction,  Math. Models Methods Appl. Sci. 27 (2017) 159–182.
  * [15] R. Glowinski,  Finite element methods for incompressible viscous flow,  in Handbook of Numerical Analysis,  P. G. Ciarlet and J. L. Lions eds., Vol. 9  (North-Holland, 2003).
  * [16] M. Haghani and M. Sarvi,  Heterogeneity of decision strategy in collective escape of human crowds: On identifying the optimum composition,  Int. J Disast. Risk Re. 35 (2019) 101064.
  * [17] M. Haghani and M. Sarvi,  Imitative (herd) behaviour in direction decision-making hinders efficiency of crowd evacuation processes,  Saf. Sci. 114 (2019) 49–60.
  * [18] D. Helbing,  Traffic and related self-driven many-particle systems,  Rev. Mod. Phys. 73 (2001) 1067–1141.
  * [19] D. Helbing, I. Farkas, and T. Vicsek,  Simulating dynamical features of escape panic,  Natures 407 (2000) 487–490.
  * [20] D. Helbing and A. Johansson,  Pedestrian, Crowd and Evacuation Dynamics,  (Springer-New York, 2009), pp. 6476–6495.
  * [21] D. Helbing, A. Johansson, and H. Z. Al-Abideen,  Dynamics of crowd disasters: An empirical study,  Phys. Rev. E 75 (2007) 046109.
  * [22] R. L. Hughes,  The flow of human crowds,  Annu. Rev 35 (2003) 169–182.
  * [23] D. Kim and A. Quaini,  A kinetic theory approach to model pedestrian dynamics in bounded domains with obstacles,  Kinet. Relat. Models 12 (2019) 1273–1296.
  * [24] A. Kirchner and A. Schadschneider,  Simulation of evacuation processes using a bionics-inspired cellular automaton model for pedestrian dynamics,  Physica A 312 (2002) 260–276.
  * [25] D. F. Llorca, M. Haghani, E. Cristiani, N. W. F. Bode, M. Boltes, and A. Corbetta,  Panic, irrationality, and herding: Three ambiguous terms in crowd dynamics research,  J Adv. Transport 2019 (2019) 9267643.
  * [26] R. Lovreglio, A. Fonzone, L. dell’Olio, and D. Borri,  A study of herding behaviour in exit choice during emergencies based on random utility theory,  Saf. Sci. 82 (2016) 421–431.
  * [27] M. Moussaïd, D. Helbing, S. Garnier, A. Johansson, M. Combe, and G. Theraulaz,  Experimental study of the behavioural mechanisms underlying self-organization in human crowds,  Proc. Royal Soc. Lond. B Biol. Sci. 276 (2009) 2755–2762.
  * [28] M. Moussaïd, M. Kapadia, T. Thrash, R.W. Sumner, M. Gross, D. Helbing, and C. Hölscher,  Crowd behaviour during high-stress evacuations in an immersive virtual environment,  J. R. Soc. Interface 13 (2016) 1320160414.
  * [29] X. Pan, C. S. Han, K. Dauber, and K. H. Law,  A multi-agent based framework for the simulation of human and social behaviors during emergency evacuations,  AI Soc. 22 (2007) 113–132.
  * [30] E. Ronchi, F. Nieto Uriz, X. Criel, and P. Reilly,  Modelling large-scale evacuation of music festivals,  Case Studies in Fire Safety 5 (2016) 11–19.
  * [31] N. Shiwakoti, M. Sarvi, G. Rose, and M. Burd,  Animal dynamics based approach for modeling pedestrian crowd egress under panic conditions,  Transport Res. B-Meth. 45 (2011) 1433–1449.
  * [32] J. Tsai, E. Bowring, S. Marsella, and M. Tambe,  Empirical evaluation of computational emotional contagion models,  in Intelligent Virtual Agents, H. H. Vilhjálmsson, S. Kopp, S. Marsella, and K. R. Thórisson, eds., (Springer, 2011), pp. 384–397.
  * [33] A.U. K. Wagoum, A. Tordeux, and W. Liao,  Understanding human queuing behaviour at exits: An empirical study,  R. Soc. open sci. 4 (2017) 160896.
  * [34] L. Wang, M. B. Short, and A. L. Bertozzi,  Efficient numerical methods for multiscale crowd dynamics with emotional contagion,  Math. Models Methods Appl. Sci. 27 (2017) 205–230.
  * [35] N. Wijermans, C. Conrado, M. Van Steen, C. Martella, and J. Li,  A landscape of crowd-management support: An integrative approach,  Saf. Sci. 86 (2016) 142–164.
  * [36] N. Wijermans, R. Jorna, W. Jager, T. van Vliet, and O. M. J. Adang,  CROSS: Modelling crowd behaviour with social-cognitive agents,  J. Artif. Soc. 16 (2013) 1–14.
  * [37] X. Zheng and Y. Cheng,  Modeling cooperative and competitive behaviors in emergency evacuation: A game-theoretical approach,  Comput. Math. with Appl. 62 (2011) 4627–4634.