Mean-field and Monte Carlo analysis of multi-species agent dynamics

We propose a mean-field (MF) approximation as a recurrence relation governing the dynamics of $m$ species of particles on a square lattice. We simultaneously perform Monte Carlo (MC) simulations under identical initial conditions to emulate the intricate motion observed in environments such as subway corridors and scramble crossings in large cities. Each species moves according to transition probabilities influenced by its respective static floor field and the state of neighboring cells. To illustrate the methodology, we analyze statistical fluctuations in the spatial distribution for $m=1$, $m=2$, and $m=4$ and for different regimes of average density and biased movement. A numerical comparison is conducted to determine the best agreement between the MC simulations and the MF approximation considering a renormalization exponent $\beta$ that optimizes the fit between methods. Finally, we report a phenomenon we term “Gaussian-to-Gaussian” behavior, in which an initially normal distribution of particles becomes distorted due to interactions among same and opposing species, passes through a transient regime, and eventually returns to a Gaussian-like profile in the steady state, after multiple rounds of motion under periodic boundary conditions.

1 Introduction

Pedestrian dynamics is one example of complex systems of self-driven agents in the realm of macro-scale physics, which give rise to a variety of emerging phenomena. Several theoretical and experimental studies have been conducted on crowd dynamics, focusing on topics ranging from city planning [1] and disaster prevention to the understanding of complex human behavior [2, 3].

Some places with crowd formation as subway corridors [4, 5], nightclub dynamics [6, 7], and scramble crossings [8] in large cities, are examples of these complex systems that exhibit collective behavior. Simply put, it features different groups of individuals aiming to reach various common destinations, for instance, the entrance/exit of a train station, reaching/leaving the bar to buy beverages, or crossing to get to another corner. Thus, the spatial and temporal context in question may involve a single group of agents moving towards a common goal, or multiple groups moving along confronting trajectories, which rapidly increases the complexity of the dynamics.

The problem of two groups of particles in counterflow has been extensively studied across various contexts, ranging from microscopic scales—such as the dynamics of charged colloids [9]—to macroscopic phenomena, such as pedestrian flows in corridors [10, 11], and even more complex situations like evacuation [12] scenarios. Notably, even in single-species models [13], intriguing phenomena such as condensate formation emerge, illustrating the richness of these systems.

An interesting approach to describe the hard-body dynamics of macroscopic systems is the lattice-gas modelling [14]. In our recent work, we used the dynamics of a lattice gas-type interaction to describe the motion of pedestrians inside nightclubs [6, 7] and in a four-way crossing walk [8], the so-called scramble crossing. In such works, we have shown that the intrinsic high correlation of this approach reflects the model’s sensitivity to jamming and condensate formation.

In that sense, Monte Carlo (MC) simulations are a natural method for studying such models, as they allow for the definition and implementation of transition rules to simulate different evolutions from specified initial conditions. These systems can also be interpreted in a broader context as multi-agent systems. Additionally, there are approaches that model these dynamics through differential equations based on Newton’s laws, within the so-called social force framework [12, 15, 16].

In that direction, the mean-field (MF) approach in particle dynamics has been extensively studied, particularly in the context of lattice gas models. These models were originally proposed to describe the behavior of specific types of fluids, such as solvents [17], electrochemical cells [18], and other systems that can be characterized by Hamiltonians derived from free energy potentials [19]. In 2015, we introduced a related model [20] based on a two-species particle framework to describe pedestrian counterflow. Unlike classical exclusion principles, our model employed density-dependent exclusion rules, where the probability of occupation varied with the local density of particles. The system’s evolution was governed by a set of coupled partial differential equations (PDEs). Other variations have been explored to model deterministic dynamics [21, 22], and further extensions of our framework have incorporated concepts from Fermi-Dirac statistics [23].

In this work, we propose a general framework for a MF approximation of the lattice gas dynamics of m-species of particles on a square lattice. We derive a generalized non-linear coupled PDE obtained by the MF approximation and compare it to MC simulations, by using a set of initial conditions for the cases $m=1$, $m=2$, and $m=4$, where each type of particle is guided by its own static floor field [24] towards its target. We study the statistical properties of the distributions of particles of the species and investigate what are the circumstances and range of parameters that make the model a suitable and efficient option for MC simulations.

Our results show excellent agreement between the PDE predictions and MC simulations for low-density regimes. However, at higher densities, jamming effects reduce the accuracy of this agreement. We present a detailed numerical study to delineate the conditions under which our method remains valid compared to the computationally more expensive MC simulations.

In the next section, we describe the method for deriving the PDEs based on a model where particles exhibit preferential directed motion combined with random movements dependent on the occupation of neighboring sites. Through a mean-field approximation, we successfully derive a generalized PDE that describes scenarios involving an arbitrary number m of species.

Subsequently, in the third section, we present our results comparing the spatial distribution of particles on a bidimensional lattice for different model parameters and initial particle concentrations. We study three different cases —$m=1$, $m=2$, and $m=4$—, analyzing the evolution of key distribution features over time under periodic boundary conditions. These features reflect the interactions between particles across the various scenarios. Additional details—including equation derivations, supplementary plots, and a discussion of numerical instabilities related to the PDEs—are provided in the supplemental material accompanying this text.

2 Methods

Our model consists of a system of $n$ particles that move along an underlying square lattice of side $l$. Each particle in the system belong to one of $m$ species and can hop to one of its first neighboring cells at each time step according to a set of transition probabilities which depends on its species “social fields” (static floor field), denoted as ${\mathbf{u}}^{(k)}={\mathbf{u}}^{(k)}(\mathbf{r})$, where $k=1,\dots ,m$. For instance, a given particle belonging to the species $q$ at cell $\mathbf{r}=(i,j)$ hops to one of its neighboring cells $\mathbf{r}\prime =(i\prime ,j\prime )$ in the unit time interval $t\to t+1$ according to the transition probability density

$P{r}_{\mathbf{r}\to {\mathbf{r}}^{\prime }}^{\left(q\right)}=p+\alpha \mathrm{\Delta }\mathbf{r}•{\mathbf{û}}^{\left(q\right)},(1)$

where $p$ is a constant probability to hop to one of the four neighboring cells, $\alpha$ is the coefficient that measure the bias towards the direction of its species static floor field, ${\mathbf{û}}^{(q)}\equiv {\mathbf{u}}^{(q)}/\Vert {\mathbf{u}}^{(q)}\Vert$ is the normalized static floor field of species $q$, and $\mathrm{\Delta }\mathbf{r}\equiv {\mathbf{r}}^{\prime }-\mathbf{r}$ is the displacement vector of the calculated transition probability. As one may expect, normalization constraints reflect on the possible values of $p$, which in a Cartesian regular lattice is $0<p\le 1/4$, whilst $\alpha$ satisfies the condition $0<\alpha \le p$.

