Abstract
We are interested in modeling and simulating the dynamics of human crowds, where the spreading of an emotion (specifically fear) influences the pedestrians’ behavior. Our focus is on crowd dynamics in venues where dense aggregations might occur within a rarefied crowd (e.g., an airport terminal) and emotional states evolve in space and time as the result of a threat (e.g., a gunshot). In the parts of the venue where crowd density is low, we consider a microscopic, individual-based model inspired by Newtonian mechanics. In this model, the fear level of each pedestrian influences their walking speed and is affected by the fear levels of the people in their vicinity. The mesoscopic model is derived from the microscopic model via a mean-field limit approach. This ensures that the two types of models are based on the same principles and analogous parameters. The mesoscopic model is adopted in the parts of the venue where crowd density is higher, i.e., we use the crowd density as a regime indicator. We propose interface conditions to be imposed at the boundary between the regions of the domain where microscopic and mesoscopic models are used. We note that we do not consider dangerously high-density crowd scenarios, for which a macroscopic (continuum) model would be more appropriate. We test our microscopic-to-mesoscopic model on problems involving a crowd walking through a corridor or evacuating from a square.
1 Introduction
Crowds have the ability to act, think, and adapt their own movement strategies in response to environmental cues. This feature, possessed by living systems in general, sets them apart from inert matter, which follows deterministic laws of physics and behaves predictably under given conditions. It also makes the mathematical modeling of crowd dynamics particularly challenging. In fact, to be accurate even in situations when people alter their typical walking strategy in response to an external cue (e.g., a gunshot or a propagating fire), a model of crowd dynamics needs to go beyond the principles of classical mechanics. Specifically, it needs to incorporate mechanisms to account for heterogeneous behavior of individual entities and their interactions, together with the effect of an emotional state (e.g., fear) on both individuals and interactions. See, e.g., [1–7]. Such a model could be instrumental in improving crowd management in venues like airport terminals or concert arenas.
Microscopic, such as agent-based models, effectively capture individual variability since they allow each agent to adopt a unique walking strategy. See, e.g., [8–12]. For a comprehensive overview of agent-based models that implement emotional contagion in crowd simulations, we refer the reader to [13]. Unfortunately, microscopic models struggle with scalability. Specifically, they cannot accurately represent complex, non-local, and nonlinearly additive interactions as crowd size increases. Mesoscopic or kinetic models (see, e.g. [5–7, 14–18]) offer an alternative that can reproduce complex interactions in crowds. Inspired by the kinetic theory of gases, these models consider individuals as active particles, instead of passive particles, i.e., particles that interact in a conservative and reversible manner as in traditional gas dynamics. The consequence is that interactions in kinetic models for crowd dynamics are irreversible, non-conservative and, in some cases, nonlocal and nonlinearly additive [6, 7, 15]. Therefore, kinetic models offer a powerful tool to overcome the limitations of agent-based models in capturing complex interactions. However, these models, which are Boltzmann-type evolution equations for the statistical distribution of pedestrian positions and velocities, become computationally expensive when the size of the computational domain is large.
This paper proposes a hybrid approach to the study of crowds under the influence of a spreading emotion aiming to combine the advantages of microscopic and kinetic models, while lessening the respective limitations. The starting point for our hybrid model is a microscopic model from [19] in which the emotional state (i.e., the fear level) of each person influences their walking speed and is affected by the emotional state of the people in their vicinity. Note that the emotional state is introduced as a variable that, in response to interactions with other people, can change in space and time and, in turn, alter the walking strategy, as it happens in real-life situations [4]. Through a mean-field limit approach, from the microscopic model we derived a (Bhatnagar–Gross–Krook-like) kinetic model. This operation in 1D is presented in [19], while here we extend it to 2D and propose how the kinetic and microscopic models can be interfaced, so that each model can be used where its use is most appropriate. Specifically, the microscopic model is adopted in the parts of the domain where the crowd density is lower than a prescribed value since in this case such model is computationally cheap and accurate. Where the density is above the prescribed value and the microscopic model loses accuracy and computational efficiency, we adopt the kinetic model. Since we limit the size of the domain where the kinetic model holds, our hybrid microscopic-to-mesoscopic approach remains computationally efficient. Figure 1 shows a picture of a real crowd with mixed density: the vertical black line separates the higher density crowd (left), which would be simulated with the kinetic model, from the lower density crowd (right), which would be simulated with the microscopic model.
FIGURE 1
In addition to the advantages discussed so far, our hybrid model stems from a multiscale vision, i.e., the dynamics at the microscopic scale define the conceptual basis toward the derivation of the models at a higher scale (mesoscopic). The reader interested in the further derivation of a macroscopic model from the kinetic model is referred to [19], where results are limited to one space dimension. At the macroscopic level, the crowd is seen as a continuum and thus macroscopic models are suited for densely packed crowds. It is shown in [19] that the kinetic description provides better resolution than the macroscopic model whose viscosity solution becomes incorrect when the characteristics at the particle level cross. Because of this, we do not consider the macroscopic model and how it interfaces to the models at the lower scales.
Following [19], for simplicity we assume that the walking direction is prescribed. This means that our model cannot reproduce the tendency of an individual in a stressful situation to follow the stream of people, also referred to as herding in the literature [20]. The reader interested in kinetic models that allow for a change in walking direction as a results of interactions with other people and the environment (i.e., walls and other kinds of obstacle) is referred to [16–18, 21–25]. We remark that only in [23] the emotional state is treated as a variable, while it is parameterized as a constant in space and time in [16, 18, 21, 22, 25] and as a function of space and time in [17, 24].
We assess our hybrid method through test cases in one and two dimensions. The tests in 1D represent one-directional motion through a corridor, where half of the people have a high fear level and the remaining half have no fear at all. This situation leads to a localized high crowd density as the panicked people walk faster than the non-panicked. The 2D tests also create a localized high-density region, resulting from a scenario similar to a group of panicked people trying to evacuate a square.
The rest of the paper is organized as follows. Section 2 briefly recalls the derivation of the 1D mesoscopic model from the 1D microscopic model from [19], extends the derivation to 2D and then presents the discretization of the 2D microscopic and mesoscopic models. The interface conditions that allow to use microscopic and mesoscopic models in different parts of a given domain are discussed in Section 3. Section 4 presents the numerical results and conclusions are drawn in Section 5.
2 From microscopic to mesoscopic
2.1 In one dimension
We consider a system of agents that are described by their positions in spatial domain and fear levels for time belonging to a time interval of interest . Each agent moves with speed equal to the agent’s fear level and the fear level of an agent tends to the local average of fear levels of nearby agents. The microscopic model, also called agent-based, is described by the following system of ordinary differential equations:for , where , the average fear level for agent , defined as: being an interaction kernel. One of the commonly used formulas for this kernel issee for example, [19]. Parameter in Equation 3 characterizes the radius of interactions. Parameter in Equation 1 describes the interaction strength. For simplicity, we assume that the interaction strengths of all agents are equal, although for one test in Section 4.1 we will make its value variable from one agent to the other.
We note that Equation 1 is simply Newton’s law, assuming that fear is the only driver for motion. In this simple model, every agent has the same mass, which then gets normalized to 1.
The kinetic formulation is derived from Equations 1, 2 via a mean-field limit approach. For this, we denote the empirical distribution density as:where is as in Equation 1, and represents the Dirac delta function. To obtain the kinetic model, let be a test function, multiply Equation 4 by , integrate with respect to space and fear level over the respective domain and differentiate in time:
By manipulating the above equation as explained in [19], one gets the following integral equation:whereand
Since is an arbitrary test function, Equation 5 leads to the following equation:which can be rewritten as:where the local average fear level is given by:
Under appropriate compactness assumptions, the empirical distribution converges weakly to a limiting distribution as , resulting in the following kinetic equation:
2.2 In two dimensions
The agents are now described by their positions within domain and fear levels . The time interval of interest is . Each agent moves with speed equal to their fear level and with a prescribed direction . This could be, e.g., the direction to the nearest exit if we are considering an evacuation scenario. So, we obtain the following microscopic model:for , where , the average fear level for agent , defined as: being an interaction kernel. We consider again the kernel in Equation 3.
For the derivation of the kinetic model, we follow the same steps used in 1D. We denote the empirical distribution density as:where is as in Equation 6 and is the Dirac delta function. We multiply Equation 7 by a test function , integrate with respect to space and fear level and differentiate in time:
Since , then
By plugging Equation 6 into the above equation, we obtain:
Using Equation 7, we get
We have:where
Similarly, we have:where
Using Equations 9, 10, the second term at the right-hand side of Equation 8 can be written as:
Thus, Equation 8 becomes
Integrating the terms on the right-hand side of this equation by parts, assuming that vanishes at the boundaries, and rearranging the terms, we get:
Since is an arbitrary test function, this leads to:which can be rewritten as:where the local average fear level is given by:
Under appropriate compactness assumptions, the empirical distribution converges weakly to a limiting distribution as , resulting in the following kinetic equation:
2.3 Discretization of the 2D models
We use the explicit Euler scheme for the system of ordinary differential Equation 8. Let be a time step and , . Denote by the values of the position and the fear level of agent at time . The values at time are given bywith
Next, we consider the discretization of Equation 11 in space, time, and fear level using a finite difference method. We assume that , for given . We partition using mesh size along the axis and along the axis, so to obtain cells:where the centers of the cells are given by , for and , where , and . The fear level varies in interval , with indicating no stress and indicating the highest stress level. Interval [0,1] is partitioned into intervals of equal length :where the centers of the cells are given by , for , with . These partitions of and [0,1] induce a partition of the 3D domain into 3D cells.
Let be the time interval of interest and set the time step toto satisfy the Courant–Friedrichs–Lewy (CFL) condition (see, e.g. [26]). Further, we denote , , and for any given quantity , its approximation at a specific time point is denoted with .
The finite difference method produces an approximation of the average over the cell
To write the finite difference discretization of Equation 11, we introduce
Additionally, let be slope limiter function. For the numerical results in Section 4, we will consider and compare two slope limiters: the van Leer limiterand the minmod slope limiter
See, e.g., [26] for more details.
A finite difference scheme for Equation 11 is as follows: At time step , find for , and , such that:where the flux in variable iswithand the flux in variable iswith
The flux and flux limiter in iswith
The subscript in Equation 15 is if and if .
Scheme Equation 14 is second-order in the fear level variable when the solution is smooth. If one sets , for all in Equation 14, a first-order upwind numerical scheme [26] is obtained.
Finally, we remark that by setting to 0 all terms in one can easily write down the 1D counterparts of Equations 12, 14.
3 Hybrid microscopic-to-mesoscopic approach
The main idea of the hybrid scheme is to use different models depending on the local crowd density. In regions of the domain where the density is lower than a critical threshold , we will use the agent-based formulation and in the regions where the density is larger than a certain threshold, we will use the mesoscopic or kinetic formulation. For clarity, we present the hybrid scheme in one dimension first, and then extend it to two dimensions.
3.1 In one dimension
For simplicity of description, we assume that the kinetic domain, i.e., the part of the computational domain where the kinetic model holds, is an interval , where and belong to the spatial mesh and can change from one time step to another. In other words, the locations of the interfaces between the agent-based and the kinetic regimes change over time. Because of this, a cell that belongs to the kinetic domain at time might be located in the agent-based domain at the next time . This happens if or . Additionally, the agents’ trajectories can cross the interface boundaries of the kinetic domain. Thus, our algorithm must describe how to transfer information from the kinetic formulation to the agent-based formulation, and vice versa. In general, the cells that belong to the kinetic domain do not have to be adjacent, i.e., we can have multiple microscopic-kinetic interfaces. So, in general, would denote the collection of cells in the kinetic regime.
In order to conserve the total “mass” of people, i.e., total number of people in the system when both probabilistic and deterministic models are used, we introduce a generalized microscopic model in which agents can have different, and in general, non-integer masses. To describe how it works, let be a set of integers representing the labels of agents present in the model at time step and let denote the position, mass, and fear level of agent . Position and fear level come from the discrete model, i.e., the 1D counterpart of Equation 12, and initially the masses of all agents are equal to 1, as mentioned in Section 2.1. Let be the discrete probability density function at time computed with the 1D counterpart of Equation 14. We assume that is defined at all mesh points regardless of whether they belong to or not, with for . Then, initially the total density of the system at the mesh points is given bywhere is an approximation of the Dirac delta function . The values of are used to determine if mesh point belong to the kinetic domain at next time step , i.e., . Specifically, if , then will belong to . Let and be the position of these interfaces at time .
When agents cross an interface and pass from the microscopic domain to the kinetic domain, we have an inflow of mass into the kinetic domain. We handle this as follows: if agent crosses the microscopic-to-kinetic interface, this agent is added to the probability distribution function :where is a piece-wise approximation of the delta function:
Note that in this case we choose a simpler approximation for the Dirac delta to have a more direct control on the support. The updated in Equation 17 is plugged into Equation 14 to find and index is removed from set , i.e., .
Next, we consider crossing the interface in the opposite direction, i.e., from the kinetic domain to the microscopic domain. Suppose , i.e., the microscopic-to-kinetic interface moves to the right from time to . Let be the set of indexes such that . Then, the mass of people that occupy a position in at time is given by:
Since at time , interval belongs to the microscopic domain (and hence a deterministic model holds there), we create a new agent that occupies the midpoint of . For conservation of mass, agent has mass equal to the mass in , and fear level equal to the average in . This means:
Then, index gets added to to form . This procedure is repeated at the other (kinetic-to-microscopic) interface and, in general, at all other interfaces.
To avoid loss or creation of mass (i.e., people), if Equation 19 is smaller than 1, a new agent is not created and the position of the interface is not changed. However, over a certain number of time steps the mass leaving the kinetic domain for the microscopic domain will amount to 1 or more. We note that in our simplified model, all pedestrians have non-negative fear levels and move with the assigned walking direction, which is aligned with the positive -axis. Thus, in general, an out-flux of mass happens at the interface positioned at . Since the walking speed equals the fear level, the mass of agents with fear level per unit of time flowing out of the kinetic region through the interface at equals
Note that, since is speed, is length and thus Equation 20 is indeed a mass. Over the time interval , the number of agents with fear level that leave the kinetic domain equalsand the total number of agents leaving the kinetic domain equals
So, when this mass is greater than 1, we create a new agent with parameters
In conclusion, at time for all agents in the microscopic domain, position and fear level are given by the 1D counterpart of Equation 12, and the mass is computed as explained above. For the cells that belong to the kinetic domain , the values of the probability distribution function are determined using the 1D counterpart of Equation 14. For both models, the average fear level is computed bywhere is either the position of agent , i.e., , or mesh point . Note that Equation 21 computes the average regarding of whether a neighborhood of agent belongs entirely to the microscopic domain, the kinetic domain, or is across the two domains.
3.2 In two dimensions
For simplicity of description, we assume that the kinetic domain, i.e., the part of the computational domain where the kinetic model holds, is a rectanglewhere and belong to the spatial mesh. Again, we denote by the set of integers representing the labels of the agents present in the microscopic domain at time step . The generic agent is characterized by their mass , which is initially equal to 1, position and fear level from Equation 12. We assume that the discrete probability density from Equation 14 is defined at all mesh points regardless of whether they below to or not, with if . Then, at time the total density of the system at the mesh points is approximated bywhere is defined component-wise in Equation 16.
We use to determine the discrete positions of the interfaces between microscopic and kinetic domains at the next time step :
Then, the kinetic domain at time is . See Figure 2.
FIGURE 2
When an agent crosses an interface and pass from the microscopic domain to the kinetic domain, the agent is added to the probability distribution function :where is defined component-wise in Equation 18. This updated is plugged into Equation 14 and index is removed from set , so that .
Next, we consider crossing the interface in the opposite direction, i.e., from the kinetic domain to the microscopic domain. Suppose that , and , , as in Figure 2. Other cases are treated analogously. Then, the following region becomes part of the microscopic domain:
This is the union of the shaded rectangles in Figure 2.
We explain the procedure we follow for , with the understanding that the same procedure is repeated for . Let be the set of indexes such that and be the set of indexes such that . The mass of people that occupy a position in , which belongs to the microscopic domain at time , is given by:
Hence, we create a new agent at the midpoint of , with mass equal Equation 23 for conservation and fear level equal to the average fear level in :
Index gets added to to form . Just like in 1D, if Equation 23 is smaller than 1, a new agent is not created and the position of the interfaces are not changed. However, over a certain number of time steps, the mass leaving the kinetic domain for the microscopic domain will be equal or exceed 1, and then the interface changes.
Let us consider interface . For simplicity, we assume that corresponds to (mesh point) and that is positive for all in . The total number of agents leaving the kinetic domain through this interface over time interval equals
When this mass is greater than or equal 1, we create a new agent with parametersand is randomly chosen from the interval . A similar computation is performed at the interface .
In conclusion, at time for all agents in the microscopic domain, position and fear level are given by Equation 12, the mass is computed as explained above, and the average fear level qi*,n is computed by
For the cells that belong to the kinetic domain , the values of the probability distribution function are given by scheme Equation 14. At a mesh point of the kinetic domain, the average fear level qi,j*,n is computed by:
4 Numerical results
We assess our hybrid method through test cases in one dimension, presented in Section 4.1, and two dimensions, reported in Section 4.2.
4.1 Tests in one dimension
We consider agents that are initially uniformly distributed on the computational domain , so we have an initial uniform density of 10. We initially set the fear level of the agents in the interval to be equal to 1 and the fear level of the agents in the interval to be 0. This is an unlikely but challenging case where half of the people in the domain have no fear, and thus do not walk, while the remaining half behind them is panicking and hence getting in motion with maximum walking speed. We recall that everyone walks with a direction aligned with the positive -axis. We set the strength of fear contagion and the interaction radius in Equation 3. The time interval of interest is [0,4].
The agent based dynamics is obtained via the 1D counterpart of Equation 12 with . Figure 3 shows the plots of the agent densityand the average fear levelwith defined in Equation 16 and , over a part of the computational domain at different times. We see that a region of high density develops near the center of the domain and it propagates to the right.
FIGURE 3
Next, we apply the hybrid microscopic-kinetic methods from Section 3.1, with critical density . As can be seen from Figure 3, initially we have everywhere. Thus, our approach starts with the agent-based model everywhere in . Then, as time passes, becomes larger than over part of the domain. This is when our approach becomes indeed hybrid, with the kinetic domain moving to the right as time evolves. We consider 4 different meshes with mesh size and the time step is set according to Equation 13. We start with the first order numerical scheme, i.e., no limiter is used. Figure 4 compares the crowd density obtained by the hybrid scheme with the density given by the microscopic model at the end of the time interval of interest, i.e., at . We see that for coarser meshes the results given by hybrid model and the microscopic model differ significantly. However, when the mesh is refined the densities computed by the 2 methods become closer and closer, as one would expect.
FIGURE 4
Now, let us consider the same meshes and time steps as above, with the second-order scheme. Figures 5, 6 show the crowd density when the van Leer limiter and the minmod limiter are applied, respectively, and the comparison with the crowd density given by the microscopic model. Again, we see great agreement as the mesh is refined.
FIGURE 5
FIGURE 6
To quantify the comparisons in Figures 4–6 and Tables 1, 2 report the corresponding differences in and norms. From these tables, we see that the second-order schemes perform marginally better than the first-order scheme and are comparable to each other. We notice that this is not a convergence test, as the solution given by the microscopic model cannot be considered the exact solution for the kinetic model. It is just showing that when the crowd is rather dense and the mesh is fine enough, the microscopic model and the hybrid approach are in agreement.
TABLE 1
| Method | ||||
|---|---|---|---|---|
| 1st order | 36.1 (3.6%) | 21.1 (2.1%) | 7.9 (0.8%) | 4.2 (0.4%) |
| van Leer | 35.3 (3.5%) | 24.5 (2.5%) | 6.8 (0.7%) | 3.5 (0.4%) |
| minmod | 32.2 (3.22%) | 21.8 (2.2%) | 6.8 (0.7%) | 3.5 (0.4%) |
Difference between the solution given by the microscopic model and the solution given by the hybrid approaches in norm at time , with relative difference included in parenthesis.
TABLE 2
| Method | ||||
|---|---|---|---|---|
| 1st order | 21.9 (20.3%) | 16.1 (14.9%) | 6.6 (6.1%) | 3.6 (3.3%) |
| van Leer | 20.1 (18.6%) | 22.7 (21%) | 5.8 (5.4%) | 2.9 (2.7%) |
| minmod | 17.9 (16.6%) | 18.2 (16.9%) | 6.1 (5.4%) | 2.9 (2.7%) |
Difference between the solution given by the microscopic model and the solution given by the hybrid approaches in norm at time , with relative difference included in parenthesis.
Finally, we present some results that show the importance of model parameter . Recall that in Equation 1 is the interaction strength, i.e., how much an agent is willing to adjust their fear level, and hence walking speed, to equilibrate it with the local average. The results presented thus far were obtained with . Figure 7 shows the crowd density given by the hybrid method for with the minmod limiter. All the other settings are the same used so far. By comparing Figure 7 with Figure 3 (left), we see that the crowd density is very sensitive to the value of , both in terms of magnitude of the peaks and overall distribution.
FIGURE 7
While constant is convenient, it is not very realistic since different people react differently to a spreading emotion. Figure 8 reports the results obtained with the microscopic model and random values of assigned to the agents, shown in the leftmost panel. Notice how the heterogeneity in the values of leads to local oscillations in the density (central panel) and average fear level (rightmost panel).
FIGURE 8
4.2 Tests in two dimensions
We consider agents initially placed uniformly in computational domain , so that the initial density is equal to 9 everywhere in the domain. The initial fear level is set to 1 inside the circle of radius 3 centered at the origin of the axes, while everywhere else it is set to 0. This set up is meant to reproduce a case where a group of people inside a uniformly dense crowd gets scared (by, e.g., the onset of a fire) and panics, while all other people are not aware of the threat and thus do not move. We assume that all moving agents choose walking direction , possibly the direction where an exit is located. We set and in Equation 3. The time interval of interest is [0,5].
The agent based dynamics is obtained by using Equation 12 with . Figure 9 shows the plots of the agent densityand the average fear levelat , with defined component-wise in Equation 16 and .
FIGURE 9
Next, we apply the hybrid microscopic-kinetic methods from Section 3.2, with critical density . Since initially everywhere in the domain, our scheme starts with the microscopic model everywhere. As time passes, becomes larger than over parts of the domain and our approach becomes truly hybrid. We consider 3 different meshes with mesh size . The time step is set according to Equation 13. We start with the first order numerical scheme, i.e., no limiter is used. Figure 10 compares the crowd density obtained by the microscopic model (top left panel) with the density computed by the hybrid approach with the three meshes under consideration at time . As in the one-dimensional case, the results from the hybrid and microscopic models differ significantly when a coarse mesh is used. However, as the mesh is refined, the densities computed by the two methods become increasingly closer.
FIGURE 10
We repeat the comparison using a second-order schemes, with the same meshes and time step as above. Figures 11, 12 report the microscopic vs. hybrid comparison in terms of crowd density when the van Leer limiter and the minmod limiter are applied, respectively. Again, we see that the agreement improves as the computational mesh is refined.
FIGURE 11
FIGURE 12
To quantify the level of agreement shown in Figures 10–12, we report in Tables 3, 4 the difference in absolute value between the solution given by the microscopic model and the solution given by the hybrid approach in the and norms. In this comparison, which we recall is not a convergence test, we do not observe a significant difference when using limiters compared to the first-order scheme. Like in the 1D case, we see that when the crowd is rather dense and the mesh is fine enough, the microscopic model and the hybrid approach are in agreement.
TABLE 3
| Method | |||
|---|---|---|---|
| 1st order | 90 (10%) | 46.4 (5%) | 29.2 (3%) |
| van Leer | 90.7 (10.4%) | 47.8 (5.1%) | 27.7 (3.2%) |
| minmod | 89.5 (10%) | 47.5 (5.4%) | 25.9 (2.9%) |
Difference in absolute value between the solution given by the microscopic model and the solution given by the hybrid approaches in the norm at time , with relative difference included in parenthesis.
TABLE 4
| Method | |||
|---|---|---|---|
| 1st order | 16.9 (37%) | 10.6 (23%) | 6.2 (13%) |
| van Leer | 17.4 (37.7%) | 11.2 (24%) | 6.2 (13%) |
| minmod | 17.1 (36.9%) | 11 (23.8%) | 5.9 (12.8%) |
Difference in absolute value between the solution given by the microscopic model and the solution given by the hybrid approaches in the norm at time , with relative difference included in parenthesis.
We conclude by considering the more realistic case where the value of the interaction strength varies from one agent to the other. Like for the corresponding 1D test, we assign a random value of in interval (0,1) to each agent. See the assigned values in Figure 13. Figure 14 reports the crowd density and average fear level at obtained with the microscopic model and these random values of . In both quantities, we observe local oscillations that were absent when the value of was constant in space (see Figure 9).
FIGURE 13
FIGURE 14
5 Conclusions and future perspectives
We developed and numerically analyzed a hybrid modeling approach for simulating crowd dynamics under the influence of propagating fear in walking venues with mixed crowd densities. The hybrid nature is related to the fact that our approach uses a microscopic (agent-based) model in low-density regions, while in higher-density areas it adopts a mesoscopic (kinetic) model, derived from the microscopic model via a mean-field limit approach. The microscopic model captures individual variability and local interactions, while the mesoscopic model efficiently handles collective behaviors in denser crowds. By coupling the two models with appropriate interface conditions, we take advantages of their strengths while mitigating their respective limitations. We ensured consistency between the two models through a mean-field limit and proposed interface conditions to enable seamless coupling of regions of varying density.
We performed tests in one and two dimensional spatial domains. We showed that our hybrid approach produces results in agreement with those provided by the microscopic model when the crowd is sufficiently dense and the computational mesh is fine enough. This validates the effectiveness of our hybrid strategy in accurately capturing crowd dynamics across varying density regimes. Additionally, our tests demonstrate the sensitivity of the computed crowd density and average fear level on model parameter , which controls the strength of emotional interactions. Thus, one would have to find a way to carefully set this parameter to realistically simulate crowd responses in fear-driven scenarios.
In order to focus on the behavioral features of crowd dynamics, we have introduced a strong simplification in this work: the walking direction is prescribed. This choice represents the main limitation of the hybrid model presented in this paper. In reality, each individual in a crowd has their own purposes and abilities, which lead to the selection of a walking strategy, i.e., a trajectory to follow in order to reach the desired target(s) and a speed to move along such trajectory. The walking strategy is influenced by the emotional state (considered in this paper) and the interactions with other people in the crowd and the environment, which includes walls, exits, and obstacles (not considered in this paper). A relatively easy way to account for such interactions is with tools of game theory, as we have done in [18, 23]. However, the resulting more realistic model presents difficulties at the numerical level related to need to contain high computational costs. One possible way to meet this need is through splitting algorithms, which break the model into a set of subproblems that are easier to solve and for which practical algorithms are readily available. The inclusion of interactions with people and the environment and the development of suitable splitting algorithms will be object of future work.
Statements
Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors upon request, without undue reservation.
Author contributions
IP: Validation, Writing – original draft, Visualization, Methodology, Formal Analysis, Investigation, Software, Writing – review and editing. AQ: Conceptualization, Investigation, Methodology, Supervision, Writing – review and editing, Writing – original draft.
Funding
The author(s) declare that no financial support was received for the research and/or publication of this article.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
References
1.
HelbingDJohanssonAAl-AbideenHZ. Dynamics of crowd disasters: an empirical study. Phys Rev E (2007) 75:046109. 10.1103/PhysRevE.75.046109
2.
HelbingDJohanssonA. Pedestrian, crowd, and evacuation dynamics. New York, NY: Springer New York (2009). p. 1–28. 10.1007/978-3-642-27737-5_382-5
3.
WijermansNConradoCvan SteenMMartellaCLiJ. A landscape of crowd-management support: an integrative approach. Saf Sci (2016) 86:142–64. 10.1016/j.ssci.2016.02.027
4.
HaghaniMSarviM. Social dynamics in emergency evacuations: disentangling crowd’s attraction and repulsion effects. Physica A: Stat Mech its Appl (2017) 475:24–34. 10.1016/j.physa.2017.02.010
5.
AylajBBellomoNGibelliLRealiA. A unified multiscale vision of behavioral crowds. Math Models Methods Appl Sci (2020) 30(1):1–22. 10.1142/s0218202520500013
6.
BellomoNGibelliLQuainiARealiA. Towards a mathematical theory of behavioral human crowds. Math Models Methods Appl Sci (2022) 32(02):321–58. 10.1142/s0218202522500087
7.
BellomoNLiaoJQuainiARussoLSiettosC. Human behavioral crowds review, critical analysis and research perspectives. Math Models Methods Appl Sci (2023) 33(08):1611–59. 10.1142/s0218202523500379
8.
HelbingDFarkasIVicsekT. Simulating dynamical features of escape panic. Natures (2000) 407:487–90. 10.1038/35035023
9.
WangJZhangLShiQYangPHuX. Modeling and simulating for congestion pedestrian evacuation with panic. Physica A: Stat Mech its Appl (2015) 428:396–409. 10.1016/j.physa.2015.01.057
10.
ZouYXieJWangB. Evacuation of pedestrians with two motion modes for panic system. PLoS One (2016) 11(4):e0153388. 10.1371/journal.pone.0153388
11.
RahmatiYTalebpourA. Learning-based game theoretical framework for modeling pedestrian motion. Phys Rev E (2018) 98:032312. 10.1103/PhysRevE.98.032312
12.
XuMXieXLvPNiuJWangHLiCet alCrowd behavior simulation with emotional contagion in unexpected multihazard situations. IEEE Trans Syst Man, Cybernetics: Syst (2021) 51(3):1–15. 10.1109/TSMC.2019.2899047
13.
van HaeringenEGerritsenCHindriksK. Emotion contagion in agent-based simulations of crowds: a systematic review. Auton Agent Multi-agent Syst (2023) 37:6. 10.1007/s10458-022-09589-z
14.
BellomoNBellouquidA. On multiscale models of pedestrian crowds from mesoscopic to macroscopic. Commun Math Sci (2015) 13:1649–64. 10.4310/cms.2015.v13.n7.a1
15.
BellomoNBellouquidAGibelliLOutadaN. A quest towards a mathematical theory of living systems. In: Modeling and simulation in science, engineering and technology. Birkhäuser (2017). 10.1007/978-3-319-57436-3
16.
BellomoNGibelliL. Toward a mathematical theory of behavioral-social dynamics for pedestrian crowds. Math Models Methods Appl Sci (2015) 25(13):2417–37. 10.1142/S0218202515400138
17.
BellomoNGibelliLOutadaN. On the interplay between behavioral dynamics and social interactions in human crowds. Kinetic Relat Models (2019) 12(2):397–409. 10.3934/krm.2019017
18.
KimDQuainiA. A kinetic theory approach to model pedestrian dynamics in bounded domains with obstacles. Kinetic Relat Models (2019) 12(6):1273–96. 10.3934/krm.2019049
19.
WangLShortMBBertozziAL. Efficient numerical methods for multiscale crowd dynamics with emotional contagion. Math Models Methods Appl Sci (2017) 27(1):205–30. 10.1142/S0218202517400073
20.
LlorcaDHaghaniMCristianiEBodeNBoltesMCorbettaA. Panic, irrationality, and herding: three ambiguous terms in crowd dynamics research. J Adv Transportation (2019) 2019:1–58. 10.1155/2019/9267643
21.
BellomoNBellouquidAKnopoffD. From the microscale to collective crowd dynamics. SIAM Multiscale Model and Simulation (2013) 11(3):943–63. 10.1137/130904569
22.
AgnelliJPColasuonnoFKnopoffD. A kinetic theory approach to the dynamics of crowd evacuation from bounded domains. Math Models Methods Appl Sci (2015) 25(01):109–29. 10.1142/S0218202515500049
23.
KimDO’ConnellKOttWQuainiA. A kinetic theory approach for 2D crowd dynamics with emotional contagion. Math Models Methods Appl Sci (2021) 31(06):1137–62. 10.1142/S0218202521400030
24.
KimDLabateDMilyKQuainiA. Data-driven learning to enhance a kinetic model of distressed crowd dynamics. Math Models Methods Appl Sci (2025) 35(07):1609–36. 10.1142/S021820252540007X
25.
BakhdilNEl MousaouiAHakimA. A kinetic BGK model for pedestrian dynamics accounting for anxiety conditions. Symmetry (2025) 17(1):19. 10.3390/sym17010019
26.
LeVequeRJ. Numerical methods for conservation laws. In: Lectures in mathematics, ETH Zürich. Basel: Birkhäuser Verlag (1990).
Summary
Keywords
crowd dynamics, fear propagation, kinetic model, microscopic model, multiscale model, complex systems
Citation
Perepelitsa I and Quaini A (2025) Coupling microscopic and mesoscopic models for crowd dynamics with emotional contagion. Front. Phys. 13:1644470. doi: 10.3389/fphy.2025.1644470
Received
10 June 2025
Accepted
14 July 2025
Published
20 August 2025
Volume
13 - 2025
Edited by
Ryosuke Yano, Tokio Marine dR Co., Ltd., Japan
Reviewed by
Nouamane Bakhdil, Cadi Ayyad University, Morocco
Nisrine Outada, Cadi Ayyad University, Morocco
Updates
Copyright
© 2025 Perepelitsa and Quaini.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Annalisa Quaini, aquaini@central.uh.edu
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.