# A Multispecies Cross-Diffusion Model for Territorial Development

Abdulaziz Alsenafi<EMAIL_ADDRESS>Department of Mathematics,
Kuwait University, Kuwait  Alethea B. T. Barbaro<EMAIL_ADDRESS>Delft Institute of Applied Mathematics, Faculty of Electrical Engineering,
Mathematics and Computer Science, Delft University of Technology, The
Netherlands

##  1 Abstract

We develop an agent-based model on a lattice to investigate territorial
development motivated by markings such as graffiti, generalizing a previously-
published model to account for $K$ groups instead of two groups. We then
analyze this model and present two novel variations. Our model assumes that
agents’ movement is a biased random walk and that the agents’ movement is away
from rival groups’ markings. All interactions between agents is indirect,
mediated through the markings. We numerically demonstrate that in a system of
three groups, the groups segregate in certain parameter regimes. Starting from
the discrete model, we also formally derive the continuum system of $2K$
convection-diffusion equations for our model. These equations exhibit cross-
diffusion due to the avoidance of the rival groups’ markings. Linear stability
analysis is performed to study the phase transition that the system undergoes
as parameters are varied. Two new variations of the model are introduced and
numerically studied. Both exhibit a phase transition, with segregation
dynamics that are distinct both from each other and from the original model.

Keywords: phase transition, cross-diffusion, segregation model, pattern
formation, agent-based model

##  2 Introduction

Many types of organisms are known to exhibit territorially. Examples include
insects, fish, amphibians, reptiles, birds, and mammals [1], and, of course,
human beings [2]. Even plants could be considered to display this trait [3],
with some such as _Eucalyptus_ excreting a chemical that inhibits the growth
of other species [4]. Reasons for territorial behavior include protection of
breeding sites and access to resources, and territorial organisms have several
ways of claiming their territory. Two common ways are through some sort of
marking, either chemical or physical, and through direct confrontation.

In this paper, we focus on the case of territory formation for a mobile
species through territorial markings. We use the example of gangs of human
beings reacting to the graffiti of other gangs, but the model could equally
well be applied to other mobile species. In [5, 6], Moorcroft et al. modeled
how different packs of animals like coyotes and wolves base their movement on
scent marking. It was discovered that both coyotes and wolves use scent
marking to tag territories [7]. Once wolves or coyotes encounter foreign
scents, they in turn mark their territory with their scent and usually head
back to their own home territory. This dynamic was studied in detail in [8],
where it was found that different packs of wolves can live in the same region
without having contact with other packs, but each has its own territory.

Researchers studying gang dynamics based the gang movement dynamics on
existing ecological models of animal species that exhibit territorial
behavior. Smith et al. [9] combined the ecological model in [6] with Hegemann
et al.’s network model [10] to produce a model for gang territoriality. The
new model was then solved numerically, and their results were compared to real
data about gang territories in Los Angeles. In [11], Barbaro et al. used a
statistical mechanics approach to study how gang territories could be formed
based on graffiti. This also drew on ideas from coyote and wolf scent-marking.
The authors chose to use a spin system, which is similar to the Ising model
[12] that simulates ferromagnetism. A two-dimensional lattice was used, with
an agent and a graffiti spin at each site, and there were only indirect
interactions between the agent spins. The authors showed that their model
exhibits a phase transition in which gangs cluster together to form territory.

Other work on modeling criminal behavior has also been done. This work is
tangential to the work presented here, but is presented for the interested
reader. Clustering methods have been used to study gang affiliations between
gang members in Los Angeles [13]. Network models are also used to study gangs.
In [10], the authors present an agent-based model, which is coupled with a
rivalry network to explore how gang rivalries are formed. In [14], a model for
burglary was developed, and a continuum system consisting of a coupled
reaction-diffusion equations was derived. This and similar systems were
analyzed in [15, 16], and [17]. Modifications of the model were explored in
[18, 19]. The reaction diffusion-diffusion equations were analyzed further in
[20], which showed that they exhibit similar behavior to chemotactic systems
with cross-diffusion. Recently, in [21], Wang et al. extended the burglary
model by including independent Poison clocks in the time steps; a martingale
with both a deterministic and a stochastic part was derived and analyzed. Work
on riots and social segregation have also followed from this line of research
[22, 23]. For a more in-depth review of the crime modeling literature, the
reader is referred to [24].

Our paper is based on the work of [25], wherein the authors performed a
bottom-up approach similar to the Moorcroft model [6] to produce a discrete
system to describe gang territorial development and formally derive from it
the following system of convection-diffusion equations:

$\begin{cases}\displaystyle\frac{\partial\xi_{A}}{\partial
t}(x,y,t)=\gamma\rho_{A}(x,y,t)-\lambda\xi_{A}(x,y,t)\\\
\displaystyle\frac{\partial\xi_{B}}{\partial
t}(x,y,t)=\gamma\rho_{B}(x,y,t)-\lambda\xi_{B}(x,y,t)\\\
\displaystyle\frac{\partial\rho_{A}}{\partial
t}(x,y,t)=\frac{D}{4}\nabla\cdot\left[\nabla\rho_{A}(x,y,t)+2\beta\left(\rho_{A}(x,y,t)\nabla\xi_{B}(x,y,t)\right)\right]\\\
\displaystyle\frac{\partial\rho_{B}}{\partial
t}(x,y,t)=\frac{D}{4}\nabla\cdot\left[\nabla\rho_{B}(x,y,t)+2\beta\left(\rho_{B}(x,y,t)\nabla\xi_{A}(x,y,t)\right)\right],\end{cases}$
(1)

where $\xi_{i}$ is the graffiti density of gang $i$ and $\rho_{i}$ is the
agent density from gang $i$. The model undergoes a phase transition from no
territorial development to distinct territorial formation as the parameter
$\beta$ is changed. This phase transition was found both in the discrete and
continuum level. The authors there only considered the case of two gangs. In
contrast, in the present work, we will generalize the model and results of
[25] to consider any finite number of gangs. A modified version of the
convection-diffusion system in [25] was analyzed in [26], where they proved a
weak stability result and identified equilibrium solutions; interestingly,
though, they did not find segregated solutions in this modified system.

This paper’s outline is as follows: In Section 2.1, we introduce our extension
of the original two-gang agent-based model [25]. The rest of the article is
based upon this extension. We next define an order parameter in Section 2.2
that will be used to analyze our system’s different states and characterize
phase transitions. In Section 3, we present the results of a special case of
our discrete model numerical simulations as well as show our analysis of the
systems’ phase transitions. In Sections 4.1 and 4.2, we will derive the
general continuum limit from the discrete model. In Sections 4.3 and 5, we
will derive a steady-state solution for our continuum model and perform linear
stability analysis to determine whenever the well-mixed solution becomes
unstable. In Section 6, we introduce and study two variations on the model,
where parameter $\beta$ is made gang-dependent. Finally, in Section 7, we
conclude with a discussion of the results and open problems.

###  2.1 Discrete Model

In this paper, we extend and generalize the interacting particle model in [25]
to now include $K$ gangs as opposed to only two, keeping all other dynamics
similar. We shall use a square lattice $S$ of size $L\times L$, with periodic
boundary conditions and area $1$. We assume that we have $K$ gangs, $1,2,\dots
K$, and the number of agents belonging in each gang $j$ is denoted by $N_{j}$.
The systems’ total number of agents is denoted by $N$:

$N=\sum_{\begin{subarray}{c}i=1\\\ i\neq j\end{subarray}}^{K}N_{i}$

These agents are distributed over the lattice. Our model allows multiple gang
agents regardless of their gang affiliation to be on the same site. We denote
the number of agents of gang $j$ at site $(x,y)$ at time $t$ by $n_{j}(x,y,t)$
and their densities are $\rho_{j}(x,y,t)=\frac{n_{j}(x,y,t)}{l^{2}}$, where
$l=\frac{1}{L^{2}}$ is the lattice spacing. The amount of graffiti belonging
to gang $j$ at site $(x,y)$ on time $t$ is denoted by $g_{j}(x,y,t)$. We
denote the graffiti density of gang $j$ by
$\xi_{j}(x,y,t)=\frac{g_{j}(x,y,t)}{l^{2}}$.

Our model assumes that every agent has to move at every time step to one of
their four neighboring sites, which are the sites up, down, to the left, and
to the right of it. That is, an agent currently occupying site $(x,y)$ would
move to an element of the set of sites
$\\{(x+l,y),(x-l,y),(x,y+l),(x,y-l)\\}$. The neighboring sites of $(x,y)$ will
be denoted by $(\tilde{x},\tilde{y})\sim(x,y)$. In our model, each agent
performs a biased random walk, trying to avoid the opposing gangs graffiti.
Following [25], our model assumes that every agent puts down its own gang’s
graffiti on the lattice and this graffiti discourages the movement of agents
from a different gang onto that lattice site. However, now that we are
considering more than two gangs, each gang must avoid the graffiti of all
other gangs, leading us to define the _opposition sum_ for gang $j$ at site
$(x,y)$ at time $t$:

$\psi_{j}(x,y,t):=\sum_{\begin{subarray}{c}i=1\\\ i\neq
j\end{subarray}}^{K}\xi_{i}(x,y,t).$ (2)

We use this opposition sum to inform the movement dynamics of the agents. The
probability of an agent from gang $j$ to move from site
$s_{1}=(x_{1},y_{1})\in S$ to one of the neighboring sites
$s_{2}=(x_{2},y_{2})\in S$ is defined to be

$M_{j}(x_{1}\rightarrow x_{2},y_{1}\rightarrow
y_{2},t)=\frac{e^{-\beta\psi_{j}(x_{2},y_{2},t)}}{\sum\limits_{(\tilde{x},\tilde{y})\sim(x_{1},y_{1})}e^{-\beta\psi_{j}(\tilde{x},\tilde{y},t)}},$
(3)

again defined analogously to [25]. Here the parameter $\beta$ encodes the
strength of the avoidance of other gangs’ graffiti. As our model assumes that
all of the agents must move at each time step, it is easily seen that

$\sum\limits_{(\tilde{x},\tilde{y})\sim(x,y)}M_{j}(x\rightarrow\tilde{x},y\rightarrow\tilde{y},t)=1.$
(4)

The expected density of gang $j$ is therefore

$\displaystyle\rho_{j}(x,y,t+\delta t)=$
$\displaystyle\rho_{j}(x,y,t)+\sum_{(\tilde{x},\tilde{y})\sim(x,y)}\rho_{A}(\tilde{x},\tilde{y},t)M_{j}(\tilde{x}\rightarrow
x,\tilde{y}\rightarrow y,t)$
$\displaystyle-\rho_{j}(x,y,t)\sum_{(\tilde{x},\tilde{y})\sim(x,y)}M_{j}(x\rightarrow\tilde{x},y\rightarrow\tilde{y},t)$
$\displaystyle=$
$\displaystyle\sum_{(\tilde{x},\tilde{y})\sim(x,y)}\rho_{A}(\tilde{x},\tilde{y},t)M_{j}(\tilde{x}\rightarrow
x,\tilde{y}\rightarrow y,t).$ (5)

For the graffiti density update rules, each agent adds graffiti at its current
site with probability $\gamma$. It is also assumed that the graffiti decays at
every site with a rate of $\lambda$. Both the graffiti addition and decay are
scaled by the time step $\delta t$. Therefore, the graffiti density at site
$(x,y)\in S$ at time $t+\delta t$ is

$\xi_{j}(x,y,t+\delta t)=\xi_{j}(x,y,t)-\left(\delta
t\cdot\lambda\right)\xi_{j}(x,y,t)+\left(\delta
t\cdot\gamma\right)\rho_{j}(x,y,t).$ (6)

In all of our simulations, we initially randomly distribute the agents’
locations using the multivariate uniform distribution on the lattice $S$. We
also assume that the lattice is initially empty of graffiti.

###  2.2 Phases and an Order Parameter

In our simulations, we shall observe two phases, the well-mixed phase and the
segregated phase. These phases, we will see, are determined by parameter
$\beta$, introduced in (3). In the well-mixed phase, the agents are
distributed randomly throughout the lattice and their movement approximates a
random walk. However, for the segregated phase, the agents’ movement is a
biased random walk, and the agents form territories by clustering together. In
this section, we shall define an order parameter and use it to quantify these
different phases.

####  2.2.1 Expected Agent Density

We first compute the expected agent density for the well-mixed state at site
$(x,y)$ for each gang $j$. In this phase, the agents from each gang are
uniformly spread over the whole lattice $S$. Thus, the expected agent density
for gang $j$ at any given site is:

$\displaystyle E\left(\rho_{j}\right)$ $\displaystyle=\sum_{(x,y)\in
S}\rho_{j}(x,y)\times\frac{1}{L^{2}}$ $\displaystyle=\sum_{(x,y)\in
S}\frac{n_{j}(x,y)}{l^{2}}\times\frac{1}{L^{2}}$ $\displaystyle=\sum_{(x,y)\in
S}n_{j}(x,y)$ $\displaystyle=N_{j}.$ (7)

For the segregated phase, the agents are entirely separated into different
territories, each occupied by a distinct gang. To determine the expected agent
density for the segregated phase, we will need the following definitions and
assumptions. We define the territory $S_{j}$ to be the set of all sites that
are dominated by gang $j$ agents, i.e. all sites which have more gang $j$
agents than agents of another type. We also note that $S_{j}$ need not be
connected. We also define the area of sublattice $S_{j}$ by $R_{j}$, which is
the area dominated by gang $j$. We assume that the agents from each gang are
uniformly distributed in their territory, and that all sites are occupied by
agents; hence, every site is assumed to contain agents from exactly one gang.
These assumptions are validated in our simulation in Section 3. Accordingly,
$\\{S_{j}\\}_{j=1}^{K}$ form a partition for the lattice $S$ for
$j=1,2,\dots,K$.