It is important to note that, in our model, each particle perceives only the static floor field specific to its species, evaluated at its current position. Consequently, the inner product defined in Equation 1 quantifies the contribution of each possible transition probability to the particle’s movement toward its target.

As one might anticipate, Equation 1 is the general form of the four possible transition probabilities each of which have a different corresponding displacement vector that can assume the values $\mathrm{\Delta }\mathbf{r}={\mathbf{e}}_x$, $-{\mathbf{e}}_x$, ${\mathbf{e}}_y$, and $-{\mathbf{e}}_y$, with ${\mathbf{e}}_x$ and ${\mathbf{e}}_y$ denoting the unit vectors of our reference frame as depicted in Figure 1. Thus, to ensure the normalization constraint, the probability that a particle remains in its current cell is expressed as:

$P{r}_{\mathbf{r}\to \mathbf{r}}^{\left(q\right)}=1-{\displaystyle \sum _{⟨{\mathbf{r}}^{\prime }⟩}}P{r}_{\mathbf{r}\to {\mathbf{r}}^{\prime }}^{\left(q\right)}=1-4p\text{,}(2)$

where the angle brackets $⟨{\mathbf{r}}^{\prime }⟩$ indicate that the summation is carried out solely over the nearest neighbours, with $i^{\prime }=i\pm 1$ and $j^{\prime }=j\pm 1$.

Figure 1

Figure 1. Cartesian lattice depiction of a likely transition to happen of a particle at position $\mathbf{r}$ to a neighbouring cell ${\mathbf{r}}^{\prime }=\mathbf{r}+\mathrm{\Delta }\mathbf{r}$.

It is important to note that these transition probabilities define the intended particle movement, conditional on the target cell being empty, as each cell can accommodate at most one particle. Thus, the probability represents what should occur, not what will necessarily occur.

2.1 Lattice gas dynamics

We focus our study on an approximate regime of a lattice gas model. The key feature, as previously noted, is the exclusion principle, which prevents particles from overlapping or occupying the same cell, regardless of the specific rules governing the transition probabilities. This characteristic makes it particularly challenging to derive an approximate recurrence relation or equation of motion, as the occupation states of neighboring lattice cells are highly correlated. Nevertheless, we develop a recurrence relation for the cell occupation that remains consistent with the asynchronous updating scheme characteristic of lattice gas dynamics with nearest-neighbor interactions. Accordingly, we propose the following recurrence relation for a given species q:

$\begin{array}{cc}{\rho }_{t+1}^{\left(q\right)}\left(\mathbf{r}\right)\hfill & =\left(1-{\rho }_t\left(\mathbf{r}\right)\right){\displaystyle \sum _{\langle {\mathbf{r}}^{\prime }\rangle }}{\rho }_t^{\left(q\right)}\left({\mathbf{r}}^{\prime }\right)P{r}_{{\mathbf{r}}^{\prime }\to \mathbf{r}}^{\left(q\right)}\hfill \\ \hfill & +{\rho }_t^{\left(q\right)}\left(\mathbf{r}\right)\left(P{r}_{\mathbf{r}\to \mathbf{r}}^{\left(q\right)}+{\displaystyle \sum _{\langle {\mathbf{r}}^{\prime }\rangle }}{\rho }_t\left({\mathbf{r}}^{\prime }\right)P{r}_{\mathbf{r}\to {\mathbf{r}}^{\prime }}^{\left(q\right)}\right),\hfill \end{array}(3)$

where ${\rho }_t(\mathbf{r})\equiv {\displaystyle \sum _{k=1}^m}{\rho }_t^{(k)}(\mathbf{r})$ is the total density of particles at a given cell $\mathbf{r}$ at time step $t$. The proposed equation decomposes the update process into two key contributions: (1) particle inflow into the cell and (2) particle persistence or outflow from the cell. In the lattice gas dynamics, the occupation state is binary, with ${\rho }_t^{(q)}(\mathbf{r})=1$ if the cell is occupied and ${\rho }_t^{(q)}(\mathbf{r})=0$ if it is vacant. However, in the mean-field (MF) approximation, ${\rho }_t(\mathbf{r})$ represents an average density field and does not directly reflect the actual number of particles $n$ in the system. To bridge this gap and improve consistency between the MF and MC descriptions, we generalized the normalization condition of the MF recurrence relation. Specifically, we modified the global constraint so that the sum over all lattice sites scales as a power of $n$:

${\displaystyle \sum _{\mathbf{r}}}{\rho }_t\left(\mathbf{r}\right)=n^{\beta },$

where $\beta$ is a tunable exponent. The introduction of $\beta$ allows us to adjust the MF approximation to better reproduce the behavior observed in MC simulations. By treating $\beta$ as a fitting parameter, we aim to compensate for the differences between the averaged MF treatment and the discrete, stochastic nature of the MC dynamics.

The first term on the right-hand side of Equation 3, the factor $\left(1-{\rho }_t(\mathbf{r})\right)$, reflects the hard-body exclusion principle, which dictates that a particle can only occupy an empty cell. This factor multiplies the sum of possible sources of inflow of particles of species $q$ from the neighboring cells of $\mathbf{r}$, i.e., the notation $⟨{\mathbf{r}}^{\prime }⟩$. Each neighboring source of particles depends on the density of the species $q$ at the neighboring site ${\rho }_t^{(q)}({\mathbf{r}}^{\prime })$ and the transition probability $P{r}_{{\mathbf{r}}^{\prime }\to \mathbf{r}}^{(q)}$, reflecting its likelihood to hop from ${\mathbf{r}}^{\prime }$ to $\mathbf{r}$. The second term on the right-hand side of Equation 3 describes the contribution of particles already present in $\mathbf{r}$. The first factor, ${\rho }_t^{(q)}(\mathbf{r})$, accounts for the actual occupation of it, whilst the factor inside the large brackets corresponds to a composition of two possibilities. The first possibility is the probability $P{r}_{\mathbf{r}\to \mathbf{r}}^{(q)}$ of the particle at $\mathbf{r}$ remaining in place, which contrasts with the second possibility that stands as the probability of the particle at $\mathbf{r}$ trying to move to an occupied neighboring cell ${\mathbf{r}}^{\prime }$.

