We develop a stochastic model that describes the dynamics of multiple groups of agents interacting on a lattice of any dimension.
Importantly, our stochastic model combines two mechanisms:
Group-specific floor-field functions specify the global motion of the various groups, and agents belonging to different groups interact.
Different groups interact with one another via slowdown interactions, similar to those considered in~\cite{ckpt14} on a one-dimensional lattice.

For the sake of clarity and brevity, we consider two groups of agents interacting on a lattice of dimension two. 
This allows us to focus on stochastic model construction, as well as the mesoscopic and macroscopic limits, without burdening the exposition with the lengthy formulas that would result from imposing more complex interaction rules between at least three groups.

Our stochastic model is a continuous-time Markov chain on a lattice.
In particular, agents from each of the two groups move on the $N_1\times N_2$ lattice $\mathcal{L}$.
We represent agents in the first and second groups by variables $\sigma^A_{j,k}(t) \equiv \sigma^A(j,k,t)$ and  $\sigma^B_{j,k}(t) \equiv \sigma^B(j,k,t)$, respectively, 
where $(j,k)$ is the position on the lattice and $t$ denotes time:
\begin{align*}
	\sigma^A_{j,k}(t) &= \left\{\begin{array}{ll} 1,& \text{if at time }t\text{ cell }(j,k)\text{ is occupied by an agent from group A;}\\
		0, & \text{otherwise,}\end{array}\right.\\
	\sigma^B_{j,k}(t) &= \left\{\begin{array}{ll} 1,& \text{if at time }t\text{ cell }(j,k)\text{ is occupied by an agent from group B;}\\
		0, & \text{otherwise.}\end{array}\right.
\end{align*}
Let $\sigma^{A} (t)$ and $\sigma^{B} (t)$ denote the collections of occupation values $\sigma_{j,k}^{A} (t)$ and $\sigma_{j,k}^{B} (t)$, respectively, over lattice sites $(j,k) \in \mathcal{L}$.
To simplify the presentation, we notationally suppress the dependence of $\sigma^A$ and $\sigma^B$ on time in what follows, unless we wish to explicitly emphasize time.

When developing an agent-based model for directional particle movement, the key challenge is the mathematical formulation of agent velocities and agent-agent interaction rules.
Here, we adopt the following general principles.
\begin{enumerate}[leftmargin=*, labelindent=\parindent, itemsep=0ex, label=\textbf{(P\arabic*)}, ref=P\arabic*]
\item
\label{li:floor-fields}
(\textbf{Floor fields})
In isolation, agent transition probabilities are determined by a velocity field that depends on the group identity of the agent.
\item
\label{li:exclusion}
(\textbf{Exclusion principle})
Two agents from the same group cannot simultaneously occupy a single lattice site.
\item
\label{li:slowdown}
(\textbf{Slowdown interactions})
Agent movement depends on the local configuration of agents from the other group.
In particular:
\begin{enumerate}[topsep=0ex, label=\textbf{(\alph*)}]
\item
Agents from different groups can occupy the same lattice site at the same time.
\item
Agent transition probabilities decrease when at least one agent from the other group occupies the same cell, or a neighboring cell in the direction of motion.
The amount of decrease depends on the local configuration.
\end{enumerate}
\item
\label{li:no-diagonal}
(\textbf{Lattice motion})
Agent motion is horizontal or vertical.
Consequently, diagonal motion can only result from two transitions (vertical followed by horizontal or vice-versa).
\end{enumerate}

With these general principles in place, we now precisely formulate our stochastic model.
Since the model is a continuous-time Markov chain, we must describe the admissible transitions and the corresponding transition probabilities.
We do so as follows.
\begin{itemize}
\item
  \textit{Admissible transitions.}
  We assign a floor velocity field (defined on the lattice $\mathcal{L}$) to each group of agents.
  Importantly, different groups may have different floor fields.
  For a given agent located at a given lattice site, the horizontal and vertical projections of the corresponding floor field vector specify the only two possible transitions for this agent (in accord with (\ref{li:no-diagonal})).
  (Only one transition is possible if the floor field vector is itself horizontal or vertical.)
  For our stochastic model, the admissible transitions are those specified by the floor fields, subject to the additional constraints imposed by the exclusion principle (\ref{li:exclusion}).

\item
  \textit{Transition probabilities.}
  For each group of agents, we define an \textit{interaction velocity field} that specifies transition rates for the stochastic model.
  These interaction velocity fields are obtained by modifying the floor fields (depending on local system configuration) to account for slowdown interactions, in accord with (\ref{li:slowdown}).
  For each admissible horizontal or vertical transition, the length of the corresponding projection of the associated interaction velocity vector gives the transition rate.
  Multiplying the transition rates by short time intervals $\Delta t$ produces the transition probabilities.
\end{itemize}


\subsection*{Floor velocity fields.}

Our two-dimensional stochastic model is significantly more flexible than the one-dimensional model in~\cite{ckpt14}, in terms of domain geometry, number of interacting groups, and complexity of motion.
When constructing stochastic models in dimension one, specifying direction of movement is straightforward (\textit{e.g.}~\cite{ckpt14,burger2}).
Such one-dimensional models can produce PDEs for bidirectional wave propagation (\textit{e.g.}~\cite{goatin,ardm11}).
In dimension at least two, however, there exists a rich set of possibilities for agent motion.
A flexible modeling framework should allow direction of motion to vary in space, in order to handle scenarios such as obstacle avoidance and motion in complex geometries.
To achieve such flexibility, we use floor velocity fields (see \textit{e.g.}~\cite{bksz01,kirch02}) to specify direction of motion for each group in isolation (as indicated in~(\ref{li:floor-fields}).
Importantly, different groups may obey different floor velocity fields.

For the sake of clarity, we consider the situation when each group is moving toward 
its own \emph{target exit point}.
We denote these exit points as $(j_0^{A}, k_0^{A})$ and $(j_0^{B}, k_0^{B})$ for groups A and B, respectively.
Interesting dynamics occur when the target points differ, a case we examine in detail.
We assume that target exit points do not move over time, although one can easily extend our model to include moving target exit points.
Let $\phi^A(j,k, j_0^A, k_0^A)$ and $\phi^B(j,k, j_0^B, k_0^B)$ denote the floor fields defined over $(j,k) \in \mathcal{L}$.
(We often omit the explicit dependence of these fields on the target exit points.)
We assume that $\phi^{A} = (\phi_{1}^{A}, \phi_{2}^{A})$ and $\phi^{B} = (\phi_{1}^{B}, \phi_{2}^{B})$ arise from potential functions $\psi^{A}$ and $\psi^{B}$:
\begin{equation}
\label{eq:phi}
\phi(j,k) = \left\{\begin{array}{rl}-\frac{\nabla \psi(j,k)}{||\nabla\psi(j,k)||_{1}},&||\nabla\psi(j,k)||_{1} \neq 0,\\ \bm{0},& \text{otherwise.}\end{array}\right.,
\end{equation}
We focus in this work on the quadratic potential
\begin{equation}
\label{eq:quadpsi}
\psi (j,k) = (j-j_{0})^{2} + (k-k_{0})^{2}.
\end{equation}

Our theoretical framework does not require that the floor fields be gradient fields, nor does it require target exit points.
In fact, our framework places no restrictions on the structure of the floor fields.
We simply use the quadratic potential~\eqref{eq:quadpsi} for simulations.

\subsection*{Interaction velocity fields.}

With the floor fields in place, we define a velocity field for each group that captures interaction effects.
Used to determine the probabilities that an agent moves to neighboring cells, these \textit{interaction velocity fields} account for the presence or absence of agents from the other group in nearby cells (slowdown interactions). 
In particular, the interaction velocity field $V^A(j,k,\sigma^{B}) = (V^A_1(j,k,\sigma^{B}), V^A_2(j,k,\sigma^{B}))$ for agents in 
group $A$ depends on the local configuration of group B agents.
This velocity field is defined as follows:

\begin{align}
	V^A_1(j,k,\sigma^{B}) &= \begin{cases}
			c_0\phi^A_1(j,k), & \text{if } \sigma^B_{j-1,k} = \sigma^B_{j,k} = \sigma^B_{j+1,k} = 0,\\
			\phi^A_1(j,k)[c_1H(-\phi^A_1(j,k)) + c_0H(\phi^A_1(j,k))], & \text{if } \sigma^B_{j-1,k} = 1, \spp \sigma^B_{j,k} = \sigma^B_{j+1,k} = 0,\\
			\phi^A_1(j,k)[c_0H(-\phi^A_1(j,k)) + c_1H(\phi^A_1(j,k))], & \text{if } \sigma^B_{j-1,k} = \sigma^B_{j,k} = 0, \spp \sigma^B_{j+1,k} = 1,\\
			c_1\phi^A_1(j,k), & \text{if } \sigma^B_{j-1,k} = \sigma^B_{j+1,k} = 1, \spp \sigma^B_{j,k} = 0,\\
			\phi^A_1(j,k)[c_2H(-\phi^A_1(j,k)) + c_3H(\phi^A_1(j,k))], & \text{if } \sigma^B_{j-1,k} = 0, \spp \sigma^B_{j,k} = \sigma^B_{j+1,k} = 1,\\
			\phi^A_1(j,k)[c_3H(-\phi^A_1(j,k)) + c_2H(\phi^A_1(j,k))], & \text{if } \sigma^B_{j-1,k} = \sigma^B_{j,k} = 1, \spp \sigma^B_{j+1,k} = 0,\\
			c_2\phi^A_1(j,k), & \text{if } \sigma^B_{j-1,k} = \sigma^B_{j+1,k} = 0, \spp \sigma^B_{j,k} = 1,\\
			c_3\phi^A_1(j,k), & \text{if }   \sigma^B_{j-1,k} = \sigma^B_{j+1,k} = \sigma^B_{j,k} = 1,
\end{cases}\label{vel1}\\
	V^A_2(j,k,\sigma^{B}) &= \begin{cases}
                   c_0\phi^A_2(j,k), & \text{if } \sigma^B_{j,k-1} = \sigma^B_{j,k} = \sigma^B_{j,k+1} = 0,\\
			\phi^A_2(j,k)[c_1H(-\phi^A_2(j,k)) + c_0H(\phi^A_2(j,k))], & \text{if } \sigma^B_{j,k-1} = 1, \spp \sigma^B_{j,k} = \sigma^B_{j,k+1} = 0,\\
			\phi^A_2(j,k)[c_0H(-\phi^A_2(j,k)) + c_1H(\phi^A_2(j,k))], & \text{if } \sigma^B_{j,k-1} = \sigma^B_{j,k} = 0, \spp \sigma^B_{j,k+1} = 1,\\
			c_1\phi^A_2(j,k), & \text{if } \sigma^B_{j,k-1} = \sigma^B_{j,k+1} = 1, \spp \sigma^B_{j,k} = 0,\\
			\phi^A_2(j,k)[c_2H(-\phi^A_2(j,k)) + c_3H(\phi^A_2(j,k))], & \text{if } \sigma^B_{j,k-1} = 0, \spp \sigma^B_{j,k} = \sigma^B_{j,k+1} = 1,\\
			\phi^A_2(j,k)[c_3H(-\phi^A_2(j,k)) + c_2H(\phi^A_2(j,k))], & \text{if } \sigma^B_{j,k-1} = \sigma^B_{j,k} = 1, \spp \sigma^B_{j,k+1} = 0,\\
			c_2\phi^A_2(j,k), & \text{if } \sigma^B_{j,k-1} = \sigma^B_{j,k+1} = 0, \spp \sigma^B_{j,k} = 1,\\
			c_3\phi^A_2(j,k), & \text{if } \sigma^B_{j,k-1} = \sigma^B_{j,k+1} = \sigma^B_{j,k} = 1.
\end{cases}\label{vel2}
\end{align}

Expressions~\eqref{vel1} and \eqref{vel2} are rather involved, because we must explicitly 
consider all possible configurations for the second group, $B$, in the vicinity
of the lattice site $(j,k)$.
However, when only one group of pedestrians is present
(\textit{e.g.} $\sigma_{j,k}^{B} (t) = 0$ for all $(j,k)$ and all $t$), the interaction velocity field assumes the simple form
\begin{align*}
& V^A_1(j,k) = c_0\phi^A_1(j,k) = c_0 \frac{j_0 - j}{|j-j_{0}| + |k-k_{0}|}, \\
& V^A_2(j,k) = c_0\phi^A_2(j,k) = c_0 \frac{k_0 - k}{|j-j_{0}| + |k-k_{0}|},
\end{align*}
where we have used the quadratic potential~\eqref{eq:quadpsi} as a particular example.
The horizontal component of this velocity, $V^A_1(j,k)$, is positive if $j < j_0$ and negative if $j > j_0$. 
This means that the particle moves toward the target coordinate, $j_0$.
The vertical component $V^A_2(j,k)$ behaves analogously.

Returning to the general setup, we define the interaction velocity field for group B, $V^{B}$, by invoking the substitutions $A \to B$ and $B \to A$ throughout~\eqref{vel1} and \eqref{vel2}.

The functions $\phi_{1}^{A}$ and $\phi_{2}^{A}$ in~\eqref{vel1} and \eqref{vel2} give the horizontal and vertical components of the floor field, respectively, and $H(x)$ is the Heaviside function with $H(0) = 0$. 
The Heaviside function is used to determine the direction of movement for the agent in cell $(j,k)$ (\textit{i.e.} the direction in which the velocity given by $\phi_1(j,k)$ or $\phi_2(j,k)$ is positive).
Further, it is used to ensure that a slowdown occurs if and only if an agent from the other group occupies the same cell as the agent in cell $(j,k)$, or the adjacent cell in the desired direction of movement, or both. Our formulation ensures that agents from the other group positioned \emph{behind} the agent in cell $(j,k)$ do not contribute to slowdown.

The velocity scalings $c_{0}$, $c_{1}$, $c_{2}$, and $c_{3}$ in $V^{A}$ quantify slowdown linked to the local configuration around a given agent.
For instance, consider the horizontal motion of an agent from group A in cell $(j,k)$; this is described by the interaction velocity component $V^A_1(j,k,\sigma^{B})$.
If there are no agents from group B in $(j,k)$ and horizontally adjacent
cells $(j-1,k)$ and $(j+1,k)$, then no slowdown occurs in the horizontal direction - the velocity scaling is $c_{0}$. 
However, if an agent from group B also occupies cell $(j,k)$, but no agents from group B appear in horizontally adjacent cells, then a slowdown should occur because the group A agent in cell $(j,k)$ must interact with the group B agent at the same location.
In this case, the velocity scaling is $c_{2}$.
Slowdown also results from the presence of group B agents in the direction of motion.
For instance, suppose that the floor velocity component $\phi_1^A(j,k)$ is negative, meaning the group A agent in cell $(j,k)$ is moving toward $(j-1,k)$.
The presence of a group B agent in cell $(j-1,k)$ should induce a slowdown.
If cell $(j-1,k)$ is occupied by a group B agent but cell $(j,k)$ is not, then the velocity scaling is $c_{1}$.
The velocity scaling is $c_{3}$ if both $(j-1,k)$ and $(j,k)$ are occupied by agents from group B.

We assume that the velocity scalings satisfy the natural relationship
\begin{equation}\label{ccond}
c_3 < c_2 \leqslant c_1 < c_0.
\end{equation}
With respect to $V_{1}^{A} (j, k, \sigma^{B})$, inequalities~\eqref{ccond} reflect the fact that velocity scaling $c_0$ corresponds to local absence of agents from group B, velocity scalings $c_1$ and $c_2$ correspond to interaction with only one agent from group B, and velocity scaling $c_3$ corresponds to interaction with two agents from group B.

\subsection*{Transition probabilities for the stochastic model.}

The interaction velocity fields $V^A(j,k,\sigma^{B})$ and $V^B(j,k,\sigma^{A})$ incorporate~(\ref{li:floor-fields}) and~(\ref{li:slowdown}).
We use them to specify transition probabilities for an agent from group $A$ (or group $B$, respectively) at the current position $(j,k)$.
Agent motion must be horizontal or vertical~(\ref{li:no-diagonal}), and must be consistent with the corresponding interaction velocity field.
That is, the group $A$ agent can move left if $V_{1}^{A} (j, k, \sigma^{B}) < 0$, right if $V_{1}^{A} (j, k, \sigma^{B}) > 0$, down if $V_{2}^{A} (j, k, \sigma^{B}) < 0$, and up if $V_{2}^{A} (j, k, \sigma^{B}) > 0$.
We enforce the exclusion principle~(\ref{li:exclusion}) for agents in the same group.

Assuming this framework, the probability of transition $(j,k) \to (j \pm 1,k)$ for a member of group A during a small time 
interval $\Delta t$ is given by
\begin{align}
\label{PA}
	P^A_{(j,k) \to (j \pm 1,k)} &= \pm\Delta t \phi^A_1(j,k)H(\pm\phi^A_1(j,k))\sigma^A_{j,k}(1-\sigma^A_{j\pm 1, k}) \times \\
\nonumber
		&\quad \left[c_0(1-\sigma^B_{j,k})(1-\sigma^B_{j\pm 1,k}) + c_1(1-\sigma^B_{j,k})\sigma^B_{j\pm 1,k} + c_2\sigma^B_{j,k}(1-\sigma^B_{j\pm 1,k}) + c_3\sigma^B_{j,k}\sigma^B_{j\pm 1,k}\right],
\end{align}
while the probability of transition $(j,k) \to (j,k\pm 1)$ for a member of group A is given by
\begin{align}
\label{PB}
	P^A_{(j,k) \to (j,k \pm 1)} &= \pm\Delta t \phi^A_2(j,k)H(\pm\phi^A_2(j,k))\sigma^A_{j,k}(1-\sigma^A_{j, k \pm 1}) \times \\
\nonumber
		&\quad \left[c_0(1-\sigma^B_{j,k})(1-\sigma^B_{j,k\pm 1}) + c_1(1-\sigma^B_{j,k})\sigma^B_{j,k\pm 1} + c_2\sigma^B_{j,k}(1-\sigma^B_{j,k\pm 1}) + c_3\sigma^B_{j,k}\sigma^B_{j,k\pm 1}\right].
\end{align}
The diagonal entries of the transition probability matrix are defined in a standard manner as
\begin{equation*}
  P^{A}_{(j,k) \to (j,k)} = 1 - P^{A}_{(j,k) \to (j + 1,k)} - P^{A}_{(j,k) \to (j - 1,k)} - P^{A}_{(j,k) \to (j,k + 1)} - P^{A}_{(j,k) \to (j,k - 1)},
\end{equation*}
so that the entries in each row sum to $1$.
In~\eqref{PA} and \eqref{PB}, the Heaviside functions ensure that the probability of transition to a neighboring horizontal (vertical) cell is nonzero only if the horizontal (vertical) projection of the floor field points in the direction of the neighboring cell.
The $\sigma^A_{j,k}(1-\sigma^A_{j\pm 1, k})$ and $\sigma^A_{j,k}(1-\sigma^A_{j, k \pm 1})$ terms express the exclusion principle: The probability of transition is nonzero only if the current cell is occupied and the target cell is not occupied by a member of the same group. 

We assume periodic boundary conditions for simplicity.
It would be easy to specify, and interesting to analyze, more exotic boundary conditions.

Our stochastic model is $\bsigma(t) = (\sigma^{A} (t), \sigma^{B} (t))$, the finite-state, continuous-time Markov chain with transition probabilities given by~\eqref{PA} and \eqref{PB}.
Note that the dimension of the state space is quite high, as is typical for spatially extended stochastic processes.

\subsection*{Remarks.}

As explained before, agent motion must be consistent with the corresponding interaction velocity field - the angle between the random direction of motion and this field must be strictly less than $\pi / 2$.
When the floor fields assume the gradient form~\eqref{eq:phi} associated with potential~\eqref{eq:quadpsi}, each agent must move toward its target exit point (down the gradient).
It would be interesting to allow motion \textit{away} from the target exit point.
Such motion could be useful for obstacle avoidance and evacuation scenarios.

Our exclusion principle~(\ref{li:exclusion}) applies only to agents from the same group. 
This modeling issue has been treated various ways in the pedestrian dynamics literature.
When modeling car traffic, the exclusion principle can be applied to all agents (even from different groups) because all agents move in the same direction.
However, vehicle traffic and pedestrian traffic differ fundamentally.
Two groups of pedestrians can move toward each other and should be allowed to pass through one another.
Such behavior occurs at complex intersections and crossings, for example.

Existing work has focused on two pass-through mechanisms - (i) the side-stepping (or avoidance) mechanism~\cite{bksz01,kirch02,kirch03,guhu081}, and (ii) the slowdown mechanism~\cite{ckpt14}.
Both mechanisms are probably relevant for pedestrian flow, especially at intermediate and high densities.
Look-ahead mechanisms would add realism as well (\textit{e.g.}~\cite{soka06,degond13,hst14}).
In the present work, we have opted to focus on the slowdown mechanism.
Note that setting $c_1=c_2=c_3=0$ in our stochastic model would invoke the complete exclusion principle with respect to agents from both groups, thereby essentially disallowing the groups to pass through one another.

\subsection*{Some extensions.}

Before examining mesoscopic and macroscopic descriptions, we discuss generalizations.
First, our modeling framework naturally extends to lattices of any dimension and to at least three groups of interacting agents.
Second, the interaction mechanisms between different groups can be more complex.
For instance, agents can be allowed to switch between groups, either probabilistically or deterministically.
This mechanism is relevant when modeling obstacle avoidance, or evacuation scenarios with multiple exits. 
Third, our model can be extended to include \lq\lq{}chemical reaction\rq\rq{} mechanisms, as is done in the reaction diffusion master equation framework (see \textit{e.g.}~\cite{rdma1,rdma2,rdma3} for recent work).
This is relevant for problems wherein transport of reacting particles is of interest. 
We intend to explore these generalizations both theoretically and computationally in future work.
{"config": "arxiv", "file": "1806.01365/sec2stochmodel-rev-v2.tex"}