Under these assumptions, we now calculate the expected agent density for the
agent density for gang $j$ in the segregated state by splitting the lattice
$S$ into two disjoint territories, $S_{j}$ and its complement ${S_{j}}^{c}$.
Calculating the expected agent density within the complement of gang $j$’s
territory easily finds:

$\displaystyle E\left(\rho_{j}\right)$
$\displaystyle=\sum_{(x,y)\in{S_{j}}^{c}}\frac{\rho_{j}(x,y)}{{R_{j}}^{c}}=0$

This last equality follows because all agents for gang $j$ are assumed to be
in $S_{j}$ in the (perfectly) segregated phase, hence none are in
${S_{j}}^{c}$. Next, calculating the expected agent density within $S_{j}$
gives us:

$\displaystyle E\left(\rho_{j}\right)$ $\displaystyle=\sum_{(x,y)\in
S_{j}}\frac{\rho_{j}(x,y)}{R_{j}}$ $\displaystyle=\sum_{(x,y)\in
S_{j}}\frac{n_{j}(x,y)}{l^{2}\times R_{j}}$
$\displaystyle=\frac{N_{j}}{l^{2}R_{j}}$

Thus, the expected agents density within $S_{j}$ is

$\displaystyle E\left(\rho_{j}(x,y)\right)$
$\displaystyle=\begin{cases}\frac{N_{j}}{l^{2}R_{j}},&(x,y)\in S_{j}\\\
0,&(x,y)\in{S_{j}}^{c}.\end{cases}$ (8)

If we further assume that the areas dominated by each gang are almost equal
and that there are an equal number of agents in each gang, we deduce that
$R_{j}=\frac{1}{K}$, where $K$ is the number of gangs. Under these
assumptions, the expected density in a segregated state for an agent from gang
$j$ is

$\displaystyle E\left(\rho_{j}(x,y)\right)$
$\displaystyle=\begin{cases}\frac{N_{j}}{l^{2}R_{j}},&(x,y)\in S_{j}\\\
0,&(x,y)\in{S_{j}}^{c}\end{cases}$
$\displaystyle=\begin{cases}\frac{KN_{j}}{l^{2}},&(x,y)\in S_{j}\\\
0,&(x,y)\in{S_{j}}^{c}.\end{cases}$

####  2.2.2 An Order Parameter

To investigate the phase transition, we define the following order parameter:

$\displaystyle\mathcal{E}(t)=~{}\frac{1}{4(K-1)}\left(\frac{1}{LN}\right)^{2}\sum_{j=1}^{K}~{}~{}\sum_{i>j}^{K}\sum_{(x,y)\in
S}\sum_{(\tilde{x},\tilde{y})\sim(x,y)}$
$\displaystyle[\left(\rho_{j}(x,y,t)-\rho_{i}(x,y,t)\right)\times$
$\displaystyle\left(\rho_{j}(\tilde{x},\tilde{y},t)-\rho_{i}(\tilde{x},\tilde{y},t)\right)].$
(9)

This order parameter is modeled after the order parameter in [25] and the
Hamiltonian function for the Ising Model [27, 12]. In terms of our model, our
order parameter becomes more positive if a site and its neighbors are
dominated by the same gang and becomes more negative when a site and its
neighbors are dominated by different gangs. It is approximately zero if none
of the gangs are dominating territory.

This is due to the fact that in the segregated, phase the agents cluster
together and this leads to there being only one gang present at site $(x,y)$.
This makes the term inside the sum in equation (9) to have a large magnitude;
if the same is true at the neighboring site, the second term inside the sum is
identical and once we multiply the two together, the resulting value would be
a positive number. However, whenever agents from all gangs become uniformly
distributed throughout the lattice, this results that the two sets of
parenthesis tend to be very small, and sometimes positive and sometimes
negative. Thus, after summing over the whole lattice and all the gangs, the
order parameter ends up near zero.

We now calculate an approximation for the order parameter when the phases are
well-mixed and when they are segregated. For simplicity, in this subsection
and all of our simulations, we consider the special case of three gangs and we
assume that the number of agents in each gang is equal, so that
$N_{j}=\frac{N}{3}$ for $j$ in $\\{1,2,3\\}$. The order parameter for this
special case is

$\displaystyle\mathcal{E}=\frac{1}{8}\left(\frac{1}{LN}\right)^{2}\sum_{(x,y)\in
S}\sum_{(\tilde{x},\tilde{y})\sim(x,y)}\biggr{[}$
$\displaystyle\left(\rho_{1}(x,y)-\rho_{2}(x,y)\right)\left(\rho_{1}(\tilde{x},\tilde{y})-\rho_{2}(\tilde{x},\tilde{y})\right)$
$\displaystyle+$
$\displaystyle\left(\rho_{1}(x,y)-\rho_{3}(x,y)\right)\left(\rho_{1}(\tilde{x},\tilde{y})-\rho_{3}(\tilde{x},\tilde{y})\right)$
$\displaystyle+$
$\displaystyle\left(\rho_{2}(x,y)-\rho_{3}(x,y)\right)\left(\rho_{2}(\tilde{x},\tilde{y})-\rho_{3}(\tilde{x},\tilde{y})\right)\biggr{]}.$
(10)

In the well-mixed state, based on equation (7) and on our assumptions that the
agents from all gangs are uniformly distributed and that each lattice site has
four neighbors, our equation simplifies to

$\displaystyle\mathcal{E}=\frac{1}{8}\left(\frac{1}{LN}\right)^{2}\sum_{(x,y)\in
S}\biggr{[}$ $\displaystyle
4\left(N_{1}-N_{2}\right)\left(N_{1}-N_{2}\right)+4\left(N_{1}-N_{3}\right)\left(N_{1}-N_{3}\right)$
$\displaystyle+$ $\displaystyle
4\left(N_{2}-N_{3}\right)\left(N_{2}-N_{3}\right)\biggr{]}$

Simplifying the terms in the brackets yields that

$\displaystyle\mathcal{E}$
$\displaystyle=\frac{4}{8}\left(\frac{1}{LN}\right)^{2}\sum_{(x,y)\in
S}\biggr{[}\left(N_{1}-N_{2}\right)^{2}+\left(N_{1}-N_{3}\right)^{2}+\left(N_{2}-N_{3}\right)^{2}\biggr{]}$
$\displaystyle=\frac{1}{2}\left(\frac{1}{N}\right)^{2}\biggr{[}\left(N_{1}-N_{2}\right)^{2}+\left(N_{1}-N_{3}\right)^{2}+\left(N_{2}-N_{3}\right)^{2}\biggr{]}.$

However, since we assumed that the number of agents from each gang
$N_{1},N_{2},$ and $N_{3}$ are equal, it follows easily that the order
parameter for the agents in a well-mixed phase is

$\mathcal{E}\approx 0.$ (11)

We will next calculate the order parameter for the segregated phase. Here, we
assume a perfectly segregated phase and split the lattice $S$ into the three
regions $S_{1}$, $S_{2}$ and $S_{3}$ belonging to each gang, which gives us:

$\displaystyle\mathcal{E}=\frac{1}{8}\left(\frac{1}{LN}\right)^{2}\Biggr{[}\sum_{(x,y)\in
S_{1}}\sum_{(\tilde{x},\tilde{y})\sim(x,y)}\biggr{[}$
$\displaystyle\left(\rho_{1}(x,y)-\rho_{2}(x,y)\right)\left(\rho_{1}(\tilde{x},\tilde{y})-\rho_{2}(\tilde{x},\tilde{y})\right)$
$\displaystyle+$
$\displaystyle\left(\rho_{1}(x,y)-\rho_{3}(x,y)\right)\left(\rho_{1}(\tilde{x},\tilde{y})-\rho_{3}(\tilde{x},\tilde{y})\right)$
$\displaystyle+$
$\displaystyle\left(\rho_{2}(x,y)-\rho_{3}(x,y)\right)\left(\rho_{2}(\tilde{x},\tilde{y})-\rho_{3}(\tilde{x},\tilde{y})\right)\biggr{]}.$
$\displaystyle+\sum_{(x,y)\in
S_{2}}\sum_{(\tilde{x},\tilde{y})\sim(x,y)}\biggr{[}$
$\displaystyle\left(\rho_{1}(x,y)-\rho_{2}(x,y)\right)\left(\rho_{1}(\tilde{x},\tilde{y})-\rho_{2}(\tilde{x},\tilde{y})\right)$
$\displaystyle+$
$\displaystyle\left(\rho_{1}(x,y)-\rho_{3}(x,y)\right)\left(\rho_{1}(\tilde{x},\tilde{y})-\rho_{3}(\tilde{x},\tilde{y})\right)$
$\displaystyle+$
$\displaystyle\left(\rho_{2}(x,y)-\rho_{3}(x,y)\right)\left(\rho_{2}(\tilde{x},\tilde{y})-\rho_{3}(\tilde{x},\tilde{y})\right)\biggr{]}$
$\displaystyle+\sum_{(x,y)\in
S_{3}}\sum_{(\tilde{x},\tilde{y})\sim(x,y)}\biggr{[}$
$\displaystyle\left(\rho_{1}(x,y)-\rho_{2}(x,y)\right)\left(\rho_{1}(\tilde{x},\tilde{y})-\rho_{2}(\tilde{x},\tilde{y})\right)$
$\displaystyle+$
$\displaystyle\left(\rho_{1}(x,y)-\rho_{3}(x,y)\right)\left(\rho_{1}(\tilde{x},\tilde{y})-\rho_{3}(\tilde{x},\tilde{y})\right)$
$\displaystyle+$
$\displaystyle\left(\rho_{2}(x,y)-\rho_{3}(x,y)\right)\left(\rho_{2}(\tilde{x},\tilde{y})-\rho_{3}(\tilde{x},\tilde{y})\right)\biggr{]}\Biggr{]}$

Using equation (8), we substitute the expectation of each $\rho_{i}$ for each
region, which gives us the following approximation:

$\displaystyle\mathcal{E}\approx\frac{1}{8}\left(\frac{1}{LN}\right)^{2}\Biggr{[}$
$\displaystyle\sum_{(x,y)\in
S_{1}}\biggr{[}4\left(\frac{N_{1}}{l^{2}R_{1}}\right)^{2}+4\left(\frac{N_{1}}{l^{2}R_{1}}\right)^{2}\biggr{]}$
$\displaystyle+$ $\displaystyle\sum_{(x,y)\in
S_{2}}\biggr{[}4\left(\frac{-N_{2}}{l^{2}R_{2}}\right)^{2}+4\left(\frac{N_{2}}{l^{2}R_{2}}\right)^{2}\biggr{]}$
$\displaystyle+$ $\displaystyle\sum_{(x,y)\in
S_{3}}\biggr{[}4\left(\frac{-N_{3}}{l^{2}R_{3}}\right)^{2}+4\left(\frac{-N_{3}}{l^{2}R_{3}}\right)^{2}\biggr{]}\Biggr{]}.$

Further simplifying yields:

$\displaystyle\mathcal{E}$
$\displaystyle\approx\frac{1}{8}\left(\frac{1}{LN}\right)^{2}\left(\frac{8}{l^{4}}\right)\left[\sum_{(x,y)\in
S_{1}}\left(\frac{N_{1}}{R_{1}}\right)^{2}+\sum_{(x,y)\in
S_{2}}\left(\frac{N_{2}}{R_{2}}\right)^{2}+\sum_{(x,y)\in
S_{3}}\left(\frac{N_{3}}{R_{3}}\right)^{2}\right]$
$\displaystyle=\left(\frac{1}{LN}\right)^{2}\left(\frac{1}{l^{4}}\right)\left[R_{1}\left(\frac{N_{1}}{R_{1}}\right)^{2}+R_{2}\left(\frac{N_{2}}{R_{2}}\right)^{2}+R_{3}\left(\frac{N_{3}}{R_{3}}\right)^{2}\right]$
$\displaystyle=\left(\frac{L}{N}\right)^{2}\left[\frac{N_{1}^{2}}{R_{1}}+\frac{N_{2}^{2}}{R_{2}}+\frac{N_{3}^{2}}{R_{3}}\right].$

By further assuming that all gangs have the same number of agents $N_{i}=N/3$
and that the regions have the same area size $R_{i}=\frac{L^{2}}{3}$, the
previous equation can be further simplified to

$\mathcal{E}\approx\left(\frac{L}{N}\right)^{2}\left(\frac{\left(N/3\right)^{2}}{L^{2}/3}+\frac{\left(N/3\right)^{2}}{L^{2}/3}+\frac{\left(N/3\right)^{2}}{L^{2}/3}\right)$

By further simplifying the equation, we easily see that

$\mathcal{E}\approx 1.$ (12)

Therefore, the order parameter for the perfectly segregated system is
approximately equal to one.

##  3 Simulations of the Discrete Model

We now will present the results of the simulations of our discrete model. For
simplicity, unless otherwise in our simulations we assume we only have $three$
gangs $1,2$ and $3$, and that all gangs are assumed to have $50,000$ agents.
We shall also assume that the lattice size $L\times L$ is $100\times 100$ with
lattice spacing 2 $l=1$, and we will use $100,000$ time steps with each step
size $\delta t=1$.

###  3.1 Well-Mixed State