With this framework, we explicitly separate the exclusion interaction, reflected in the terms $\left(1-{\rho }_t(\mathbf{r})\right)$ and ${\rho }_t(\mathbf{r})$, from the transition probabilities, which can incorporate additional interaction mechanisms such as long-range fields or cooperative dynamics. As an example, Equation 3 can reproduce even deterministic models of particles, such as proposed in [21], which could allow us to describe the model proposed by Cividini et al., for instance, by a proper choice of parameters such as fixing $m=2$ with one species moving eastbound and the other species moving northbound with transition probabilities equal to one towards their respective target directions.

We can rewrite Equation 3 as follows

$\begin{array}{cc}{\rho }_{t+1}^{\left(q\right)}\left(\mathbf{r}\right)\hfill & ={\displaystyle \sum _{⟨{\mathbf{r}}^{\prime }⟩}}\left[\left({\omega }_t^{\left(q\right)}\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)-{\omega }_t^{\left(q\right)}\left({\mathbf{r}}^{\prime },\mathbf{r}\right)\right)\right.\hfill \\ \hfill & -\left.\left({\rho }_t\left(\mathbf{r}\right){\omega }_t^{\left(q\right)}\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)-{\rho }_t\left({\mathbf{r}}^{\prime }\right){\omega }_t^{\left(q\right)}\left({\mathbf{r}}^{\prime },\mathbf{r}\right)\right)\right]\hfill \\ \hfill & +{\rho }_t^{\left(q\right)},\hfill \end{array}(4)$

where ${\omega }_t^{(q)}(\mathbf{r},{\mathbf{r}}^{\prime })\equiv {\rho }_t^{(q)}({\mathbf{r}}^{\prime })P{r}_{{\mathbf{r}}^{\prime }\to \mathbf{r}}^{(q)}$ is the flow rate of particles of species $q$ entering the cell at $\mathbf{r}$ from the cell ${\mathbf{r}}^{\prime }$. Alternatively, we could write Equation 4 as

${\rho }_{t+1}^{\left(q\right)}\left(\mathbf{r}\right)={\mathrm{\Phi }}_t^{\left(q\right)}\left(\mathbf{r}\right)-{\mathrm{\Psi }}_t^{\left(q\right)}\left(\mathbf{r}\right)+{\rho }_t^{\left(q\right)},(5)$

where

${\mathrm{\Phi }}_t^{\left(q\right)}\left(\mathbf{r}\right)\equiv {\displaystyle \sum _{⟨{\mathbf{r}}^{\prime }⟩}}\left({\omega }_t^{\left(q\right)}\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)-{\omega }_t^{\left(q\right)}\left({\mathbf{r}}^{\prime },\mathbf{r}\right)\right)(6)$

is the net flow rate of particles of species $q$ into the cell at $\mathbf{r}$ and

${\mathrm{\Psi }}_t^{\left(q\right)}\left(\mathbf{r}\right)\equiv {\displaystyle \sum _{⟨{\mathbf{r}}^{\prime }⟩}}\left({\rho }_t\left(\mathbf{r}\right){\omega }_t^{\left(q\right)}\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)-{\rho }_t\left({\mathbf{r}}^{\prime }\right){\omega }_t^{\left(q\right)}\left({\mathbf{r}}^{\prime },\mathbf{r}\right)\right)(7)$

is the frustrated net flow rate of particles of species $q$ into the cell at $\mathbf{r}$. The frustrated net flow rate measures the particle interaction, i.e., measures the amount of $q$ particles unable to access or leave the cell at $\mathbf{r}$ due to the its occupation and of its neighbouring cells Equations 5-7 provide an alternative representation that offers a more physical perspective on the mean-field dynamics.

2.2 Mean-field regime

When applying Equation 1 into the lattice gas recurrence relation given by Equation 3, we end up with the specific recurrence relation

$\begin{array}{cc}{\rho }_{t+1}^{\left(q\right)}\left(\mathbf{r}\right)\hfill & =\left(1-{\rho }_t\left(\mathbf{r}\right)\right){\displaystyle \sum _{\langle {\mathbf{r}}^{\prime }\rangle }}{\rho }_t^{\left(q\right)}\left({\mathbf{r}}^{\prime }\right)\left(p-\alpha \left(\mathrm{\Delta }\mathbf{r}•{\mathbf{û}}^{\left(q\right)}\left({\mathbf{r}}^{\prime }\right)\right)\right)\hfill \\ \hfill & +{\rho }_t^{\left(q\right)}\left(\mathbf{r}\right)\left[\left(1-4p\right)+{\displaystyle \sum _{\langle {\mathbf{r}}^{\prime }\rangle }}{\rho }_t\left({\mathbf{r}}^{\prime }\right)\left(p+\alpha \left(\mathrm{\Delta }\mathbf{r}•{\mathbf{û}}^{\left(q\right)}\left(\mathbf{r}\right)\right)\right)\right].\hfill \end{array}(8)$

After some algebra (see Supplementary Material A) we obtain the following partial differential equation

$\begin{array}{cc}\frac{\partial {\rho }^{\left(q\right)}}{\partial t}\hfill & =c_1^{\prime }{\nabla }^2{\rho }^{\left(q\right)}-c_2^{\prime }\nabla •\left({\rho }^{\left(q\right)}{\mathbf{û}}^{\left(q\right)}\right)-c_3^{\prime }{\displaystyle \sum _{k\ne q}}\left({\rho }^{\left(k\right)}{\nabla }^2{\rho }^{\left(q\right)}-{\rho }^{\left(q\right)}{\nabla }^2{\rho }^{\left(k\right)}\right)\hfill \\ \hfill & +c_4^{\prime }{\displaystyle \sum _{k=1}^m}\left[{\rho }^{\left(k\right)}\nabla •\left({\rho }^{\left(q\right)}{\mathbf{û}}^{\left(q\right)}\right)+{\rho }^{\left(q\right)}\nabla {\rho }^{\left(k\right)}•{\mathbf{û}}^{\left(q\right)}\right],\hfill \end{array}(9)$

where $c_1^{\prime },c_3^{\prime }\propto \underset{\tau ,ϵ\to 0}{\lim }\frac{p{ϵ}^2}{\tau }$ and $c_2^{\prime },c_4^{\prime }\propto \underset{\tau ,ϵ\to 0}{\lim }\frac{\alpha ϵ}{\tau }$.

In Equation 9, we notice that the constants $c_3^{\prime }$ and $c_4^{\prime }$ are species’ coupling factors and the mathematical path we took from Equation 3 makes them not to depend on each possible pair of species. However, a more general approach would consider these constants as dependent on each possible coupling pair of species, which we will not going to explore in this work.

To study how good a fit is the mean-field approximation is to the particle dynamics, we investigate how the spatial distributions of both methods differ from each other for different sets of the model parameters such as the total number of particles $(n)$, the total number of species $(m)$, and the bias level of movement $(\alpha )$.

3 Discussion

As mean-field approaches approximate the microscopic behavior to an average behavior, their fidelity to the model tends to be limited to a certain range of parameters. This restriction is potentiated when considering cases of highly correlated systems of particles, such as the lattice gas dynamics we study in this work. In particular, we addressed how a system of an arbitrary number of species $m$, each carrying a different static and uniform social field $\mathbf{u}$, influences emerging phenomena. As stated previously, the social field’s direction is the only cause of difference between particles of different species, so that a system of $m>1$, but with ${\mathbf{u}}^{(1)}={\mathbf{u}}^{(2)}=\cdots ={\mathbf{u}}^{(m)}$ is not different from the case $m=1$. In that sense, the complexity of the model comes from different configurations of static floor fields.

Thus, to simplify things a little bit, we use a relatively simple rule to define social fields ${\mathbf{û}}^{(q)}$, where each species has a guaranteed different field from each other in the two-dimensional Cartesian lattice:

${\mathbf{û}}^{\left(q\right)}=\left(\cos \left(\frac{2\pi \left(q-1\right)}{m}\right),\sin \left(\frac{2\pi \left(q-1\right)}{m}\right)\right).(10)$

By defining the static floor fields in Equation 10 as such, we took the first species to have its preferential direction of movement parallel to the x-axis, whilst all other species’ directions are defined in regularly divided angles according to the total number of particles $(m)$.

As initial conditions, we considered each species concentrated in separate regions of the lattice as initial “wave packages” to describe what would be a real physical confrontation scenario, such as found in pedestrian crossings [8, 21]. More specifically, we considered each species’ initial distribution as a bivariate uncorrelated normal distribution with mean values defined in a circular fashion analogous to the static floor fields previously defined with a light difference. To observe confrontation between species, we defined the initial distributions of each as ${\mu }_x^{(q)}=\frac{l}{2}\left[1-\frac{1}{2}\cos \left(\frac{2\pi (q-1)}{m}\right)\right]$ and ${\mu }_y^{(q)}=\frac{l}{2}\left[1-\frac{1}{2}\sin \left(\frac{2\pi (q-1)}{m}\right)\right]$. That particular choice guarantees that agents will have to cross each other at the lattice center, simulating scenarios where species have well-defined positions and confrontation between groups of agents with different objectives is bound to occur as depicted in Figure 2.

Figure 2

Figure 2. Schematic lattice (not to scale) illustrating a system with $m$ species in a situation involving confrontation. Each point marks the mean $\left({\mu }_x^{(q)},{\mu }_y^{(q)}\right)$ of a species’ initial distribution, with arrows indicating the direction of its static floor field. Angular separations are illustrative; actual values depend on $m$. The figure conveys the structure of initial conditions and static fields.

For all species, we define the standard deviation of the initial distributions to be ${\sigma }_x^{(q)}={\sigma }_y^{(q)}=l/32$. The reason behind this specific choice for the standard deviation and not a smaller value is the constraint of the exclusion feature of the lattice model. For instance, a Dirac’s delta function as initial distribution, where ${\sigma }_x^{(q)}={\sigma }_y^{(q)}\to 0$, can only be implemented on MC simulations if we put only one particle in the species mean initial position and we would have to make the system large enough so to have an approximate distribution of choice, which would be impractical computationally. In that sense, to study systems with more than one particle per species, we have to choose initial distributions not so localized.

Specifically, we have examined the cases of $m=1$, $m=2$, and $m=4$, which represent simpler scenarios where the species’ static floor fields align with the $x$- and $y$-axes. This alignment with the Cartesian axes provides a symmetric structure that facilitates a clearer comparison between the two methods via the marginalized distribution relative to the species’ static floor field. For instance, if a species has a preferred direction $\mathbf{û}=\pm {\mathbf{e}}_x$, we compute the marginalized density $\rho (x)$ by summing $\rho (x,y)$ over all $y$-values. This reduction to a one-dimensional distribution simplifies the analysis.

In all cases, we studied a system of $l=128$ with periodic boundary conditions on both directions (toroid). In a previous study, we performed a finite-size scaling analysis of the three-species version of this model using an indirect observable—the bar profits (see Ref. [6]). The results indicated that a lattice size of $l=128$ was sufficiently large to suppress boundary effects. We also used $p=1/4$ in this work, which reflects a scenario of high noise dynamics (see Equation 2) where particles only stay still if the target cell is occupied.

We simulated the particle dynamics of our model via MC algorithm with a shuffled asynchronous update scheme. This means that at any given MC step, we form a list of $n$ particles chosen at random (uniform distributed pseudo-random number generator [25]) from the population of $n$ that compose the system. We then update the position of each sequentially, according to the transition probabilities given by Equations 1, 2 and its neighboring cell states. For the initial conditions, we used a Box-Muller algorithm [25] to generate two independent and normally distributed pseudo-random variables ($X$ and $Y$) to define each particle’s position according to their species. Throughout this work, variables originating from MC simulations were obtained by averaging the time series of each run at each timestep over $n_{\mathit{run}}=10^6/n$ samples generated with same parameter, but different seeds for the pseudo-random number generators.

3.1 $m=1$ case

The simplest scenario we examined involves only a single species. As we defined previously by Equation 10, a single species tends to move according to the static floor field $\mathbf{û}={\mathbf{e}}_x$, with its initial distribution having ${\mu }_x^{(1)}=l/4$ and ${\mu }_y^{(1)}=l/2$ with ${\sigma }_x^{(q)}={\sigma }_y^{(q)}=l/32$. In these simulations, we used $\beta =1$.

In Figure 3, we present plots of the marginal distribution of particles along $x$, where circles represent the Monte Carlo (MC) simulation results and solid lines correspond to predictions from the recurrence relation (mean-field approximation). We examine these distributions for various particle counts, ranging from $n=64$ to $n=448$ in increments of $\mathrm{\Delta }n=64$ and shown at four distinct time steps: (a) $t=50$, (b) $t=250$, (c) $t=450$, and (d) $t=4950$. For this study, we used $\alpha =0.15$, which corresponds to agents with mid range level of impetus as we already shown in [7, 8].

Figure 3

Figure 3. Evolution of normalized marginal distributions of $x$ of MC simulations (circles) and MF (lines) for different time steps (a,b) $t=50$, (b,f) $t=250$, (c,g) $t=450$, and (d,h) $t=4950$ considering one-species $(m=1)$. Plots labeled from (a-d) illustrate the influence of different values of $n$ for a fixed $\alpha =0.15$, whereas plots from (e-h) show the effect of varying $\alpha$ for a fixed $n=192$.