We start our simulations with $\beta=5\times 10^{-6}$ and the resulting
lattice simulations are visualized in Figure 1. The first two lattices in
Figure 1 represent the time evolution of agent density, whereas the last two
lattices represent the graffiti density over time. We assign the colors red,
green and blue for gangs $1$, $2$ and $3$ respectively. The color white is
used if there are the same number of agents or the same amount of graffiti
from all gangs at a site. The colors cyan, magenta and yellow are used if a
site has two gangs (blue and green, blue and red, or red and green,
respectively). Finally, if the site is empty, then it will be assigned the
color black.

Figure 1: Agent (left two) and graffiti (right two) densities’ temporal
evolution for a well-mixed state. Here we have $N_{1}=N_{2}=N_{3}=50,000$,
with $\lambda=\gamma=0.5$, $\beta=5\times 10^{-6}$, $\delta t=1$ and the
lattice size is $100\times 100$. Note that the initial graffiti lattice
appears black because it is empty. The final graffiti lattice appears white
because all sites have (almost) the same graffiti densities from all three
gangs. It is clear from this figure that the agents remain well mixed over
time.

From Figure 1, we clearly see that the gangs remain well mixed over time for
$\beta=5\times 10^{-6}$. We do not see any patterns being formed for the
graffiti, with the initial graffiti lattice black and the final graffiti
lattice white, and the gang agents’ movement is in essence a two-dimensional
random walk, hardly taking the opposing gangs graffiti into consideration due
to the low $\beta$ value. This is due to the way the agents are allowed to
move in equation (3), where a very small $\beta$ values give the agent a
probability of nearly $0.25$ to move to one of the four neighboring sites.

###  3.2 Segregated State

The value of $\beta$ is now increased so that it is equal to $3\times 10^{-5}$
and we keep all other parameters the same. The resulting lattice is visualized
in Figure 2. The top row illustrates the time evolution of agent density,
whereas the bottom row shows the graffiti density over time.

Figure 2: Agent (top) and graffiti (bottom) densities temporal evolution for a
segregated state. Here ,we have $N_{1}=N_{2}=N_{3}=50,000$, with
$\lambda=\gamma=0.5$, $\beta=3\times 10^{-5}$, $\delta t=1$ and the lattice
size is $100\times 100$. We see that the agents segregate into distinct
territories, coarsening over time.

From Figure 2, we see that initially the agents are well-mixed. However, as
time evolves, we see that agents from each gang cluster together to form all-
red, all-green and all-blue territories. As time increases, the patterns in
both the agent and graffiti densities coarsen. From the same figure, we
clearly see that the graffiti density is similar to the agents density and the
agents movements in the state is based on the other gangs graffiti. That is,
in this state, the $\beta$ value is large enough that the agents are reacting
to the opposing gang graffiti and the agents movement is no longer an unbiased
random walk. We can also observe from this figure that the areas with more
than one gang’s graffiti lie at the boundaries of the territories dominated by
each gang. Similarly, this is where we observe the agents overlapping, though
to a lesser extent. Presumably, this overlap enables the coarsening see in the
figure.

###  3.3 System Parameters and the Discrete Phase Transition

####  3.3.1 Effects of $\beta$

In Sections 3.1 and 3.2, we saw that changing the value of the parameter
$\beta$ could lead to a phase transition. In order for us to study the phase
transitions, we use the concept of order parameter that we introduced in
Section 2.2.2. In equation (10), we defined an order parameter for a system of
three gangs to be

$\displaystyle\mathcal{E}=\frac{1}{8}\left(\frac{1}{LN}\right)^{2}\sum_{(x,y)\in
S}\sum_{(\tilde{x},\tilde{y})\sim(x,y)}\biggr{[}$
$\displaystyle\left(\rho_{1}(x,y)-\rho_{2}(x,y)\right)\left(\rho_{1}(\tilde{x},\tilde{y})-\rho_{2}(\tilde{x},\tilde{y})\right)$
$\displaystyle+$
$\displaystyle\left(\rho_{1}(x,y)-\rho_{3}(x,y)\right)\left(\rho_{1}(\tilde{x},\tilde{y})-\rho_{3}(\tilde{x},\tilde{y})\right)$
$\displaystyle+$
$\displaystyle\left(\rho_{2}(x,y)-\rho_{3}(x,y)\right)\left(\rho_{2}(\tilde{x},\tilde{y})-\rho_{3}(\tilde{x},\tilde{y})\right)\biggr{]}.$
(13)

This order parameter is defined to have a low value for a well-mixed phase and
high value for a segregated phase; in Section 2.2.2, we saw that the order
parameter $\mathcal{E}\approx 0$ for a well mixed state and
$\mathcal{E}\approx 1$ for a fully segregated state. For our simulations, we
graphed the order parameter over the course of the simulation for different
values of $\beta$, visualizing the output in Figure 3.

Figure 3: How changing the $\beta$ parameter affects the system. Here we have
$N_{1}=N_{2}=N_{3}=50,000$, with $\lambda=\gamma=0.5$ and the lattice size is
$100\times 100$. (Left) It is seen that for a small $\beta$ value the system
remains well-mixed and the order parameter is almost zero over all time steps.
For larger $\beta$ values, we see that the order parameter increases quickly
as the system segregates and levels off to around one. (Right) At the final
time step, we take the order parameter value for different $\beta$ values. We
clearly see that as the $\beta$ value increases there is a critical $\beta$
value at which a phase transition occurs.

We see in the left plot in Figure 3, the time evolution of the order parameter
for different $\beta$ values. Here, we easily see that given enough time
steps, the order parameter levels off to a certain value, presumably its
asymptotic value. We see that for $\beta=0$ and $\beta=0.000005$, the order
parameter remained approximately zero throughout all time steps. This is
expected as the system remains well-mixed for these relatively small $\beta$
values. However, we see that once we increase the values of $\beta$, then the
order parameter starts to increase. For instance, if $\beta=0.000025$ or
$\beta=0.000030$, then the order parameter increases fairly quickly in the
first $10,000$ time steps before leveling off to just under the fully-
segregated value of $1$ for the remaining time steps. This shows us that for
these relatively large $\beta$ values, the system segregates fairly quickly
and remains segregated throughout the simulation. Finally, we also see that if
we choose $\beta=0.00001$ then the order parameter does increase and the
system does exhibit some segregation, but this is not perfect segregation as
the order parameter levels off to around $0.4$.

It is evident from the left plot of Figure 3 that there is a critical $\beta$
in which the system undergoes a phase transition. We define the critical
$\beta$ to be the value where the order parameter is equal to $0.01$, and
denote it by $\beta^{*}$. To find the value of $\beta^{*}$, we take the final
value of the order parameter and plot it against different $\beta$ values. The
output is then visualized on the right plot of Figure 3. From that plot, we
can see that the phase transition occurs when
$\beta^{*}\in(0.000005,0.000006)$.

####  3.3.2 Effects of Other Parameters

In order to investigate how other parameters such as system mass, time step,
lattice size, graffiti rate and decay rate affect the system phase transition,
we vary one parameter at a time while keeping all other parameters fixed. This
is important since in the derivation of the continuum equations for our
system, we will assume that both the time step $\delta t$ and the lattice
spacing $l$ approach zero. It is therefore essential to know if a finer grid
affects our discrete model as opposed to a coarser grid. We also would like to
know if taking smaller or bigger time steps might affect the rate of
segregation and if it has any effect on the phase transition.

We begin by studying how the time step might affect the system. To do that, we
keep all our system parameters constant and decrease the time step $\delta t$
from $1$ to $0.1$; we then plot the final order parameter value for different
$\beta$ values. The results are visualized on the right plot in Figure 3. In
the plot, it is clear that the smaller time step does not affect the rate of
segregation, nor does it affect where the phase transition occurs.

We were also interested in how the mass might affect our system. In Figure 4,
we see in the first plot that when the mass is $75,000$ the critical $\beta$
at which the phase transition occurs is about $1.8\times 10^{-5}$. However, in
the middle plot, the mass is increased to $150,000$ and this time the phase
transition occurs around $0.9\times 10^{-5}$. Thus, we notice that as the
systems’ mass increases, the resulting phase transition happens at a smaller
$\beta$ value. Physically, this makes sense, since having a larger number of
agents implies that there will be more graffiti being added at each site and
thus a smaller $\beta$ value should be sufficient for the agents to react to
the graffiti field.

Figure 4: The Effect of Parameters on the Phase Transition. Here we have
$N_{1}=N_{2}=N_{3}=50,000$, with $\lambda=\gamma=0.5,\delta t=1$ and the
lattice size is $100\times 100$. In all three plots, it is seen that for small
$\beta$ values, the system remains well-mixed and the resulting order
parameter values is approximately zero over time. However, for larger $\beta$
values we see that the order parameter increases quickly as the system
segregates and levels off to around one.The order parameter value is taken at
the final time step for different $\beta$ values. We clearly see that as the
$\beta$ value increases there is a critical $\beta$ value in which a phase
transition occurs. The three plots show the effects of changing the mass and
the ratio $\gamma/\lambda$.

We also investigated how the ratio $\frac{\gamma}{\lambda}$ might change where
the phase transition occurs. Again, we kept all other parameters fixed and
changed the value of the ratio by altering the decay rate $\lambda$. Having a
higher decay rate means that the graffiti is decaying more quickly and thus
each site would have less graffiti. We found that by decreasing the
$\frac{\gamma}{\lambda}$ ratio, a higher $\beta$ value is needed for
segregation. This is evident in the middle and right plots in Figure 4. There,
we clearly see that when the $\frac{\gamma}{\lambda}=1$, the critical $\beta$
is around $0.9\times 10^{-5}$, whereas when the $\frac{\gamma}{\lambda}$ is
decreased to $0.5$, the critical $\beta$ is around $1.8\times 10^{-5}$.
Physically, this is due to the fact that less graffiti on a site means that a
larger $\beta$ value is necessary for the agents to react to it.

We were also interested in investigating how changing the grid size might
affect the segregation when we alter the other parameters. In Figures 3 and 4,
we see that increasing the number of sites from $L=50$ to $L=100$ does not
have any noticeable effect on our discrete model.

##  4 Deriving the Convection-Diffusion System

In this section, we will formally derive the continuum equations of our system
and will prove that the limiting system of convection-diffusion equations is

$\begin{cases}\displaystyle\frac{\partial\xi_{j}}{\partial
t}(x,y,t)=\gamma\rho_{j}(x,y,t)-\lambda\xi_{j}(x,y,t)\\\
\displaystyle\frac{\partial\rho_{j}}{\partial
t}(x,y,t)=\frac{D}{4}\nabla\cdot\left[\nabla\rho_{j}(x,y,t)+2\beta\left(\rho_{j}(x,y,t)\nabla\left(\sum_{\begin{subarray}{c}i=1\\\
i\neq j\end{subarray}}^{K}\xi_{i}(x,y,t)\right)\right)\right]\end{cases}$ (14)

with periodic boundary conditions, where $j\in\\{1,2,\dots,K\\}$. Since our
discrete model is a multiple-gang extension of the two-species model in [25],
we proceed with finding the continuum equations by following the steps of the
derivation of the continuum model therein. With minor modifications, the same
derivation goes through for this multiple-gang case.

Deriving the continuum limits from discrete models is of great interest to the
mathematical community; for example, even just from the crime modeling
literature, we can refer you to papers [28, 18, 29, 14]. These continuum
equations are often formally derived by assuming appropriate smoothness of the
gang density and graffiti density and taking both the grid spacing and time
step to zero, as we will do here. The continuum partial differential equations
give us more tools for understanding the macroscopic behavior of the model.

###  4.1 Continuum Graffiti Density

We start by formally deriving the continuum equations for graffiti densities,
recalling that for $j\in\\{1,2,\dots,K\\}$, the discrete model (6):

$\xi_{j}(x,y,t+\delta t)=\xi_{j}(x,y,t)-\delta
t\cdot\lambda\cdot\xi_{j}(x,y,t)+\delta t\cdot\gamma\cdot\rho_{j}(x,y,t).$

Rearranging the equation and dividing by $\delta t$ gives us:

$\frac{\xi_{j}(x,y,t+\delta t)-\xi_{j}(x,y,t)}{\delta
t}=\gamma\cdot\rho_{j}(x,y,t)-\lambda\cdot\xi_{j}(x,y,t).$

This is now in the form of a difference equation. Assuming sufficient
smoothness of the agent and graffiti densities $\rho_{j}$ and $\xi_{j}$, we
take $\delta t\rightarrow 0$. This gives us the final form of the graffiti
continuum equation for gang $j$:

$\frac{\partial\xi_{j}}{\partial
t}(x,y,t)=\gamma\rho_{j}(x,y,t)-\lambda\xi_{j}(x,y,t).$ (15)

###  4.2 Continuum Agent Density

####  4.2.1 Tools for the Derivation

Deriving the continuum equations for the agent densities is more complex, so
before we begin, we define several quantities that will be useful in the
derivation. We will need equation (2), which we recall here:

$\psi_{j}(x,y,t):=\sum_{\begin{subarray}{c}i=1\\\ i\neq
j\end{subarray}}^{K}\xi_{i}(x,y,t).$

Employing this notation, the first quantity we define is

$T_{j}(x,y,t):=\frac{e^{\beta\psi_{j}(x,y,t)}}{4+l^{2}\left(\beta^{2}\left(\nabla\psi_{j}(x,y,t)\right)^{2}-\beta\Delta\psi_{j}(x,y,t)\right)}.$
(16)

We will use $T_{j}$ to account for the influences of the neighbors and
neighbor’s neighbors in the discrete model.

Next, we derive approximations to $\nabla T_{j}$ and $\Delta T_{j}$, which we
will use later in this section. For simplicity, the notation $(x,y,t)$ will be
dropped as there will be no neighbors $(\tilde{x},\tilde{y},t)$ in the
derivation of these quantities. We start by simplifying $T_{j}$ using Taylor
series approximations. Recall the Taylor expansion

$\frac{1}{x+h}=\frac{1}{x}-\frac{h}{x^{2}}+\mathcal{O}(h^{2}).$

We apply Taylor this expansion to $T_{j}$ with $x=4$ and
$h=l^{2}\left(\beta^{2}\left(\nabla\psi_{j}\right)^{2}-\beta\Delta\psi_{j}\right)$:

$\displaystyle T_{j}$
$\displaystyle=\frac{e^{\beta\psi_{j}}}{4}\left(1-\frac{l^{2}}{4}\left(\beta^{2}\left(\nabla\psi_{j}\right)^{2}-\beta\Delta\psi_{j}\right)\right)+\mathcal{O}(l^{4}).$
(17)

Note that here, we are depending on the smoothness of $\psi_{j}$. Then, by
taking the gradient of (17), we find that

$\displaystyle\nabla T_{j}$
$\displaystyle=\frac{e^{\beta\psi_{j}}}{4}\left(\beta\nabla\psi_{j}-\frac{l^{2}}{4}\left(\beta^{3}(\nabla\psi_{j})^{3}+\beta^{2}\nabla\psi_{j}\Delta\psi_{j}-\beta\nabla^{3}\psi_{j}\right)\right)+\mathcal{O}(l^{4}).$
(18)

We also can find $\Delta T_{j}$: $\displaystyle\Delta T_{j}$
$\displaystyle=\nabla\cdot(\nabla T_{j})$ $\displaystyle\Rightarrow\Delta
T_{j}$
$\displaystyle=\frac{e^{\beta\psi_{j}}}{4}\bigg{(}\left(\beta^{2}(\nabla\psi_{j})^{2}+\beta\Delta\psi_{j}\right)-\frac{l^{2}}{4}\Big{(}4\beta^{3}(\nabla\psi_{j})^{2}\Delta\psi_{j}$
$\displaystyle+\beta^{2}(\Delta\psi_{j})^{2}+\beta^{4}(\nabla\psi_{j})^{4}-\beta\nabla^{4}\psi_{j}\Big{)}\bigg{)}+\mathcal{O}(l^{4}).$

We now focus our attention on the movement probability. Starting from
definition (3), recall that the probability that an agent from gang $j$ moves
from site $(x,y)$ to a neighboring site $(x_{1},y_{1})$ is

$M_{j}(x\rightarrow x_{1},y\rightarrow
y_{1},t)=\frac{e^{-\beta\psi_{j}(x_{1},y_{1},t)}}{\sum\limits_{(\tilde{x},\tilde{y})\sim(x,y)}e^{-\beta\psi_{j}(\tilde{x},\tilde{y},t)}},$

where $(\tilde{x},\tilde{y})$ are the neighbors of site $(x,y)$. We now
slightly modify the above definition so that we evaluate the probability that
an agent at a neighboring site $(\tilde{x},\tilde{y})$ moves to site $(x,y)$:

$M_{j}(\tilde{x}\rightarrow x,\tilde{y}\rightarrow
y,t)=\frac{e^{-\beta\psi_{j}(x,y,t)}}{\sum\limits_{(\tilde{\tilde{x}},\tilde{\tilde{y}})\sim(\tilde{x},\tilde{y})}e^{-\beta\psi_{j}(\tilde{\tilde{x}},\tilde{\tilde{y}},t)}},$
(21)

where $(\tilde{\tilde{x}},\tilde{\tilde{y}})$ are the four neighbors of site
$(\tilde{x},\tilde{y})$.

To remove the presence of the neighbors’ neighbors
$(\tilde{\tilde{x}},\tilde{\tilde{y}})$ from the denominator, we apply the
discrete Laplacian to find thats

$\sum\limits_{(\tilde{\tilde{x}},\tilde{\tilde{y}})\sim(\tilde{x},\tilde{y})}e^{-\beta\psi_{j}(\tilde{\tilde{x}},\tilde{\tilde{y}},t)}=4e^{-\beta\psi_{j}(\tilde{x},\tilde{y},t)}+l^{2}\Delta\left(e^{-\beta\psi_{j}(\tilde{x},\tilde{y},t)}\right)+\mathcal{O}(l^{4}).$
(22)

Noting that

$\displaystyle\Delta e^{-\beta\psi_{j}(\tilde{x},\tilde{y},t)}$
$\displaystyle=\nabla\cdot\nabla\left(e^{-\beta\psi_{j}(\tilde{x},\tilde{y},t)}\right)$
$\displaystyle=\nabla\cdot\left(-\beta\nabla\psi_{j}(\tilde{x},\tilde{y},t)e^{-\beta\psi_{j}(\tilde{x},\tilde{y},t)}\right)$
$\displaystyle=\left[\beta^{2}\left(\nabla\psi_{j}(\tilde{x},\tilde{y},t)\right)^{2}-\beta\Delta\psi_{j}(\tilde{x},\tilde{y},t)\right]e^{-\beta\psi_{j}(\tilde{x},\tilde{y},t)},$
(23)

Combining equations (22) and (23) gives us

$\displaystyle\sum\limits_{(\tilde{\tilde{x}},\tilde{\tilde{y}})\sim(\tilde{x},\tilde{y})}e^{-\beta\psi_{j}(\tilde{\tilde{x}},\tilde{\tilde{y}},t)}=$
$\displaystyle
e^{-\beta\psi_{j}(\tilde{x},\tilde{y},t)}\left(4+l^{2}\left(\beta^{2}\left(\nabla\psi_{j}(\tilde{x},\tilde{y},t)\right)^{2}-\beta\Delta\psi_{j}(\tilde{x},\tilde{y},t)\right)\right)$
$\displaystyle+\mathcal{O}(l^{4}).$

Substituting it back into equation (21), and replacing the denominator gives
us

$\displaystyle M_{j}(\tilde{x}\rightarrow x,\tilde{y}\rightarrow y,t)$
$\displaystyle=\frac{e^{-\beta\psi_{j}(x,y,t)}}{\left[4+l^{2}\left(\beta^{2}\left(\nabla\psi_{j}(\tilde{x},\tilde{y},t)\right)^{2}-\beta\Delta\psi_{j}(\tilde{x},\tilde{y},t)\right)\right]e^{-\beta\psi_{j}(\tilde{x},\tilde{y},t)}+\mathcal{O}(l^{4})}$
$\displaystyle\approx
e^{-\beta\psi_{j}(x,y,t)}\left[\frac{e^{\beta\psi_{j}(\tilde{x},\tilde{y},t)}}{4+l^{2}\left(\beta^{2}\left(\nabla\psi_{j}(\tilde{x},\tilde{y},t)\right)^{2}-\beta\Delta\psi_{j}(\tilde{x},\tilde{y},t)\right)}\right].$

The term inside the large brackets in the equation above takes the form of
(16) where $(x,y,t)$ is replaced with $(\tilde{x},\tilde{y},t)$, yielding the
following approximation:

$M_{j}(\tilde{x}\rightarrow x,\tilde{y}\rightarrow y,t)\approx
e^{-\beta\psi_{j}(x,y,t)}T_{j}(\tilde{x},\tilde{y},t).$ (24)

####  4.2.2 The Derivation

We now have all the tools needed to formally derive the agent density
continuum equations. We will be using the discrete Laplacian approximation in
order to approximate the influence of the neighbors of site $(x,y)$. We will
also be using equation (24) to simplify the discrete model.

Starting from the discrete model, we recall equation (5):

$\begin{split}\rho_{j}(x,y,t+\delta
t)=&\rho_{j}(x,y,t)+\sum_{(\tilde{x},\tilde{y})\sim(x,y)}\rho_{j}(\tilde{x},\tilde{y},t)M_{j}(\tilde{x}\rightarrow
x,\tilde{y}\rightarrow y,t)\\\
&-\rho_{j}(x,y,t)\sum_{(\tilde{x},\tilde{y})\sim(x,y)}M_{j}(x\rightarrow\tilde{x},y\rightarrow\tilde{y},t).\end{split}$

Rearranging the equation and dividing both sides by $\delta t$ gives us

$\begin{split}\frac{\rho_{j}(x,y,t+\delta t)-\rho_{j}(x,y,t)}{\delta
t}=\frac{1}{\delta
t}\left[\sum_{(\tilde{x},\tilde{y})\sim(x,y)}\rho_{j}(\tilde{x},\tilde{y},t)M_{j}(\tilde{x}\rightarrow
x,\tilde{y}\rightarrow y,t)\right.\\\
\left.-\rho_{j}(x,y,t)\sum_{(\tilde{x},\tilde{y})\sim(x,y)}M_{j}(x\rightarrow\tilde{x},y\rightarrow\tilde{y},t)\right].\end{split}$

By equation (24), and noting that each agent has to move to one of the
neighboring sites,

$\displaystyle\frac{\rho_{j}(x,y,t+\delta t)-\rho_{j}(x,y,t)}{\delta
t}=\frac{1}{\delta t}\Bigg{[}$ $\displaystyle
e^{-\beta\psi_{j}(x,y,t)}\sum_{(\tilde{x},\tilde{y})\sim(x,y)}\rho_{j}(\tilde{x},\tilde{y},t)T_{j}(\tilde{x},\tilde{y},t)$
$\displaystyle-$ $\displaystyle\rho_{j}(x,y,t)+\mathcal{O}(l^{4})\Bigg{]}.$
(25)

Using the discrete Laplacian technique, we can approximate the contribution of
the neighboring sites, giving on the right-hand side

$\displaystyle\frac{1}{\delta t}\Bigg{[}$ $\displaystyle
e^{-\beta\psi_{j}(x,y,t)}\bigg{(}4\rho_{j}(x,y,t)T_{j}(x,y,t)+l^{2}\Delta\Big{(}\rho_{j}(x,y,t)T_{j}(x,y,t)\Big{)}\bigg{)}$
$\displaystyle-\rho_{j}(x,y,t)+\mathcal{O}(l^{4})\Bigg{]}.$

The notation $(x,y,t)$ is again dropped as there are no longer any neighbors
$(\tilde{x},\tilde{y},t)$ remaining in this derivation. Hence, we write,

$\frac{1}{\delta
t}\Bigg{[}e^{-\beta\psi_{j}}\bigg{(}4\rho_{j}T_{j}+l^{2}\Delta\Big{(}\rho_{j}T_{j})\Big{)}\bigg{)}-\rho_{j}+\mathcal{O}(l^{4})\Bigg{]}.$
(26)

From definition (16), $T_{j}(x,y,t)$ is substituted back into the first term
of (26):

$\displaystyle\frac{1}{\delta
t}\Bigg{[}4\rho_{j}e^{-\beta\psi_{j}}\left(\frac{e^{\beta\psi_{j}}}{4+l^{2}\left(\left(\beta\nabla\psi_{j}\right)^{2}-\beta\Delta\psi_{j}\right)}\right)+l^{2}e^{-\beta\psi_{j}}\Delta\Big{(}\rho_{j}T_{j}\Big{)}-\rho_{j}+\mathcal{O}(l^{4})\Bigg{]}.$
Simplifying the expression yields $\displaystyle\frac{1}{\delta
t}\Bigg{[}4\rho_{j}\left(\frac{1}{4+l^{2}\left(\left(\beta\nabla\psi_{j}\right)^{2}-\beta\Delta\psi_{j}\right)}\right)-\rho_{j}+l^{2}e^{-\beta\psi_{j}}\Delta\Big{(}\rho_{j}T_{j}\Big{)}+\mathcal{O}(l^{4})\Bigg{]}.$
(27)

Using a Taylor series expansion on the first term within the brackets yields,

$\left(\frac{1}{4+l^{2}\left(\left(\beta\nabla\psi_{j}\right)^{2}-\beta\Delta\psi_{j}\right)}\right)=\frac{1}{4}-\frac{l^{2}\left(\left(\beta\nabla\psi_{j}\right)^{2}-\beta\Delta\psi_{j}\right)}{4^{2}}+\mathcal{O}(l^{4}).$

Expression (27) thus becomes,

$\frac{1}{\delta
t}\left[4\rho_{j}\left(\frac{1}{4}-\frac{l^{2}\left(\left(\beta\nabla\psi_{j}\right)^{2}-\beta\Delta\psi_{j}\right)}{4^{2}}\right)-\rho_{j}+l^{2}e^{-\beta\psi_{j}}\Delta\Big{(}\rho_{j}T_{j}\Big{)}+\mathcal{O}(l^{4})\right].$

Simplifying the expression yields

$\begin{split}\frac{\rho_{j}(x,y,t+\delta t)-\rho_{j}(x,y,t)}{\delta
t}=&\frac{l^{2}}{\delta
t}\left[-\frac{\rho_{j}}{4}\left(\left(\beta\nabla\psi_{j}\right)^{2}-\beta\Delta\psi_{j}\right)+e^{-\beta\psi_{j}}\Delta\Big{(}\rho_{j}T_{j}\Big{)}\right]\\\
&+\mathcal{O}\left(\frac{l^{4}}{\delta t}\right).\end{split}$ (28)

However, we can further simplify this by noting that

$\Delta\Big{(}\rho_{j}T_{j}\Big{)}=\Big{(}T_{j}\Delta\rho_{j}+2\nabla
T_{j}\nabla\rho_{j}+\rho_{j}\Delta T_{j}\Big{)}.$

From (17) through (4.2.1), we have