As a general remark, we notice a qualitatively good agreement between the two methods and that the distributions in more densely populated systems shows more positive skewness of the probability density function (PDF) in comparison to systems with fewer particles, revealing the exclusion effects. These exclusion effects are a result of agents leading the pack acting as obstacles for the particles in the bulk. As time passes, we notice that the distributions present increasing dispersion with the “front tail” maintaining synchronicity for different values of $n$. It is important to disclose that in the study that produced Figure 3, the recurrence relation exhibited numerical instability for $n\ge 512$, which the corresponding data is omitted here and throughout the manuscript for conciseness. To understand the source of the instabilities, note that Equation 8 includes the term $1-{\rho }_t(\mathbf{r})$, which, in the context of discrete particle dynamics, only takes values 0 or 1, as previously explained. In this regime, the factor $1-{\rho }_t(\mathbf{r})$ behaves in a stable and well-defined manner. However, under the mean-field approximation, increasing either the total density $N/V$ or the parameter $\alpha$ can lead to locally high average densities. As a result, $1-{\rho }_t(\mathbf{r})$ may become negative, causing ${\rho }_{t+1}(\mathbf{r})$ to also take on unphysical negative values. When this occurs, the normalization of the distribution becomes unstable, as the sum involves both large positive and negative terms—sometimes yielding normalized values with magnitudes as high as $10^{20}$. Thus, providing a quantitative indication that there exists a certain range of the parameters of the model beyond which the approximate method becomes unreliable and serves as justification for the assumption that a tuning parameter, $\beta$, might be necessary to also preserve the mean-field stability. The reader is advised to check the supplementary material to check the form of marginal distributions of particles right before numerical collapse.

Another key parameter we investigated is the bias movement level, $\alpha$, which quantifies the particles’ tendency to move in their species’ preferred direction. For this analysis, we fixed $n=192$, corresponding to a generally low-density scenario as shown, while ensuring sufficient particle interactions to produce the positive skewness observed in the distribution.

In the lower panels of Figure 3 (plots e to h), we present the marginal distributions at time steps (e) $t=50$, (f) $t=250$, (g) $t=450$, and (h) $t=4950$ for $\alpha =0.050,0.100,0.150$, and 0.200. As expected, distributions with higher $\alpha$ values exhibit faster drift velocities in the $+x$ direction. Additionally, increasing $\alpha$ amplifies the skewness, similar to the effect observed with higher $n$. Notably, at longer simulation times, larger $\alpha$ values lead to greater discrepancies between the two methods, as seen in Figure 3h.

Given the phenomena observed in Figure 3, we notice that systems with a greater total number of particles $(n)$ tend to show less agreement, even if qualitatively their shapes follow the same pattern. However, this disparity is not observed when changing the level of movement bias $(\alpha )$ shown in Figure 3. With this in mind, we show, in Figure 4, how the exponent of normalization, $\beta$, influences the shape of the marginal distribution of $x$.

Figure 4

Figure 4. (a) Influence of $\beta$ on the distribution profile for a system with parameters $m=1$, $n=128$, and $\alpha =0.249$. (b) Variation of the average squared deviation $\overline{\mathrm{\Delta }{S}^2}$ as a function of $\beta$, for a single-species system $(m=1)$ with $\alpha =0.15$. (c) Completing the analysis for $m=1$, we show the linear dependence of the optimal value ${\beta }_c$ on the number of agents. The data is fitted with a linear function $ax+b$, yielding $a=-0.000588\pm 5.3\times 10^{-5}$ and $b=1.17\pm 0.017$. (d) By extending the same calculations to systems with multiple species ($m=2$ and $m=4$), we observe that the results are well described by power function fits for ${\beta }_c$, using the form $a{x}^{-b}$. The fitted parameters are: $m=2$: $a=3.3\pm 0.11$, $b=0.193\pm 0.0065$; $m=4$: $a=3.04\pm 0.23$, $b=0.176\pm 0.014$.

Figure 4a illustrates the influence of $\beta$ on the distribution profile for a system with parameters $m=1$, $n=128$, and $\alpha =0.249$. As discussed in the previous section, the normalization constraint in the recurrence relation can be interpreted through an ansatz in which the number of particles is raised to a power $\beta$, with $\beta$ expected to be close to 1, if not exactly 1. As such, this parameter serves to tune the agreement between methods at the same time, it allows for systems with number of particles expected to show numerical collapse to evolve in a well-behaved manner.

To better measure the agreement between the two methods and the range of parameters where the mean-field shows numerical stability, we compared both methods using the time average of the squared difference of the spatial entropies of each method, i.e.,

$\overline{{\mathrm{\Delta }S}^2}\equiv \frac{1}{t_{\max }}{\displaystyle \sum _{t=1}^{t_{\max }}}{\left(S_{t_{\text{MC}}}-S_{t_{\text{MF}}}\right)}^2,(11)$

where $S_t={\sum }_{\mathbf{r}}{\rho }_t(\mathbf{r})\ln {\rho }_t(\mathbf{r})$. The loss function introduced in Equation 11 is one among several possible choices for quantifying the discrepancy between the MC and MF results. However, as we demonstrate below, it proved to be particularly suitable for optimization purposes. For comparison, we also tested an alternative loss function based on the squared differences of local densities, as illustrated in Supplementary Figure S2 and defined in its caption in the Supplementary Material. This alternative metric failed to produce a well-behaved optimization landscape, which made it ineffective for identifying optimal values of $\beta$. For this reason, we opted not to pursue it further.

In Figure 4b, we show $\overline{{\mathrm{\Delta }S}^2}$ as a function of the exponent $\beta$ for different values of $n$ and for a fixed simulation time of $t_{\max }=10^4$ time steps. We studied the range of $\beta \in [0.5:1.5]$. In this range, we observed numerical instability of the recurrence relation in all curves shown for values of $\beta$ greater than certain values, which are reflected by the discontinuity of the curves. More specifically, greater $n$ shows lower $\beta$’s where the numerical collapse occurs. However, within the range of numerical stability, we observe that, for $n\le 256$, methods show an optimal agreement at $\beta \equiv {\beta }_c>1$, shown by the local minima which become closer to 1 as $n$ increases. For $n\ge 256$, we notice that the numerical instability of the MF makes the lowest stable $\beta$ assume the stance of optimal value, i.e., it becomes the value that minimizes $\overline{{\mathrm{\Delta }S}^2}$.

Figure 4c completes the analysis for $m=1$, by showing the linear dependence of the optimal value ${\beta }_c$ on the number of agents. The data is fitted with a linear function $ax+b$, yielding $a=-0.000588\pm 5.3\times 10^{-5}$ and $b=1.17\pm 0.017$. Such behavior will be different of upper values of $m$ as it is previously presented in Figure 4d. However, we must first study the particle distribution for these cases under varying parameters, which will be performed in the next sections.