$\displaystyle T_{j}$
$\displaystyle=\frac{e^{\beta\psi_{j}}}{4}+\mathcal{O}(l^{2}),$
$\displaystyle\nabla T_{j}$ $\displaystyle=\frac{\beta
e^{\beta\psi_{j}}}{4}\nabla\psi_{j}+\mathcal{O}(l^{2}),$ $\displaystyle\Delta
T_{j}$
$\displaystyle=\frac{e^{\beta\psi_{j}}}{4}\left(\beta\Delta\psi_{j}+(\beta\nabla\psi_{j})^{2}\right)+\mathcal{O}(l^{2}).$

Therefore,

$\displaystyle\Delta\Big{(}\rho_{j}T_{j}\Big{)}$
$\displaystyle=\frac{e^{\beta\psi_{j}}}{4}\Delta\rho_{j}+\frac{2\beta
e^{\beta\psi_{j}}}{4}\nabla\psi_{j}\nabla\rho_{j}+\frac{e^{\beta\psi_{j}}}{4}\rho_{j}\left(\beta\Delta\psi_{j}+(\beta\nabla\psi_{j})^{2}\right)+\mathcal{O}(l^{2})$
$\displaystyle=\frac{e^{\beta\psi_{j}}}{4}\left[\Delta\rho_{j}+2\beta\nabla\psi_{j}\nabla\rho_{j}+\rho_{j}\left(\left(\beta\nabla\psi_{j}\right)^{2}+\beta\Delta\psi_{j}\right)\right]+\mathcal{O}(l^{2}).$
(29)

Substituting (29) back into (28) gives us

$\displaystyle\frac{\rho_{j}(x,y,t+\delta t)-\rho_{j}(x,y,t)}{\delta t}$
$\displaystyle=\frac{l^{2}}{4\delta
t}\Bigg{[}-\rho_{j}\left((\beta\nabla\psi_{j})^{2}-\beta\Delta\psi_{j}\right)+\Delta\rho_{j}+2\beta\nabla\psi_{j}\nabla\rho_{j}$
$\displaystyle\qquad+\rho_{j}\left(\left(\beta\nabla\psi_{j}\right)^{2}+\beta\Delta\psi_{j}\right)\Bigg{]}+\mathcal{O}\left(\frac{l^{4}}{\delta
t}\right)$ $\displaystyle=\frac{l^{2}}{4\delta
t}\Bigg{[}\Delta\rho_{j}+2\beta\nabla\rho_{j}\nabla\psi_{j}+2\beta\rho_{j}\Delta\psi_{j}\bigg{]}+\mathcal{O}\left(\frac{l^{4}}{\delta
t}\right)$ $\displaystyle=\frac{l^{2}}{4\delta
t}\Bigg{[}\Delta\rho_{j}+2\beta\left(\nabla\rho_{j}\nabla\psi_{j}+\rho_{j}\Delta\psi_{j}\right)\bigg{]}+\mathcal{O}\left(\frac{l^{4}}{\delta
t}\right)$ $\displaystyle=\frac{l^{2}}{4\delta
t}\Bigg{[}\Delta\rho_{j}+2\beta\nabla\cdot\Big{(}\rho_{j}\nabla\psi_{j}\Big{)}\Bigg{]}+\mathcal{O}\left(\frac{l^{4}}{\delta
t}\right)$ $\displaystyle=\frac{l^{2}}{4\delta
t}\nabla\cdot\Bigg{[}\nabla\rho_{j}+2\beta\Big{(}\rho_{j}\nabla\psi_{j}\Big{)}\Bigg{]}+\mathcal{O}\left(\frac{l^{4}}{\delta
t}\right).$

Assuming that the agent density $\rho_{j}$ is sufficiently smooth and the
following limits

$\begin{split}l&\rightarrow 0,\\\ \delta t&\rightarrow 0,\\\
\frac{l^{2}}{\delta t}&\rightarrow D,\end{split}$ (30)

gives us the final form for the continuum equations for the density of gang
$j$ agents:

$\frac{\partial\rho_{j}}{\partial
t}=\frac{D}{4}\nabla\cdot\Bigg{[}\nabla\rho_{j}+2\beta\Big{(}\rho_{j}\nabla\psi_{j}\Big{)}\Bigg{]}.$
(31)

Finally, from (15) and (31), and using equation (2) to express everything in
terms of agent density and graffiti density, the limiting convection-diffusion
system for our model is

$\begin{cases}\displaystyle\frac{\partial\xi_{j}}{\partial
t}(x,y,t)=\gamma\rho_{j}(x,y,t)-\lambda\xi_{j}(x,y,t)\\\
\displaystyle\frac{\partial\rho_{j}}{\partial
t}(x,y,t)=\frac{D}{4}\nabla\cdot\left[\nabla\rho_{j}(x,y,t)+2\beta\left(\rho_{j}(x,y,t)\nabla\left(\sum_{\begin{subarray}{c}i=1\\\
i\neq j\end{subarray}}^{K}\xi_{i}(x,y,t)\right)\right)\right]\end{cases}$ (32)

for $j=1,2,\dots K$ with periodic boundary conditions.

###  4.3 Steady-State Solutions

Considering steady-state solutions for the graffiti density, we find from the
evolution equations for for the graffiti density that

$\displaystyle\frac{\partial\xi_{j}}{\partial t}(x,y,t)$ $\displaystyle=0$
$\displaystyle\Rightarrow\gamma\rho_{j}(x,y,t)-\lambda\xi_{j}(x,y,t)$
$\displaystyle=0$ $\displaystyle\Rightarrow\xi_{j}$
$\displaystyle=\frac{\gamma}{\lambda}\rho_{j}.$ (33)

We now focus our attention on the steady-state solutions for the agent density
of gang $j$:

$\displaystyle\frac{\partial\rho_{j}}{\partial
t}(x,y,t)=\frac{D}{4}\nabla\cdot\left[\nabla\rho_{j}(x,y,t)+2\beta\left(\rho_{j}(x,y,t)\nabla\left(\sum_{\begin{subarray}{c}i=1\\\
i\neq j\end{subarray}}^{K}\xi_{i}(x,y,t)\right)\right)\right]$
$\displaystyle=0,$

considering solutions of the form

$\displaystyle\nabla\rho_{j}(x,y,t)+2\beta\left(\rho_{j}(x,y,t)\nabla\left(\sum_{\begin{subarray}{c}i=1\\\
i\neq j\end{subarray}}^{K}\xi_{i}(x,y,t)\right)\right)$ $\displaystyle=c.$

Using the steady-state graffiti density derived in equation (33), we find that

$\displaystyle\nabla\rho_{j}(x,y,t)+\frac{2\beta\gamma^{2}}{\lambda^{2}}\left(\rho_{j}(x,y,t)\nabla\cdot\left(\sum_{\begin{subarray}{c}i=1\\\
i\neq j\end{subarray}}^{K}\rho_{i}(x,y,t)\right)\right)$ $\displaystyle=c.$
(34)

Any form of $\rho_{j}(x,y,t)$ satisfying the above equation is a steady-state
solution of our system. For simplicity, we suppose that $\rho_{j}$ is a
constant for all $j$. In this case, the steady-state solution of our problem
takes the form:

$\displaystyle\begin{cases}\xi_{j}&=\frac{\gamma}{\lambda}\rho_{j},\\\
\rho_{j}&=c_{j},\end{cases}$ (35)

for $j=1,2,\dots,K$ and $c_{j}$ is a positive constant.

To test whether these steady-state solutions of the continuum system (35)
agree with the discrete model, we first start our simulations in a steady-
state solution from (35). Here we will start the simulations with the agents
from all three gangs are completely segregated. The results of the simulation
are visualized in Figure 5. In that figure, we clearly see that the agents
remain segregated and the system does not deviate from the steady-state
solutions over time.

Figure 5: Temporal evolution of a steady-state solution for the agent density.
Here we have $N_{1}=N_{2}=N_{3}=50,625$, with $\lambda=\gamma=0.5$,
$\beta=3\times 10^{-4}$ and the lattice size is $75\times 75$. It is shown
that if the system started initially at a steady-state then it would remain
there.

From the steady-state solutions of the continuum system in (35) we see that
the graffiti density is $\xi=\frac{\gamma}{\lambda}\rho$, which implies that
the steady-state solution of the graffiti and agent densities at a site are
equal up to scaled amount of $\frac{\gamma}{\lambda}.$ To check whether this
generally holds for our discrete model in Section 3, we test our discrete
system with several ratio values that are different. We also take a cross-
sectional slice over the lattice at the first and final time steps. This would
allow us to see if the steady-state of the discrete model and the continuum
system agree. We visualize the results in Figure 6.

Figure 6: Cross-sectional slices of the agent and graffiti densities for
different $\frac{\gamma}{\lambda}$ at the final time step for a segregated
state. Here we have $N_{1}=N_{2}=N=3=50,000$ with $\delta t=1$, $\beta=3\times
10^{-4}$ and the lattice size is $100\times 100$. In the top row we have
$\gamma=0.25$ and $\lambda=0.5$, for the middle row $\gamma=0.5$ and
$\lambda=0.5$ and $\gamma=0.5$ and $\lambda=0.25$ for the bottom row. From
looking at the cross-sectional slices we clearly that
$\xi_{j}\approx\frac{\gamma}{\lambda}\rho_{j},j=1,2,3$ and this agrees with
our steady-state solution which was $\xi_{j}=\frac{\gamma}{\lambda}\rho_{j}$.

##  5 Linear Stability Analysis

To have a better understanding of our system we linearize our model by
considering a perturbation of the equilibrium solution (35) to the well-mixed
state. We assume that our perturbations are of the form $\epsilon=\delta
e^{\alpha t}e^{ikx}$, with $\delta<<1$, and in that case our solution will
take the following form:

$\displaystyle\begin{cases}\xi_{j}&=\bar{\xi_{j}}+\delta_{\xi_{j}}e^{\alpha
t}e^{ikx}\\\ \rho_{j}&=\bar{\rho_{j}}+\delta_{\rho_{j}}e^{\alpha
t}e^{ikx}\end{cases}$ (36)

Here, $e^{ikx}=\cos(kx)+i\sin(kx)$, where $k$ represents the wave number of
the spatial wave. In order for the equilibrium solution to be stable, $\alpha$
must be negative so that it forces the perturbations to decay as time
increases. For more examples of this kind of perturbation being used to study
the stability of equilibrium solutions, the interested reader is referred to
[30, 18, 14, 31].

To analyze the dynamics of these solutions, we now substitute (36) into the
evolution equations (32). We start with the first equation:

$\frac{\partial\xi_{j}}{\partial t}=\gamma\rho_{j}-\lambda\xi_{j}.$ (37)

Substituting (36) into (37) yields

$\displaystyle\frac{\partial}{\partial
t}\left(\bar{\xi_{j}}+\delta_{\xi_{j}}e^{\alpha t}e^{ikx}\right)$
$\displaystyle=\gamma(\bar{\rho_{j}}+\delta_{\rho_{j}}e^{\alpha
t}e^{ikx})-\lambda(\bar{\xi_{j}}+\delta_{\xi_{j}}e^{\alpha t}e^{ikx}).$ Since
we assumed $\bar{\xi_{j}}$ to be an equilibrium solution, its derivative with
respect to time is zero, $\displaystyle\alpha\delta_{\xi_{j}}e^{\alpha
t}e^{ikx}$
$\displaystyle=(\gamma\bar{\rho_{j}}-\lambda\bar{\xi_{j}})+(\gamma\delta_{\rho_{j}}-\lambda\delta_{\xi_{j}})e^{\alpha
t}e^{ikx},$
$\displaystyle=(\gamma\delta_{\rho_{j}}-\lambda\delta_{\xi_{j}})e^{\alpha
t}e^{ikx}.$ Since
$\gamma\bar{\rho_{j}}-\lambda\bar{\xi_{j}}=\frac{\partial\bar{\xi_{j}}}{\partial
t}=0$. Hence $\displaystyle\alpha\delta_{\xi_{j}}$
$\displaystyle=(\gamma\delta_{\rho_{j}}-\lambda\delta_{\xi_{j}}),\text{ for
j=1,2, \ldots, K.}$ (38)

Next, we substitute (36) into the evolution equation

$\frac{\partial\rho_{j}}{\partial
t}=\frac{D}{4}\nabla\cdot\Big{[}\nabla\rho_{j}+2\beta\big{(}\rho_{j}\nabla\psi_{j}\big{)}\Big{]}$

giving us

$\displaystyle\frac{\partial}{\partial
t}\left(\bar{\rho_{j}}+\delta_{\rho_{j}}e^{\alpha t}e^{ikx}\right)$
$\displaystyle=\frac{D}{4}\Delta\left(\bar{\rho_{j}}+\delta_{\rho_{j}}e^{\alpha
t}e^{ikx}\right)$
$\displaystyle\quad+\frac{D\beta}{2}\nabla\cdot\left((\bar{\rho_{j}}+\delta_{\rho_{j}}e^{\alpha
t}e^{ikx})\nabla\left(\sum_{\begin{subarray}{c}l=1\\\ l\neq
j\end{subarray}}^{K}(\bar{\xi_{l}}+\delta_{\xi_{l}}e^{\alpha
t}e^{ikx})\right)\right).$

Since we are working in one dimension and our equilibrium solution is a
constant in both space and time,

$\displaystyle\alpha\delta_{\rho_{j}}e^{\alpha t}e^{ikx}$
$\displaystyle=\frac{-D|k|^{2}}{4}\delta_{\rho_{j}}e^{\alpha
t}e^{ikx}+\frac{D\beta}{2}\frac{d}{dx}\left((\bar{\rho_{j}}+\delta_{\rho_{j}}e^{\alpha
t}e^{ikx})(ik\sum_{\begin{subarray}{c}l=1\\\ l\neq
j\end{subarray}}^{K}\delta_{\xi_{l}}e^{\alpha t}e^{ikx})\right)$
$\displaystyle=\frac{-D|k|^{2}}{4}\delta_{\rho_{j}}e^{\alpha
t}e^{ikx}+\frac{D\beta}{2}\frac{d}{dx}\left(ik\bar{\rho_{j}}\sum_{\begin{subarray}{c}l=1\\\
l\neq j\end{subarray}}^{K}\delta_{\xi_{l}}e^{\alpha
t}e^{ikx}\right)+\mathcal{O}(\delta_{\rho_{j}}\sum_{\begin{subarray}{c}l=1\\\
l\neq j\end{subarray}}^{K}\delta_{\xi_{l}})$
$\displaystyle=\frac{-D|k|^{2}}{4}\delta_{\rho_{j}}e^{\alpha
t}e^{ikx}-\frac{D\beta|k|^{2}}{2}\bar{\rho_{j}}\sum_{\begin{subarray}{c}l=1\\\
l\neq j\end{subarray}}^{K}\delta_{\xi_{l}}e^{\alpha
t}e^{ikx}+\mathcal{O}(\delta_{\rho_{j}}\sum_{\begin{subarray}{c}l=1\\\ l\neq
j\end{subarray}}^{K}\delta_{\xi_{2}})$
$\displaystyle=\frac{-D|k|^{2}}{4}\left(\delta_{\rho_{j}}+2\beta\bar{\rho_{j}}\sum_{\begin{subarray}{c}l=1\\\
l\neq j\end{subarray}}^{K}\delta_{\xi_{l}}\right)e^{\alpha
t}e^{ikx}+\mathcal{O}(\delta_{\rho_{j}}\sum_{\begin{subarray}{c}l=1\\\ l\neq
j\end{subarray}}^{K}\delta_{\xi_{l}}).$

We can safely neglect the term
$\displaystyle\mathcal{O}(\delta_{\rho_{j}}\sum_{\begin{subarray}{c}l=1\\\
l\neq j\end{subarray}}^{K}\delta_{\xi_{l}})$, since
$|\delta_{\rho_{j}}|,|\delta_{\xi_{l}}|<<1$; therefore,

$\alpha\delta_{\rho_{j}}=\frac{-D|k|^{2}}{4}\left(\delta_{\rho_{j}}+2\beta\bar{\rho_{j}}\sum_{\begin{subarray}{c}l=1\\\
l\neq j\end{subarray}}^{K}\delta_{\xi_{l}}\right),\text{ for j=1, 2, $\dots$,
K.}$ (39)

Next, we write the equations from (38) and (39) in a systems form:

$\displaystyle(\gamma\delta_{\rho_{j}}-\lambda\delta_{\xi_{j}})$
$\displaystyle=\alpha\delta_{\xi_{j}}$
$\displaystyle\frac{-D|k|^{2}}{4}\left(\delta_{\rho_{j}}+2\beta\bar{\rho_{j}}\sum_{\begin{subarray}{c}l=1\\\
l\neq j\end{subarray}}^{K}\delta_{\xi_{l}}\right)$
$\displaystyle=\alpha\delta_{\rho_{j}},\text{ where }j=1,\dots,K.$

For simplicity we consider the case where $K=3$, with all gangs having the
same $\beta$ parameter and write the system in matrix vector format, giving
us:

$\begin{bmatrix}-\lambda&0&0&\gamma&0&0\\\ 0&-\lambda&0&0&\gamma&0\\\
0&0&-\lambda&0&0&\gamma\\\ 0&\frac{-\beta
D\bar{\rho_{1}}|k|^{2}}{2}&\frac{-\beta
D\bar{\rho_{1}}|k|^{2}}{2}&\frac{-D|k|^{2}}{4}&0&0\\\ \frac{-\beta
D\bar{\rho_{2}}|k|^{2}}{2}&0&\frac{-\beta
D\bar{\rho_{2}}|k|^{2}}{2}&0&\frac{-D|k|^{2}}{4}&0\\\ \frac{-\beta
D\bar{\rho_{3}}|k|^{2}}{2}&\frac{-\beta
D\bar{\rho_{3}}|k|^{2}}{2}&0&0&0&\frac{-D|k|^{2}}{4}\\\
\end{bmatrix}\begin{bmatrix}\delta_{\xi_{1}}\\\ \delta_{\xi_{2}}\\\
\delta_{\xi_{3}}\\\ \delta_{\rho_{1}}\\\ \delta_{\rho_{2}}\\\
\delta_{\rho_{3}}\end{bmatrix}=\alpha\begin{bmatrix}\delta_{\xi_{1}}\\\
\delta_{\xi_{2}}\\\ \delta_{\xi_{3}}\\\ \delta_{\rho_{1}}\\\
\delta_{\rho_{2}}\\\ \delta_{\rho_{3}}\end{bmatrix}.$

This gives us

$\displaystyle F\vec{\delta}$ $\displaystyle=\alpha\vec{\delta}$
$\displaystyle\iff\left(F-\alpha I_{4}\right)\vec{\delta}$ $\displaystyle=0,$

which reduces to an eigenvalue problem for matrix $F$. For the problem to have
a non trivial solution (i.e. $\vec{\delta}\neq 0$), the determinant of
$(F-\alpha I_{4})$ must be zero. Therefore,

$\begin{vmatrix}-(\lambda+\alpha)&0&0&\gamma&0&0\\\
0&-(\lambda+\alpha)&0&0&\gamma&0\\\ 0&0&-(\lambda+\alpha)&0&0&\gamma\\\
0&\frac{-\beta D\bar{\rho_{1}}|k|^{2}}{2}&\frac{-\beta
D\bar{\rho_{1}}|k|^{2}}{2}&-\left(\frac{D|k|^{2}}{4}+\alpha\right)&0&0\\\
\frac{-\beta D\bar{\rho_{2}}|k|^{2}}{2}&0&\frac{-\beta
D\bar{\rho_{2}}|k|^{2}}{2}&0&-\left(\frac{D|k|^{2}}{4}+\alpha\right)&0\\\
\frac{-\beta D\bar{\rho_{3}}|k|^{2}}{2}&\frac{-\beta
D\bar{\rho_{3}}|k|^{2}}{2}&0&0&0&-\left(\frac{D|k|^{2}}{4}+\alpha\right)\\\
\end{vmatrix}=0,$

giving us the following characteristic polynomial

$\displaystyle f(\alpha)$ $\displaystyle=\frac{1}{64}\Big{[}4\alpha
D^{2}(\alpha+\lambda)(3\alpha^{2}+6\alpha\lambda+3\lambda^{2}-4\beta^{2}\gamma^{2}(\bar{\rho_{1}}\bar{\rho_{2}}+\bar{\rho_{1}}\bar{\rho_{3}}+\bar{\rho_{2}}\bar{\rho_{3}}))k^{4}$
$\displaystyle+D^{3}\Big{(}\alpha^{3}+3\alpha^{2}\lambda+\lambda^{3}+16\beta^{3}\gamma^{3}\bar{\rho_{1}}\bar{\rho_{2}}\bar{\rho_{3}}-4\beta^{2}\gamma^{2}\lambda(\bar{\rho_{1}}\bar{\rho_{2}}+\bar{\rho_{1}}\bar{\rho_{3}}+\bar{\rho_{2}}\bar{\rho_{3}})$
$\displaystyle+\alpha(3\lambda^{2}-4\beta^{2}\gamma^{2}(\bar{\rho_{1}}\bar{\rho_{2}}+\bar{\rho_{1}}\bar{\rho_{3}}+\bar{\rho_{2}}\bar{\rho_{3}}))\Big{)}k^{6}+64\alpha^{3}(\alpha+\lambda)^{3}+48\alpha^{2}D(\alpha+\lambda)^{3}k^{2}\Big{]}$
$\displaystyle=0.$

Solving the characteristic polynomial gives six eigenvalues, which we solve
numerically using Mathematica and plot the results for different values of
$\beta$ in Figure 7.

Figure 7: The six eigenvalues versus the wave number k plotted for different
$\beta$ values. Here we have $D=0.0001$, mass $=150,000$ and the
$\frac{\gamma}{\lambda}$ ratio $=1$.

To determine the stability of our system, we recall that the system becomes
linearly unstable when any eigenvalue has a positive real part. We see in the
top two rows of Figure 7 that for small $\beta$ values, none of the
eigenvalues had a positive real part, and thus the system remained stable. In
terms of our model, this makes sense since in our discrete simulations, the
system remains well-mixed for these $\beta$ values. However, when we increase
the value of $\beta$ to $0.00005$, we see in the bottom three figures that the
second and sixth eigenvalues have a positive real part, and the system has
thus become linearly unstable. This again agrees with our physical intuition.
We note that the values $\beta$ from the discrete model matches those of the
linearized system of partial differential equations.

##  6 Variations of the model: varying $\beta$ by gang

In Section 2.1 equation (3), we defined the probability that an agent from
gang $j$ moves from site $s_{1}=(x_{1},y_{1})\in S$ to one of the neighboring
sites $s_{2}=(x_{2},y_{2})\in S$ to be

$M_{j}(x_{1}\rightarrow x_{2},y_{1}\rightarrow
y_{2},t)=\frac{e^{-\beta\psi_{j}(x_{2},y_{2},t)}}{\sum\limits_{(\tilde{x},\tilde{y})\sim(x_{1},y_{1})}e^{-\beta\psi_{j}(\tilde{x},\tilde{y},t)}},$

with

$\psi_{j}(x,y,t):=\sum_{\begin{subarray}{c}i=1\\\ i\neq
j\end{subarray}}^{K}\xi_{i}(x,y,t)$

from equation (2). The parameter $\beta$ then controls how strongly each gang
reacts to the graffiti of the other gangs. However, it is reasonable to
consider that this parameter $\beta$ might vary by gang. Here, we explore
variations of the model incorporating this idea. In this section, we will make
two different modifications of (3) and explore how these modifications affect
the system of PDEs and the segregation behavior of the model.

###  6.1 Timidity Model (Variation 1)

In the first modification of the model, instead of having identical $\beta$
values for all gangs, we change it so that gang $j$ has a distinct
corresponding $\beta$ value, denoted by $\beta_{j}$. This $\beta_{j}$
determines how much attention gang $j$ places on the graffiti of the other
gangs. In essence, this $\beta_{j}$ encodes the timidity of gang $j$, with
higher $\beta_{j}$ corresponding to higher timidity, causing gang $j$ to more
strongly avoid other gangs’ graffiti. Hence, the modified definition for
movement becomes,

$M_{j}(x_{1}\rightarrow x_{2},y_{1}\rightarrow
y_{2},t)=\frac{e^{-\beta_{j}\psi_{j}(x_{2},y_{2},t)}}{\sum\limits_{(\tilde{x},\tilde{y})\sim(x_{1},y_{1})}e^{-\beta_{j}\psi_{j}(\tilde{x},\tilde{y},t)}}.$
(40)

In this variation of the model, gang $j$ avoids all other gangs’ graffiti with
rate $\beta_{j}$. All of the graffiti from other gangs count equally and are
identically avoided. For example, let us consider the case of three gangs $1$,
$2$, and $3$ such that gang $2$ has a relatively large $\beta_{2}$ value, gang
$3$ has a relatively small $\beta_{3}$ value, and gang $1$ has an intermediate
$\beta_{1}$ value. Then gang $2$’s agents would strongly avoid areas where the
other two gangs, $1$ and $3$, have tagged. Gang $3$’s agents, on the other
hand, would more freely on the lattice, as the small $\beta_{3}$ value leads
it to not place much importance on other gangs’ graffiti. Gang $1$’s agents’
movement dynamics would lie somewhere in between.

If one follows the derivation of the continuum equations in Section 4 but
replacing (3) with (40), it can be easily shown that the resulting system of
equations for $j=1,2,\dots,K$ are

$\begin{cases}\displaystyle\frac{\partial\xi_{j}}{\partial
t}(x,y,t)=\gamma\rho_{j}(x,y,t)-\lambda\xi_{j}(x,y,t)\\\
\displaystyle\frac{\partial\rho_{j}}{\partial
t}(x,y,t)=\frac{D}{4}\nabla\cdot\left[\nabla\rho_{j}(x,y,t)+2\beta_{j}\left(\rho_{j}(x,y,t)\nabla\left(\sum_{\begin{subarray}{c}i=1\\\
i\neq j\end{subarray}}^{K}\xi_{i}(x,y,t)\right)\right)\right]\end{cases}$ (41)

with periodic boundary conditions. We can see that the $\beta_{j}$ values will
then affect the balance between the diffusion and the advection terms
differently depending on the gang affiliation, making diffusion relatively
stronger for those gangs with lower $\beta_{j}$ values.

To test how these changes affect our discrete model, we ran our simulations
with three gangs $1,2$ and $3$; all gangs are assumed to have the identical
number of agents $N=50,000$. We also assume that the lattice size $L\times L$
is equal to $100\times 100$, and will use $100,000$ time steps with each step
size $\delta t=1$. We assigned $\beta$ values as described above, so that the
first gang has $\beta_{1}=2\times 10^{-5}$, whereas the second gang $2$ was
assigned a larger value of $\beta_{2}=3.5\times 10^{-5}$ and the third gang
was assigned a low value of $\beta_{3}=0.5\times 10^{-5}$. The results of the
simulations are presented in Figures 8, 9, 10, and 11, and also in Table 1.

Figure 8: Top row: Agent and graffiti densities’ temporal evolution for the
Timidity Model. Here $\beta_{1}=2\times 10^{-5},~{}\beta_{2}=3.5\times
10^{-5}$ and $\beta_{3}=0.5\times 10^{-5}$. We also have
$N_{1}=N_{2}=N_{3}=50,000$, with $\lambda=\gamma=0.5$, $\delta t=1$ and the
lattice size is $100\times 100$. It is clearly seen that the agents segregate
over time. Bottom row: The densities for gangs 1 (left), 2 (middle), and 3
(right) can be seen after $100,000$ time steps.