3.2 $m=2$ case

The case $m=2$ is the simplest case within our proposed framework of initial conditions (see Figure 2) where two species of particles move in counterflow, which is a scenario vastly studied in the literature, once it can describe a two species systems with opposing responses to an external field such as oppositely charged coloids [9, 20] (external electrical field) or chromatographic columns [10, 26–29] (gravitational field), to cite a few. For this case, initial distributions assume ${\mu }_x^{(1)}=l/4$ and ${\mu }_y^{(1)}=l/2$, ${\mu }_x^{(2)}=3l/4$ and ${\mu }_y^{(2)}=l/2$ as we defined in the previous section. Static floor fields are then ${\mathbf{û}}^{(1)}={\mathbf{e}}_x$ and ${\mathbf{û}}^{(2)}=-{\mathbf{e}}_x$ as stated by Equation 10, which still carries a rather trivial symmetry in the $y$-direction (diffusive profile) making observations of the marginal distribution on the $x$ coordinate rather more interesting because is where the confrontation occurs.

In Figure 5, we show the marginal distributions of $x$ of the MC simulations with circles and lines for the MF approach. We studied two values of the total number of particles shown simultaneously: $n=128$ shown in purple and $n=256$ shown in green for the time steps (a) $t=50$, (b) $t=250$, (c) $t=450$, and (d) $t=4950$, where both species appear as different peaks with same color for each case. In this first case, we used $\alpha =0.15$ and $\beta =1.0$, and we observe that, similarly to the $m=1$ case, a greater total number of particles makes the distributions more asymmetric due to the exclusion interaction of particles of the same species. We do not observe, however, a clear influence of the opposing species on the shape of the distribution. We can interpret this “lack” of interaction deformity occurring because of the exclusion effects of the nearest neighbor interaction, suggesting that a sort of shield is formed by the same-species particles in a low-density regime. This warrants further investigation in future work.

Figure 5

Figure 5. Evolution of normalized marginal distributions of $x$ of the MC simulations (circles) and MF (lines) for different time steps (a,e) $t=50$, (b,f) $t=250$, (c,g) $t=450$, and (d,h) $t=4950$. Plots (a-d) show the influence of different values of $n$ for $\alpha =0.15$, whilst plots (e-h) show different $\alpha$’s for a fixed $n=192$.

While both methods show reasonably good agreement in Figure 5, we find that the agreement can be further improved by optimizing $\beta$. Revisiting Figure 4d and extending our earlier analysis for $m=1$ to $m=2$, we observe that the critical value ${\beta }_c$ minimizing $\overline{{\mathrm{\Delta }S}^2}$ follows a power-law dependence, $a{x}^{-b}$. For $m=2$, the fit yields $a=3.3\pm 0.11$ and $b=0.193\pm 0.0065$, contrasting with the linear behavior found for $m=1$. This indicates that increasing the model complexity (from $m=1$ to $m=2$) modifies the parameter dependence. However, no significant changes are observed beyond $m=2$ up to $m=4$, as seen in Figure 4D—a point we will address in detail later. Importantly, we note that lower-density cases do not provide the best agreement between the two methods, revealing non-trivial model-specific characteristics that lack complete theoretical understanding.

3.3 $m=4$ case

The last case we discuss in this work is the case with four species moving along the $x-$ and $y-$directions. More precisely, their initial distributions have mean value: ${\mu }_x^{(1)}=l/4$ and ${\mu }_y^{(1)}=l/2$, ${\mu }_x^{(2)}=3l/4$ and ${\mu }_y^{(2)}=l/2$, ${\mu }_x^{(3)}=l/2$ and ${\mu }_y^{(3)}=l/4$, and ${\mu }_x^{(4)}=l/2$ and ${\mu }_y^{(4)}=3l/2$ as we expected to be by our previous definition. Similarly, the static floor fields assume the form ${\mathbf{û}}^{(1)}={\mathbf{e}}_x$, ${\mathbf{û}}^{(2)}={\mathbf{e}}_y$, ${\mathbf{û}}^{(3)}=-{\mathbf{e}}_x$, and ${\mathbf{û}}^{(4)}=-{\mathbf{e}}_y$ (see Equation 10). Similarly to the two prior cases, having an initial condition with radial symmetry about the lattice center, also makes the four species present a symmetric “cross section” of the distribution, i.e., the distribution on the direction perpendicular to ${\mathbf{û}}^{(q)}$ shows no asymmetry. Thus, studying the marginal distribution of $x$ of the species $q=1$ or $q=3$ would be no different from studying the marginal distribution of $y$ of the species $q=2$ and $q=4$ given the system’s initial arrangement and static floor fields (see Figure 2).

With this in mind, Figure 6 displays the particle distribution ${\rho }^{(1)}(x)$ for different values of $n$ using $\alpha =0.249$ at four time points: (a) $t=50$, (b) $t=250$, (c) $t=450$, and (d) $t=4950$. Our analysis reveals two key observations:

Figure 6

Figure 6. Evolution of normalized marginal distributions of $x$ of the MC simulations (circles) and MF (lines) for different time steps (a,e) $t=50$, (b,f) $t=250$, (c,g) $t=450$, and (d,h) $t=4950$. Plots (a-d) show the influence of different values of $n$ for $\alpha =0.249$, whilst plots (e-h) show different $\alpha$’s for a fixed $n=256$.

(1) First, the discrepancy between methods becomes more pronounced for larger $n$ values in this regime of highly driven agents. Despite this disagreement, the MC simulations demonstrate the formation of transient condensates for $n=320$, as evidenced by the emerging peak at $t=250$.

(2) Second, the mean-field solution exhibits anomalous behavior at $t=250$, particularly in the left tail of ${\rho }^{(1)}(x)$, which appears bent near $x\approx 100$. Notably, however, both methods show qualitatively good agreement in the front tail of ${\rho }^{(1)}(x)$, suggesting synchronized drifting profiles in the steady state.

Building on our analysis of the $m=4$ case, Figure 6 presents the evolution of ${\rho }^{(1)}(x)$ distributions for various $\alpha$ values in a system with $n=256$ particles. We also examine the same four time points: (e) $t=50$, (f) $t=250$, (g) $t=450$, and (h) $t=4950$. The results reveal that while the wave fronts remain synchronized between MC simulations and mean-field theory, the MC distributions exhibit significantly greater dispersion. Returning to Figure 4d, we apply the same minimization procedure for $\overline{{\mathrm{\Delta }S}^2}$ used for the m = 1 and m = 2 cases. For m = 4, we obtain an excellent power-law fit identical in form to the m = 2 case, with fitted parameters a = 3.04 $\pm$ 0.23 and b = 0.176 $\pm$ 0.014.