From Figure 8, which shows the temporal evolution of the agent and graffiti
densities, we can see that the system does segregate over time, however the
segregation differs from the original discrete model simulations in Section
3.2. We can see that the agents from the gang with the largest $\beta_{j}$
value, gang $2$, cluster tightly together into small, highly dense spots and
do not venture outside these spots. This is because they are the most strongly
avoidant of the other gangs’ graffiti, so they are the most timid. Most of the
agents from gang $1$, which has the next highest $\beta_{j}$ value, also
gather into fairly dense groups, motivated by avoiding the graffiti of Gangs
$2$ and $3$. However, because $\beta_{1}$ is less strong than $\beta_{2}$, a
smattering of gang $1$ agents can also be seen spreading roughly evenly over
the whole domain aside from the area occupied by gang $2$. The area occupied
by gang $2$ is avoided by all other gangs because of the high concentration of
graffiti laid down by the strongly localized agents. Gang $3$’s agents wander
more freely but still avoid the areas with denser graffiti, avoiding gang
$2$’s area more strongly than gang $1$’s area due to the higher concentration
of graffiti there. But gang $3$’s low $\beta_{3}$ allows them to spread over
much more of the territory, hence dominating more of the lattice than the
other two gangs.

Figure 11 shows cross-sectional slices of the lattice, in order to more
clearly show the agent and graffiti density for each gang. On the left, we see
the agent (top) and graffiti (bottom) densities for the Timidity Model. We can
again observe that the gang with the highest $\beta_{j}$ value, gang $2$, has
the smallest and densest territory, with a high density of graffiti and little
interference from the other gangs inside this territory. Gang $1$, with the
next-largest $\beta_{j}$ value, has a larger and less distinct territory, with
a medium graffiti density, while gang $3$, with the smallest $\beta_{j}$, is
dominating a very large but fairly mixed territory. We can see agents from all
gangs coexisting at different densities in the area dominated by gang $3$ due
to the lower graffiti concentration there.

We also use the same order parameter that we employed with the regular
discrete model to evaluate this variation on the model, and the results are
presented in the plot on the left in Figure 10. Based on our order parameter,
we find that the system does indeed show signs of segregation. We also note
that the order parameter does not scale to the value of one; this is because
the definition was based on all gangs having approximately the same area in
the final segregated state.

In Table 1, we consider three-gang simulations with six different sets of
parameters and tabulate how much of the territory at equilibrium is dominated
by each of the gangs. The $\beta_{j}$ values are listed in the third column
and we focus here on the percentage of the territory is listed in the fourth
column (the fifth column contains information on the percentage of territory
at equilibrium for the second variation of the model, discussed in the
subsequent subsection). We can see from the table that the percentage of
dominated territory has an inverse relationship with the value of $\beta_{j}$.

To better examine this relationship, in Figure 9, we plot the $\beta_{j}$
values against the percentage of territory dominated by the corresponding
gang. We can see that the territory percentage is roughly inversely
proportional to the $\beta_{j}$ value, meaning that, in the parameter regime
where territories form, one can expect this model to produce larger
territories for those gangs with smaller $\beta_{j}$. This is an important
feature of this variation at an ecological level.

Parameter set | Gang | Value of $\beta_{j}$ | % Territory, Model 1 | % Territory, Model 2
---|---|---|---|---
Set 1 | Gang $1$ | $\beta_{1}=0.000005$ |  $55.02$% |  $11.27$%
Gang $2$ | $\beta_{2}=0.00002$ |  $28.23$% |  $32.48$%
Gang $3$ | $\beta_{3}=0.000035$ |  $10.20$% |  $54.10$%
Set 2 | Gang $1$ | $\beta_{1}=0.000015$ |  $41.50$% |  $25.95$%
Gang $2$ | $\beta_{2}=0.00002$ |  $31.28$% |  $32.62$%
Gang $3$ | $\beta_{3}=0.000025$ |  $25.19$% |  $39.70$%
Set 3 | Gang $1$ | $\beta_{1}=0.00001$ |  $55.13$% |  $16.27$%
Gang $2$ | $\beta_{2}=0.00002$ |  $28.40$% |  $28.63$%
Gang $3$ | $\beta_{3}=0.00004$ |  $14.15$% |  $53.92$%
Set 4 | Gang $1$ | $\beta_{1}=0.000012$ |  $51.45$% |  $19.39$%
Gang $2$ | $\beta_{2}=0.000024$ |  $27.02$% |  $34.85$%
Gang $3$ | $\beta_{3}=0.000032$ |  $19.99$% |  $44.58$%
Set 5 | Gang $1$ | $\beta_{1}=0.000022$ |  $33.11$% |  $32.72$%
Gang $2$ | $\beta_{2}=0.000022$ |  $32.71$% |  $32.84$%
Gang $3$ | $\beta_{3}=0.000022$ |  $32.92$% |  $33.04$%
Set 6 | Gang $1$ | $\beta_{1}=0.000018$ |  $44.85$% |  $23.66$%
Gang $2$ | $\beta_{2}=0.000028$ |  $29.50$% |  $34.32$%
Gang $3$ | $\beta_{3}=0.000034$ |  $24.89$% |  $41.16$%
Table 1: Here, we see the results of both variations of the original model for
six different sets of $\beta_{j}$ in three-gang simulations. Here, Model 1
refers to the Timidity Model variation, while Model 2 refers to the Threat
Level Model variation. The $\beta_{j}$ values are listed, along with the
percentage of the lattice dominated by each gang at equilibrium. Note that the
percentages do not add to $100$% because in each simulation, a small
percentage of the lattice is not clearly dominated by any one of the gangs.

Figure 9: Here, we plot graphs of the beta values $\beta_{j}$ against the
percentage of the area dominated by gang $j$ for the Timidity Model (variation
1). We use six sets of parameters, enumerated in Table 1.

Figure 10: How changing the $\beta$ parameter for different gangs affects the
system. Here we have $N_{1}=N_{2}=N_{3}=50,000$, with $\lambda=\gamma=0.5$ and
the lattice size is $100\times 100$. We note that we cannot expect the order
parameter to approach $1$ even in the case of perfect segregation, since the
territories vary in size depending on the values of the $\beta_{j}$s. (Left):
Here, we plot the order parameter for the Timidity Model, using
$\beta_{1}=2\times 10^{-5},\beta_{2}=3.5\times 10^{-5}$ and
$\beta_{3}=0.5\times 10^{-5}$. (Right): Here, we plot the order parameter for
the Threat Level Model, where the $\beta_{i}$ parameters represent threat
levels; we use $\beta_{1}=2\times 10^{-5},\beta_{2}=3.5\times 10^{-5}$ and
$\beta_{3}=0.5\times 10^{-5}$. As in the original model, we see that our order
parameter behaves similarly, increasing quickly before leveling off in both
cases as the system segregates.

Figure 11: Cross-sectional slices of the agent and graffiti densities for
different $\beta$ extensions at the final time step for a segregated state.
Here we have $N_{1}=N_{2}=N=3=50,000$ with $\delta t=1$ and the lattice size
is $100\times 100$; in both simulations, $\beta_{1}=2\times
10^{-5},\beta_{2}=3.5\times 10^{-5}$ and $\beta_{3}=0.5\times 10^{-5}$.
(Left): Here, we consider the Timidity variation of the model. We observe that
the territories range from small and very dense, with little incursion from
the other gangs, to large and spread out, with other gang members encroaching
on the territory, as the gangs’ $\beta_{j}$ value varies from high to low.
(Right): Here, we consider the Threat Level variation of the model. We see
that the size of the territory here is correlated with the $\beta_{j}$ value
for the gang, and that all of the territories here seem well-defined, with
little of the territorial encroachment seen in the model pictured on the left.

###  6.2 Threat Level Model (Variation 2)

We now consider a different modification of movement dynamics (3). This model
is intended to apply in a situation where some gangs are more aggressive or
territorial than others. So instead of considering a $\beta$ value which is
the same for all gangs, we consider the case where the gangs have varying
threat levels. To this end, each gang $i$ has a corresponding threat level
encoded by parameter $\beta_{i}$. This means that gang $j$ will more strongly
avoid more threatening gangs, i.e. those gangs with relatively large $\beta$
values. Based on this, we must modify the opposition sum from equation (2), so
that it becomes

$\psi_{j}(x,y,t):=\sum_{\begin{subarray}{c}i=1\\\ i\neq
j\end{subarray}}^{K}\beta_{i}\xi_{i}(x,y,t).$ (42)

Note that the $\beta_{i}$ parameters can no longer pull out of the sum. The
new movement probability then becomes

$M_{j}(x_{1}\rightarrow x_{2},y_{1}\rightarrow
y_{2},t)=\frac{e^{-\psi_{j}(x_{2},y_{2},t)}}{\sum\limits_{(\tilde{x},\tilde{y})\sim(x_{1},y_{1})}e^{-\psi_{j}(\tilde{x},\tilde{y},t)}}.$
(43)

Here, every gang then avoids the graffiti of gang $i$ with rate $\beta_{i}$.
This model applies in the case where the gangs have differing threat levels,
so that some gangs are to be avoided more than others. For example, let us
suppose that gang $2$ has a large $\beta_{2}$ value, gang $3$ has a small
$\beta_{3}$ value, and gang $1$ has an intermediate $\beta_{1}$ value. As
$\beta_{2}$ is large, gang $2$’s territory will be strongly avoided by both
gangs $1$ and $3$. Furthermore, since gang $3$ has a small threat level
$\beta_{3}$, its graffiti will not be avoided as much by the other gangs and
it will need a higher graffiti density in order to claim territory for itself.

If we follow the same steps used to derive the continuum equations in Section
4, now substituting (3) with (43), it can easily be shown that the resulting
system of equations for $j=1,2,\dots,K$ are

$\begin{cases}\displaystyle\frac{\partial\xi_{j}}{\partial
t}(x,y,t)=\gamma\rho_{j}(x,y,t)-\lambda\xi_{j}(x,y,t)\\\
\displaystyle\frac{\partial\rho_{j}}{\partial
t}(x,y,t)=\frac{D}{4}\nabla\cdot\left[\nabla\rho_{j}(x,y,t)+2\left(\rho_{j}(x,y,t)\nabla\left(\sum_{\begin{subarray}{c}i=1\\\
i\neq
j\end{subarray}}^{K}\beta_{i}\xi_{i}(x,y,t)\right)\right)\right]\end{cases}$
(44)

with periodic boundary conditions. Note that the parameters $\beta_{i}$ now
cannot be pulled to the front of the second term of the second equation, and
instead must remain inside the sum.

To test these changes with our discrete model, we ran our simulations with
three gangs $1,2$ and $3$, where all gangs are assumed to have $50,000$
agents. We assume that the lattice size $L\times L$ is equal to $100\times
100$, and use $100,000$ time steps with each step size $\delta t=1$. We
assigned the first gang to have $\beta_{1}=2\times 10^{-5}$, the second gang
$2$ to have a larger value of $\beta_{2}=3.5\times 10^{-5}$, while the third
gang is assigned a low value of $\beta_{3}=0.5\times 10^{-5}$. The results of
these simulations are presented in Figures 10, 11, 12, 13, as well as Table 1.

Figure 12: Top row: Temporal evolution of the agent and graffiti densities for
the Threat Level Model. Here $\beta_{1}=2\times 10^{-5},\beta_{2}=3.5\times
10^{-5}$ and $\beta_{3}=0.5\times 10^{-5}$. We also have
$N_{1}=N_{2}=N_{3}=50,000$, with $\lambda=\gamma=0.5$, $\delta t=1$ and the
lattice size is $100\times 100$. It is clearly seen that the agents segregate
over time. Bottom row: The agent and graffiti densities for gangs 1 (left), 2
(middle), and 3 (right) can be seen after $100,000$ time steps.

From Figure 12, we can see that the system can segregate over time, in the
right parameter regime. This segregation, however, differs both from that of
the discrete model in Section 3.2 and from that of the previous subsection.
Here, we see that the gang with the largest $\beta$ value, whose territory
appears in blue in the top row of Figure 12, has the largest and least dense
territory. This is reasonable since the other gangs avoid the graffiti of gang
1 quite strongly; therefore, the gang does not need to put down as much
graffiti to maintain a territory. They can then spread over more space and
still maintain their territory. The gang with the smallest $\beta$ value, on
the other hand, whose color is green in the top row of Figure 12, clearly has
the smallest and most dense territory. This makes sense, since the other gangs
are not avoiding the territory of gang $3$ very strongly; gang $3$ then has to
put down a much higher density of graffiti to force the other gangs to avoid
it, and it can only do this by limiting its gang members to a smaller area.

We also tested segregation using the same order parameter to that we used
previously; the evolution of the order parameter for this model is presented
on the right in Figure 10. We note that, as in the previous variation of the
model, we can no longer expect the order parameter to tend to $1$ in the fully
segregated case. However, we do still see segregation over time.

Figure 11 shows cross-sectional slices of the lattice, to show the agent and
graffiti density for each gang. On the right, we see the agent (top) and
graffiti (bottom) densities. From this figure, we can see that the territories
formed in this variation are much more distinct than in the last variation;
there is very little overlap inside the territories. This is in contrast to
the first variation on the model. We can also observe that the $\beta_{j}$
value seems to be proportional to the territory size. Traveling outside an
agent’s own territory seemingly happens only along the boundaries of other
gangs’ territories.

In Table 1, as described in the previous subsection, we see the results of
this model run with three gangs. We ran the simulation with six different sets
of $\beta_{1},\beta_{2},$ and $\beta_{3}$ and, in the right-hand column of the
table, we see the percentage of the lattice occupied at steady-state by each
of the three gangs. We can see that in this variation of the model, in
contrast to the last variation, the size of the territory in each simulation
seems to be directly proportional to the values of $\beta_{i}$.

We further examine this result in Figure 13, where we plot the values of
$\beta_{i}$ for each simulation against the percentage of the lattice occupied
by each of the gangs. We see in this figure that the $\beta_{i}$ and the
percentage of occupied areas are indeed very nearly directly proportional. It
is an interesting open question why this is the case.

Figure 13: Here, we plot the $beta_{j}$ values against the percentage of the
area dominated by gang $j$ for the Threat Level variation of the model. We use
the six sets of parameters enumerated in Table 1.

##  7 Discussion

In this work, we have presented an extension of a previous agent-based system
that models gang territorial development “motivated” by graffiti tagging [25]
to now include a finite number $K$ of gangs as opposed to only two. In the
special case of three gangs, we have shown by using numerical simulation that
our model also undergoes a phase transition as we change the value of
different parameters. We formally derived the continuum limit for our model,
giving us a set of $2\times K$ convection-diffusion equations with cross-
diffusion. By using linear stability analysis on the continuum equations, we
showed that there is a bifurcation point in which the well-mixed state becomes
linearly unstable. Furthermore, we have numerically shown that the bifurcation
point matches the critical parameter found in the numerical simulations for
the case of $K=3$ for the discrete model. This generalization from two to $K$
gangs makes the model much more flexible. In the form presented in this paper,
the model can be applied to many coexisting gangs or many packs of animals,
and this is important in practice, since it can rarely be assured that there
are only two.

We have also presented two novel variations of the model, each of which
exhibits different segregation dynamics from the original model and from the
other variation. These variations allow for further flexibility. For the
Timidity model (variation 1), each gang is allowed a different value of the
$\beta$ parameter, allowing some more timid gangs (with large $\beta$) to be
more sensitive to the existence of graffiti and some (with small $\beta$) to
be less sensitive. Assuming the gangs have identical membership, this resulted
in the more timid gangs having smaller and more distinct territories, while
the less timid gangs had larger and less distinct territories where members of
other gangs were also occasionally present. For the Threat Level model
(variation 2), each gang $i$ has a threat level $\beta_{i}$ associated to
their graffiti, so that other gangs react more strongly to the graffiti of
gangs with a large $\beta_{i}$ and less strongly to those with a small
$\beta_{i}$. When gangs have identical membership, this variation results in
larger territories for gangs with higher threat level $\beta_{i}$ and smaller
territories for gangs with lower threat levels. In contrast to the Timidity
model, all of the territories are distinct, with very little overlap from
other gangs’ agents. These two variations could prove useful in ecological
applications where more is known about the traits of the groups.

The model is also intriguing from the perspective of pattern formation. The
segregation dynamics for the system with constant $\beta$ and the two
variations give three different dynamics for the territory formation. These
new models open the possibility of further studies, such as comparing pattern
formation with similarly segregating systems such as Cahn-Hilliard [32].
Additionally, this model exhibits a phase transition from non-segregating
populations to segregating populations as $\beta$ changes; it is highly likely
that a phase transition would also occur as $\lambda$ increases. An open
problem with significant ecological consequences would be to look for this
phase transition, since it would provide an indication that climate change, in
particular increased precipitation, could have an effect on the territorial
dynamics for animals such as wolves and coyotes.

The system of PDEs derived in this paper also are interesting in their own
right. The form is reminiscent of Patlak-Keller-Segel model [33, 34], with
chemo-repellent rather than chemo-attractant and no diffusion of the chemical.
The graffiti densities evolve in response only to the agent and graffiti
densities of the corresponding gang, while the agent densities evolve only in
response to the corresponding gang’s agent density and the graffiti densities
of all the other gangs. This leads to a system’s cross-diffusion form.
Originating in spatial ecology [35, 36, 37], cross-diffusion is widely
recognized as a mechanism for pattern formation [38]. Recent interest in
cross-diffusion has led to advances in analytical understanding of these
systems [39, 40, 41, 42]. Since this paper offers three variations on a novel
cross-diffusion system, new avenues are opened for further numerical and
analytical study to better understand the properties and behavior of these
systems, such as the analytical work done on the two-gang system [26].

##  8 Acknowledgements

The authors would like to thank Nancy Rodriguez, Havva Yoldas, and Nicola
Zamponi for helpful discussions of the original model upon which this paper is
based.

## References

  * [1] Ethan J Temeles.  The role of neighbours in territorial systems: when are they’dear enemies’?  Animal Behaviour, 47(2):339–350, 1994.
  * [2] Robert David Sack.  Human territoriality: its theory and history, volume 7.  CUP Archive, 1986.
  * [3] H Jochen Schenk, Ragan Morrison Callaway, and BE Mahall.  Spatial root segregation: are plants territorial?  Advances in ecological research, 28:145–180, 1999.
  * [4] FE May and JE Ash.  An assessment of the allelopathic potential of eucalyptus.  Australian journal of botany, 38(3):245–254, 1990.
  * [5] Paul R. Moorcroft, Mark A. Lewis, and Robert L. Crabtree.  Home range analysis using a mechanistic home range model.  Ecology, 80(50):1656–1665, 7 1999.
  * [6] Paul R. Moorcroft, Mark A. Lewis, and Robert L. Crabtree.  Mechanistic home range models capture spatial patterns and dynamics of coyote territories in Yellowstone.  Proceedings of The Royal Society B, 273:1651–1659, 2006.
  * [7] Roger P Peters and L David Mech.  Scent marking in wolves.  American Scientist, 63(6):628–637, 1975.
  * [8] M.A. Lewis, K.A.J. White, and J.D. Murray.  Analysis of a model for wolf territories.  Journal of Mathematical Biology, 35:749–774, 1997.
  * [9] Laura M. Smith, Andrea L. Bertozzi, P. Jeffrey Brantingham, George E. Tita, and Matthew Valasik.  Adaptation of an ecological territiorial model to street gang spatial patterns in Los Angeles.  Discrete amd Continuous Dynamical Systems, 32(9):3223–3244, 2012\.
  * [10] Rachel A. Hegemann, Laura M. Smith, Alethea B.T. Barbaro, Andrea L. Bertozzi, Shannon E. Reid, and George E. Tita.  Geographical influences of an emerging network of gang rivalries.  Physica A: Statistical Mechanics and its Applications, 390(21):3894–3914, 2011.
  * [11] Alethea B.T. Barbaro, Lincoln Chayes, and Maria R. D’Orsogna.  Territorial developments based on graffiti: A statistical mechanics approach.  Physica A, 392(1):252–270, 2013.
  * [12] Ernst Ising.  Beitrag zur theorie des ferromagnetismus.  Zeitschrift für Physik A Hadrons and Nuclei, 31(1):253–258, 1925.
  * [13] Yves van Gennip, Blake Hunter, Raymond Ahn, Peter Elliott, Kyle Luh, Megan Halvorson, Shannon Reid, Matthew Valasik, James Wo, George E. Tita, Andrea L. Bertozzi, and P. Jeffrey Brantingham.  Community detection using spectral clustering on sparse geosocial data.  SIAM Journal on Applied Mathematics, 73(1):67–83, 2013.
  * [14] Martin B. Short, Maria R. D’Orsogna, Virginia B. Pasour, George E. Tita, P.Jeffrey Brantingham, Andrea L. Bertozzi, and Lincoln B. Chayes.  A statistical model of criminal behavior.  Mathematical Models and Methods in Applied Sciences, 18(supp01):1249–1267, 2008.
  * [15] N Rodríguez.  On the global well-posedness theory for a class of pde models for criminal activity.  Physica D: Nonlinear Phenomena, 260:191–200, 2013.
  * [16] Nancy Rodriguez and Andrea Bertozzi.  Local existence and uniqueness of solutions to a pde model for criminal behavior.  Mathematical Models and Methods in Applied Sciences, 20(supp01):1425–1457, 2010.
  * [17] Henri Berestycki, Nancy Rodriguez, and Lenya Ryzhik.  Traveling wave solutions in a reaction-diffusion model for criminal activity.  Multiscale Modeling & Simulation, 11(4):1097–1126, 2013.
  * [18] Paul A. Jones, P. Jeffrey Brantingham, and Lincoln R. Chayes.  Statistical models of criminal behavior: The effects of law enforcement actions.  Mathematical Models and Methods in Applied Sciences, 20:1397–1423, 2010.
  * [19] Joseph R Zipkin, Martin B Short, and Andrea L Bertozzi.  Cops on the dots in a mathematical model of urban crime and police response.  Discrete Contin. Dyn. Syst. Ser. B, 19(5):1479–1506, 2014.
  * [20] Linfeng Mei and Juncheng Wei.  The existence and stability of spike solutions for a chemotax is system modeling crime pattern formation.  Mathematical Models and Methods in Applied Sciences, 30(9):1727–1764, 2020.
  * [21] Chuntian Wang, Yuan Zhang, Andrea L. Bertozzi, and Martin B. Short.  A stochastic-statistical residential burglary model with independent poisson clocks.  European Journal of Applied Mathematics, 32:35–38, 2021.
  * [22] H Berestyki and N Rodríguez.  Analysis of a heterogeneous model for riot dynamics: the effect of censorship of information.  European Journal of Applied Mathematics, 27(3):554, 2016.
  * [23] Nancy Rodríguez and Lenya Ryzhik.  Exploring the effects of social preference, economic disparity, and heterogeneous environments on segregation.  Communications in Mathematical Sciences, 14(2):363–387, 2016.
  * [24] Maria R D’Orsogna and Matjaž Perc.  Statistical physics of crime: A review.  Physics of life reviews, 12:1–21, 2015.
  * [25] Abdulaziz Alsenafi and Alethea B.T. Barbaro.  A convection–diffusion model for gang territoriality.  Physica A: Statistical Mechanics and its Applications, 510:765–786, 2018.
  * [26] Alethea BT Barbaro, Nancy Rodriguez, Havva Yoldaş, and Nicola Zamponi.  Analysis of a cross-diffusion model for rival gangs interaction in a city.  arXiv preprint arXiv:2009.04189, 2020.
  * [27] Rodney J Baxter.  Exactly solved models in statistical mechanics.  Courier Corporation, 2007.
  * [28] Yasmin Dolak and Christian Schmeiser.  Kinetic models for chemotaxis: Hydrodynamic limits and spatio-temporal mechanisms.  Journal of mathematical biology, 51(6):595–615, 2005.
  * [29] Martin B. Short, Andrea L. Bertozzi, and P. Jeffrey Brantingham.  Nonlinear patterns in urban crime: Hotspots, bifurcations, and suppression.  SIAM Journal on Applied Dynamical Systems, 9(2):462–483, 2010.
  * [30] Brian K. Briscoe, Mark A. Lewis, and Stephen E. Parrish.  Home range formation in wolves due to scent making.  Bulletin of Mathematical Biology, 64:261–284, 2002.
  * [31] K.A.J. White, M.A. Lewis, and J.D. Murray.  A model for wolf-pack territory formation and maintenance.  Journal of Theoretical Biology, 178(1):29–43, 1996.
  * [32] John W. Cahn and John E. Hilliard.  Free energy of a nonuniform system. I. Interfacial free energy.  The Journal of chemical physics, 28(2):258–267, 1958.
  * [33] Clifford S Patlak.  Random walk with persistence and external bias.  The bulletin of mathematical biophysics, 15(3):311–338, 1953.
  * [34] Evelyn F Keller and Lee A Segel.  Initiation of slime mold aggregation viewed as an instability.  Journal of theoretical biology, 26(3):399–415, 1970.
  * [35] M Morisita.  Population density and dispersal of a water strider. gerris lacustris: Observations and considerations on animal aggregations.  Contributions on Physiology and Ecology, Kyoto University, 65:1–149, 1950.
  * [36] Masaaki Morisita.  Habitat preference and evaluation of environment of an animal. experimental studies on the population density of an antlion, glenuroides japonicus m’l.[= correctly hagenomyia micans]. i.  Physiology and Ecology, 5(1):1–16, 1952.
  * [37] Morton E Gurtin and AC Pipkin.  A note on interacting populations that disperse to avoid crowding.  Quarterly of Applied Mathematics, 42(1):87–94, 1984.
  * [38] Vladimir K Vanag and Irving R Epstein.  Cross-diffusion and pattern formation in reaction–diffusion systems.  Physical Chemistry Chemical Physics, 11(6):897–912, 2009.
  * [39] Martin Burger, José A Carrillo, Jan-Frederik Pietschmann, and Markus Schmidtchen.  Segregation effects and gap formation in cross-diffusion models.  Interfaces and Free Boundaries, 22(2):175–203, 2020.
  * [40] Marco Di Francesco, Antonio Esposito, and Simone Fagioli.  Nonlinear degenerate cross-diffusion systems with nonlocal interaction.  Nonlinear Analysis, 169:94–117, 2018.
  * [41] José A Carrillo, Yanghong Huang, and Markus Schmidtchen.  Zoology of a nonlocal cross-diffusion model for two species.  SIAM Journal on Applied Mathematics, 78(2):1078–1104, 2018.
  * [42] Maria Bruna, Martin Burger, Helene Ranetbauer, and Marie-Therese Wolfram.  Cross-diffusion systems with excluded-volume effects and asymptotic gradient flow structures.  Journal of Nonlinear Science, 27(2):687–719, 2017.