It is important to note that, despite some differences in the evolution of particle profiles, a qualitatively similar behavior is observed in both cases when the systems are initialized with a normal distribution of particles in the environment. The interactions among particles—especially the interaction with counterflowing particles—tend to deform the initial distribution shape. Over time, however, after several rounds of confrontation and interaction, the particles begin to reorganize themselves, ultimately restoring a Gaussian-like behavior.

This transition from one Gaussian state to another is illustrated in Figure 7. In Figures 7a–e, we show the results from MC simulations, while Figure 7f corresponds to the mean field approximation. All plots refer to the case of $m=4$ species, $n=256$ particles, and interaction parameter $\alpha =0.249$.

Figure 7

Figure 7. Plots (a–e) show the particle distribution at different times obtained via Monte Carlo simulations. At $t=50$, significant deviations from a normal distribution are observed. As time progresses, the system gradually returns to a Gaussian-like profile, culminating in a noisy steady-state distribution at $t=400$, as seen in plot (e). Plot (f) summarizes the corresponding results for the mean-field approximation. Here, the quality of the Gaussian fit over time is quantified by the coefficient of determination $R^2$, represented by dark yellow circles. Initially, interactions among particles cause strong deviations, leading to a transient regime. At later times, the system converges again to a Gaussian distribution—less noisy than that of the Monte Carlo steady state. The inset in (a) compares normal and log-normal fits on a semi-logarithmic scale. The insets in (f) illustrate the particle profile at different stages of the evolution.

In Figure 7a, we observe the particle distribution at $t=50$. We observe that the initial Gaussian behavior is disrupted due to particle interactions. The red and blue curves represent normal and log-normal fits, respectively, with the latter providing a more accurate fit at this stage. The inset plot shows the same trend on a semi-log scale. By $t=150$, the distribution becomes strongly non-Gaussian, with a pronounced peak emerging, as shown in Figure 7b.

At $t=250$, the shape of the distribution becomes nearly triangular (Figure 7c), signaling the beginning of a return to Gaussianity. A good Gaussian fit is eventually recovered in Figure 7d, and a noisy but stable Gaussian distribution is observed at $t=400$, as shown in Figure 7e. This reflects the establishment of a steady state after extensive particles interactions.

Finally, in Figure 7f, we analyze the mean-field approximation. Unlike the MC simulations, which display the full evolution of the particle profile, we summarize the results using the coefficient of determination $R^2$, which quantifies the quality of a normal distribution fit over time. These values are represented by circles. Initially, we observe a sharp decline in $R^2$, indicating the breakdown of Gaussianity (illustrated in the inset). As $t$ increases, the system goes through a transient (points inside blue rectangle) phase before recovering a Gaussian profile, again shown in the inset. This confirms that the same Gaussian-to-Gaussian behavior occurs in the mean-field scenario, in agreement with the MC simulations presented in Figures 7a–e.

4 Conclusion

In this work, we proposed a mean-field approximation to describe the dynamics of an agent-based model of $m$-species of particles. Specifically, we proposed a recurrence relation that shows the evolution of lattice gas dynamics with an asynchronous updating scheme in a Von Neumann neighborhood.

We showed that the MF method shows a good agreement with the particle dynamics implemented through MC simulation using initial conditions and static floor fields carrying polar symmetry in relation to the center of the lattice. We studied the cases of one, two, and four species, which, with the initial conditions, made for a sort of “controlled” environment so that we could access more easily the agreement between methods.

The MF approach carried a normalization constant exponent $\beta$, which relied on the parameters $n$ and $m$, because the agreement of methods is strongly influenced by the complexity of the interactions. We showed that when one species is considered, its optimal value, ${\beta }_c$, has a linear dependence on the total number of particles, when $m=2$ and $m=4$ it shows to be a function of some power of $n$.

We highlighted that our proposed mean-field approach is not bulletproof by showing what specific set of parameters made it show numerical instability due to competing terms inside the recurrence relation.

Despite studying only three cases, our approach is general enough to allow an arbitrary number of particles, even for initial conditions and/or species preferential directions of motion different from those studied in this work. Not only that, our framework proposed by Equation 3 presents a general feature of application to asynchronous lattice gas models, even reproducing the recurrence relation of other models such as the deterministic two-species dynamics studied in [21].

We are currently exploring higher systems with more than four species, which means that some species are going to present static floor field at angles lower than $90°$ with the $x-$ and $y-$direction. As a consequence, we expect to observe ${\beta }_c$ presenting a more complex dependency with $n$.

Finally, our model reveals a distinctive pattern we term the “Gaussian-to-Gaussian” transition. In this phenomenon, particles initially distributed in a Gaussian configuration undergo significant deformation of this pattern, only to eventually recover a Gaussian distribution after multiple interaction cycles. Both Monte Carlo simulations and mean-field theory reproduce this effect, though they follow different temporal evolution pathways.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

ES: Data curation, Writing – original draft, Supervision, Methodology, Investigation, Conceptualization, Software, Funding acquisition, Visualization, Writing – review and editing, Resources, Project administration, Validation, Formal Analysis. RD: Methodology, Formal Analysis, Conceptualization, Project administration, Validation, Visualization, Supervision, Writing – original draft, Funding acquisition, Software, Investigation, Resources, Data curation, Writing – review and editing. SG: Project administration, Formal Analysis, Data curation, Validation, Writing – review and editing, Writing – original draft, Resources, Conceptualization, Supervision, Visualization, Funding acquisition, Software, Investigation, Methodology.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. We gratefully acknowledge the partial support provided by CNPq: E. V. ES under grant 153315/2024-5, RS under grant 304575/2022-4, and SG under grant 314738/2021-5.

Acknowledgments

We appreciate the availability of computational resources from the Lovelace cluster at IF-UFRGS.

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.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2025.1648895/full#supplementary-material

References

1. Guo X, Angulo A, Tavakoli A, Robartes E, Chen TD, Heydarian A. Rethinking infrastructure design: evaluating pedestrians and vrus’ psychophysiological and behavioral responses to different roadway designs. Scientific Rep (2023) 13:4278. doi:10.1038/s41598-023-31041-9

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Geoerg P, Schumann J, Boltes M, Kinateder M. How people with disabilities influence crowd dynamics of pedestrian movement through bottlenecks. Scientific Rep (2022) 12:14273. doi:10.1038/s41598-022-18142-7

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Nicholls VI, Wiener J, Meso AI, Miellet S. The impact of perceptual complexity on road crossing decisions in younger and older adults. Scientific Rep (2024) 14:479. doi:10.1038/s41598-023-49456-9

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Zhang H, Xu J, Jia L, Shi Y. Modelling the walking behavior of pedestrians in the junction with chamfer zone of subway station. Physica A: Stat Mech its Appl (2022) 602:127656. doi:10.1016/j.physa.2022.127656

CrossRef Full Text | Google Scholar

5. Huang D, Yang Y, Peng X, Huang J, Mo P, Liu Z, et al. Modelling the pedestrian’s willingness to walk on the subway platform: a novel approach to analyze in-vehicle crowd congestion. Transportation Res E: Logistics Transportation Rev (2024) 181:103359. doi:10.1016/j.tre.2023.103359

CrossRef Full Text | Google Scholar

6. Stock EV, da Silva R. Lattice gas model to describe a nightclub dynamics. Chaos, Solitons and Fractals (2023) 168:113117. doi:10.1016/j.chaos.2023.113117

CrossRef Full Text | Google Scholar

7. Stock EV, da Silva R, Gonçalves S. Nightclub bar dynamics: statistics of serving times. The Eur Phys J B (2024) 97:179. doi:10.1140/epjb/s10051-024-00803-3

CrossRef Full Text | Google Scholar

8. Stock EV, da Silva R. Exploring crossing times and congestion patterns at scramble intersections in pedestrian dynamics models: a statistical analysis. Physica A: Stat Mech its Appl (2024) 649:129942. doi:10.1016/j.physa.2024.129942

CrossRef Full Text | Google Scholar

9. Vissers T, van Blaaderen A, Imhof A. Band formation in mixtures of oppositely charged colloids driven by an ac electric field. Phys Rev Lett (2011) 106:228303. doi:10.1103/PhysRevLett.106.228303

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Stock EV, da Silva R, Fernandes HA. Statistics, distillation, and ordering emergence in a two-dimensional stochastic model of particles in counterflowing streams. Phys Rev E (2017) 96:012155. doi:10.1103/PhysRevE.96.012155

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Zhang J, Seyfried A. Comparison of intersecting pedestrian flows based on experiments. Physica A: Stat Mech its Appl (2014) 405:316–25. doi:10.1016/j.physa.2014.03.004

CrossRef Full Text | Google Scholar

12. Helbing D, Farkas I, Vicsek T. Simulating dynamical features of escape panic. Nature (2000) 407:487–90. doi:10.1038/35035023

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Majumdar SN, Evans MR, Zia RKP. Nature of the condensate in mass transport models. Phys Rev Lett (2005) 94:180601. doi:10.1103/physrevlett.94.180601

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Katz S, Lebowitz JL, Spohn H. Phase transitions in stationary nonequilibrium states of model lattice systems. Phys Rev B (1983) 28:1655–8. doi:10.1103/PhysRevB.28.1655

CrossRef Full Text | Google Scholar

15. Helbing D, Molnár P. Social force model for pedestrian dynamics. Phys Rev E (1995) 51:4282–6. doi:10.1103/PhysRevE.51.4282

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Oliveira CLN, Vieira AP, Helbing D, Andrade JS, Herrmann HJ. Keep-left behavior induced by asymmetrically profiled walls. Phys Rev X (2016) 6:011003. doi:10.1103/PhysRevX.6.011003

CrossRef Full Text | Google Scholar

17. Kleintjens L, Van der Haegen R, van Opstal L, Koningsveld R. Mean-field lattice-gas modelling of supercritical phase behavior. The J Supercrit Fluids (1988) 1:23–30. doi:10.1016/0896-8446(88)90006-X

CrossRef Full Text | Google Scholar

18. Bernard MO, Plapp M, Jfmc G. Mean-field kinetic lattice gas model of electrochemical cells. Phys Rev E (2003) 68:011604. doi:10.1103/PhysRevE.68.011604

CrossRef Full Text | Google Scholar

19. Dickman R. Mean-field theory of the driven diffusive lattice gas. Phys Rev A (1988) 38:2588–93. doi:10.1103/physreva.38.2588

PubMed Abstract | CrossRef Full Text | Google Scholar

20. da Silva R, Hentz A, Alves A. Stochastic model of self-driven two-species objects inspired by particular aspects of a pedestrian dynamics. Physica A: Stat Mech its Appl (2015) 437:139–48. doi:10.1016/j.physa.2015.05.104

CrossRef Full Text | Google Scholar

21. Cividini J, Hilhorst HJ, Appert-Rolland C. Crossing pedestrian traffic flows, the diagonal stripe pattern, and the chevron effect. J Phys A: Math Theor (2013) 46:345002. doi:10.1088/1751-8113/46/34/345002

CrossRef Full Text | Google Scholar

22. Appert-Rolland C, Cividini J, Hilhorst H, Degond P. Pedestrian flows: from individuals to crowds. Transportation Res Proced (2014) 2:468–76. doi:10.1016/j.trpro.2014.09.062

CrossRef Full Text | Google Scholar

23. Stock EV, da Silva R. Coexistence and crossover phenomena in a fermi-like model of particles in counterflowing streams. Phys Rev E (2020) 102:022139. doi:10.1103/PhysRevE.102.022139

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Schadschneider A, Chraibi M, Seyfried A, Tordeux A, Zhang J. Pedestrian dynamics: from empirical results to modeling. Cham: Springer International Publishing (2018). p. 63–102. doi:10.1007/978-3-030-05129-7_4

CrossRef Full Text | Google Scholar

25. Press WH, Sa T, Vetterling WT, Flannery BP. Numerical recipes in fortran 77: the art of scientific computing. vol. 1 (1996).

Google Scholar

26. Giddings JC, Byring H. A molecular dynamic theory of chromatography. The J Phys Chem (1955) 59:416–21. doi:10.1021/j150527a009

CrossRef Full Text | Google Scholar

27. da Silva R, Lamb LC, Lima EC, Dupont J. A simple combinatorial method to describe particle retention time in random media with applications in chromatography. Physica A: Stat Mech its Appl (2012) 391:1–7. doi:10.1016/j.physa.2011.08.006

CrossRef Full Text | Google Scholar

28. Stock EV, da Silva R, da Cunha CR. Numerical study of condensation in a fermi-like model of counterflowing particles via gini coefficient. J Stat Mech Theor Exp (2019) 2019:083208. doi:10.1088/1742-5468/ab333d

CrossRef Full Text | Google Scholar

29. da Silva R, Stock EV. Mobile-to-clogging transition in a fermi-like model of counterflowing particles. Phys Rev E (2019) 99:042148. doi:10.1103/PhysRevE.99.042148

PubMed Abstract | CrossRef Full Text | Google Scholar