A Multispecies Cross-Diffusion Model for Territorial Development

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.


Introduction
Many types of organisms are known to exhibit territorially. Examples include insects, fish, amphibians, reptiles, birds, and mammals [1], and, of course, hu-man 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: where ξ i is the graffiti density of gang i and ρ 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 β 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 β is made gang-dependent. Finally, in Section 7, we conclude with a discussion of the results and open problems.

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 × L, with periodic boundary conditions and area 1. We assume that we have K gangs, 1, 2, . . . 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 : 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 ρ j (x, y, t) = nj (x,y,t) l 2 , where l = 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 ξ j (x, y, t) = gj (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 (x,ỹ) ∼ (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: 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 ) ∈ S to one of the neighboring sites s 2 = (x 2 , y 2 ) ∈ S is defined to be again defined analogously to [25]. Here the parameter β 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 The expected density of gang j is therefore For the graffiti density update rules, each agent adds graffiti at its current site with probability γ. It is also assumed that the graffiti decays at every site with a rate of λ. Both the graffiti addition and decay are scaled by the time step δt. Therefore, the graffiti density at site (x, y) ∈ S at time t + δt is 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.

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 β, 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.

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: 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 } K j=1 form a partition for the lattice S for j = 1, 2, . . . , 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: 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: Thus, the expected agents density within S j is 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 = 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

An Order Parameter
To investigate the phase transition, we define the following order parameter: 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 = N 3 for j in {1, 2, 3}. The order parameter for this special case is 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 Simplifying the terms in the brackets yields that 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 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: Using equation (8), we substitute the expectation of each ρ i for each region, which gives us the following approximation: Further simplifying yields: 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 = L 2 3 , the previous equation can be further simplified to By further simplifying the equation, we easily see that Therefore, the order parameter for the perfectly segregated system is approximately equal to one.

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 × L is 100 × 100 with lattice spacing 2 l = 1, and we will use 100, 000 time steps with each step size δt = 1.

Well-Mixed State
We start our simulations with β = 5 × 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 λ = γ = 0.5, β = 5×10 −6 , δt = 1 and the lattice size is 100×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 β = 5 × 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 β value. This is due to the way the agents are allowed to move in equation (3), where a very small β values give the agent a probability of nearly 0.25 to move to one of the four neighboring sites.

Segregated State
The value of β is now increased so that it is equal to 3 × 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. 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 β value is large enough that the agents are  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 λ = γ = 0.5, β = 3 × 10 −5 , δt = 1 and the lattice size is 100 × 100. We see that the agents segregate into distinct territories, coarsening over time.
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.

Effects of β
In Sections 3.1 and 3.2, we saw that changing the value of the parameter β 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 parameter E ≈ 0 for a well mixed state and E ≈ 1 for a fully segregated state. For our simulations, we graphed the order parameter over the course of the simulation for different values of β, visualizing the output in Figure 3. (Left) It is seen that for a small β value the system remains well-mixed and the order parameter is almost zero over all time steps. For larger β 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 β values. We clearly see that as the β value increases there is a critical β 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 β 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 β = 0 and β = 0.000005, the order parameter remained approximately zero throughout all time steps. This is expected as the system remains wellmixed for these relatively small β values. However, we see that once we increase the values of β, then the order parameter starts to increase. For instance, if β = 0.000025 or β = 0.000030, then the order parameter increases fairly quickly in the first 10, 000 time steps before leveling off to just under the fullysegregated value of 1 for the remaining time steps. This shows us that for these relatively large β values, the system segregates fairly quickly and remains segregated throughout the simulation. Finally, we also see that if we choose β = 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 β in which the system undergoes a phase transition. We define the critical β to be the value where the order parameter is equal to 0.01, and denote it by β * . To find the value of β * , we take the final value of the order parameter and plot it against different β 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 β * ∈ (0.000005, 0.000006).

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 δ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 δt from 1 to 0.1; we then plot the final order parameter value for different β 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 β at which the phase transition occurs is about 1.8×10 −5 . However, in the middle plot, the mass is increased to 150, 000 and this time the phase transition occurs around 0.9 × 10 −5 . Thus, we notice that as the systems' mass increases, the resulting phase transition happens at a smaller β 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 β value should be sufficient for the agents to react to the graffiti field. Here we have N 1 = N 2 = N 3 = 50, 000, with λ = γ = 0.5, δt = 1 and the lattice size is 100 × 100. In all three plots, it is seen that for small β values, the system remains well-mixed and the resulting order parameter values is approximately zero over time. However, for larger β 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 β values. We clearly see that as the β value increases there is a critical β value in which a phase transition occurs. The three plots show the effects of changing the mass and the ratio γ/λ.
We also investigated how the ratio γ λ 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 λ. 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 γ λ ratio, a higher β value is needed for segregation. This is evident in the middle and right plots in Figure 4. There, we clearly see that when the γ λ = 1, the critical β is around 0.9 × 10 −5 , whereas when the γ λ is decreased to 0.5, the critical β is around 1.8 × 10 −5 . Physically, this is due to the fact that less graffiti on a site means that a larger β 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.

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 with periodic boundary conditions, where j ∈ {1, 2, . . . , 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.

Continuum Graffiti Density
We start by formally deriving the continuum equations for graffiti densities, recalling that for j ∈ {1, 2, . . . , K}, the discrete model (6): Rearranging the equation and dividing by δt gives us: This is now in the form of a difference equation. Assuming sufficient smoothness of the agent and graffiti densities ρ j and ξ j , we take δt → 0. This gives us the final form of the graffiti continuum equation for gang j:

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: Employing this notation, the first quantity we define is 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 ∇T j and ∆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 (x,ỹ, t) in the derivation of these quantities. We start by simplifying T j using Taylor series approximations. Recall the Taylor expansion We apply Taylor this expansion to T j with x = 4 and h = l 2 β 2 (∇ψ j ) 2 − β∆ψ j : Note that here, we are depending on the smoothness of ψ j . Then, by taking the gradient of (17), we find that We also can find ∆T j : 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 where (x,ỹ) 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 (x,ỹ) moves to site (x, y): where (x,ỹ) are the four neighbors of site (x,ỹ).
To remove the presence of the neighbors' neighbors (x,ỹ) from the denominator, we apply the discrete Laplacian to find thats Noting that Combining equations (21) and (22) gives us Substituting it back into equation (20), and replacing the denominator gives us The term inside the large brackets in the equation above takes the form of (16) where (x, y, t) is replaced with (x,ỹ, t), yielding the following approximation:

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 (23) to simplify the discrete model. Starting from the discrete model, we recall equation (5): Rearranging the equation and dividing both sides by δt gives us By equation (23), and noting that each agent has to move to one of the neighboring sites, Using the discrete Laplacian technique, we can approximate the contribution of the neighboring sites, giving on the right-hand side The notation (x, y, t) is again dropped as there are no longer any neighbors (x,ỹ, t) remaining in this derivation. Hence, we write, From definition (16), T j (x, y, t) is substituted back into the first term of (25): Simplifying the expression yields Using a Taylor series expansion on the first term within the brackets yields, Expression (26) thus becomes, Simplifying the expression yields However, we can further simplify this by noting that From (17) through (19), we have Therefore,

Substituting (28) back into (27) gives us
Assuming that the agent density ρ j is sufficiently smooth and the following limits gives us the final form for the continuum equations for the density of gang j agents: Finally, from (15) and (30), and using equation (2) to express everything in terms of agent density and graffiti density, the limiting convection-diffusion system for our model is for j = 1, 2, . . . K with periodic boundary conditions.

Steady-State Solutions
Considering steady-state solutions for the graffiti density, we find from the evolution equations for for the graffiti density that We now focus our attention on the steady-state solutions for the agent density of gang j: considering solutions of the form Using the steady-state graffiti density derived in equation (32), we find that Any form of ρ j (x, y, t) satisfying the above equation is a steady-state solution of our system. For simplicity, we suppose that ρ j is a constant for all j. In this case, the steady-state solution of our problem takes the form: for j = 1, 2, . . . , K and c j is a positive constant.
To test whether these steady-state solutions of the continuum system (34) agree with the discrete model, we first start our simulations in a steady-state solution from (34). 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.

Agent
Here we have N 1 = N 2 = N 3 = 50, 625, with λ = γ = 0.5, β = 3 × 10 −4 and the lattice size is 75 × 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 (34) we see that the graffiti density is ξ = γ λ ρ, which implies that the steady-state solution of the graffiti and agent densities at a site are equal up to scaled amount of γ λ . 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. From looking at the cross-sectional slices we clearly that ξ j ≈ γ λ ρ j , j = 1, 2, 3 and this agrees with our steady-state solution which was ξ j = γ λ ρ j .

Linear Stability Analysis
To have a better understanding of our system we linearize our model by considering a perturbation of the equilibrium solution (34) to the well-mixed state. We assume that our perturbations are of the form ǫ = δe αt e ikx , with δ << 1, and in that case our solution will take the following form: 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, α 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].
Next, we substitute (35) into the evolution equation Since we are working in one dimension and our equilibrium solution is a constant in both space and time, We can safely neglect the term O(δ ρj K l=1 l =j δ ξ l ), since |δ ρj |, |δ ξ l | << 1; therefore, Next, we write the equations from (37) and (38) in a systems form: For simplicity we consider the case where K = 3, with all gangs having the same β parameter and write the system in matrix vector format, giving us:

This gives us
which reduces to an eigenvalue problem for matrix F . For the problem to have a non trivial solution (i.e. δ = 0), the determinant of (F − αI 4 ) must be zero. Therefore, giving us the following characteristic polynomial Solving the characteristic polynomial gives six eigenvalues, which we solve numerically using Mathematica and plot the results for different values of β in Figure 7.
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 β 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 β values. However, when we increase the value of β 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 β from the discrete model matches those of the linearized system of partial differential equations. 6 Variations of the model: varying β 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 ) ∈ S to one of the neighboring sites s 2 = (x 2 , y 2 ) ∈ S to be from equation (2). The parameter β then controls how strongly each gang reacts to the graffiti of the other gangs. However, it is reasonable to consider that this parameter β 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.

Timidity Model (Variation 1)
In the first modification of the model, instead of having identical β values for all gangs, we change it so that gang j has a distinct corresponding β value, denoted by β j . This β j determines how much attention gang j places on the graffiti of the other gangs. In essence, this β j encodes the timidity of gang j, with higher β j corresponding to higher timidity, causing gang j to more strongly avoid other gangs' graffiti. Hence, the modified definition for movement becomes, In this variation of the model, gang j avoids all other gangs' graffiti with rate β 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 β 2 value, gang 3 has a relatively small β 3 value, and gang 1 has an intermediate β 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 β 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 (39), it can be easily shown that the resulting system of equations for j = 1, 2, . . . , K are with periodic boundary conditions. We can see that the β 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 β 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 × L is equal to 100 × 100, and will use 100, 000 time steps with each step size δt = 1. We assigned β values as described above, so that the first gang has β 1 = 2 × 10 −5 , whereas the second gang 2 was assigned a larger value of β 2 = 3.5 × 10 −5 and the third gang was assigned a low value of β 3 = 0.5 × 10 −5 . The results of the simulations are presented in Figures 8, 9, 10, and 11, and also in Table 1.
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 also have N 1 = N 2 = N 3 = 50, 000, with λ = γ = 0.5, δt = 1 and the lattice size is 100 × 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.
We can see that the agents from the gang with the largest β 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 β j value, also gather into fairly dense groups, motivated by avoiding the graffiti of Gangs 2 and 3. However, because β 1 is less strong than β 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 β 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 β 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 β j value, has a larger and less distinct territory, with a medium graffiti density, while gang 3, with the smallest β 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 β 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 β j .
To better examine this relationship, in Figure 9, we plot the β 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 β 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 β j . This is an important feature of this variation at an ecological level.

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 β 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 β i . This means that gang j will more strongly avoid more threatening gangs, i.e. those gangs with relatively large β values. Based on this, we must modify the opposition sum from equation (2), so that it becomes Note that the β i parameters can no longer pull out of the sum. The new movement probability then becomes Here, every gang then avoids the graffiti of gang i with rate β 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  Percentage of Area Dominated

Beta vs. Area
Parameter Set 6 Figure 9: Here, we plot graphs of the beta values β 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 β parameter for different gangs affects the system. Here we have N 1 = N 2 = N 3 = 50, 000, with λ = γ = 0.5 and the lattice size is 100 × 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 β j s. (Left): Here, we plot the order parameter for the Timidity Model, using β 1 = 2×10 −5 , β 2 = 3.5×10 −5 and β 3 = 0.5×10 −5 . (Right): Here, we plot the order parameter for the Threat Level Model, where the β i parameters represent threat levels; we use β 1 = 2 × 10 −5 , β 2 = 3.5 × 10 −5 and β 3 = 0.5 × 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. gang 2 has a large β 2 value, gang 3 has a small β 3 value, and gang 1 has an intermediate β 1 value. As β 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 β 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 (42), it can easily be shown that the resulting system of equations for j = 1, 2, . . . , K are (43) with periodic boundary conditions. Note that the parameters β 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×L is equal to 100×100, and use 100, 000 time steps with each step size δt = 1. We assigned the first gang to have β 1 = 2 × 10 −5 , 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' β 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 β 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.
the second gang 2 to have a larger value of β 2 = 3.5 × 10 −5 , while the third gang is assigned a low value of β 3 = 0.5 × 10 −5 . The results of these simulations are presented in Figures 10, 11, 12, 13, as well as Table 1.
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 β 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 We also have N 1 = N 2 = N 3 = 50, 000, with λ = γ = 0.5, δt = 1 and the lattice size is 100 × 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.
territory. The gang with the smallest β 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 β 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 β 1 , β 2 , and β 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

Beta vs. Area
Parameter Set 6 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.
the territory in each simulation seems to be directly proportional to the values of β i . We further examine this result in Figure 13, where we plot the values of β i for each simulation against the percentage of the lattice occupied by each of the gangs. We see in this figure that the β i and the percentage of occupied areas are indeed very nearly directly proportional. It is an interesting open question why this is the case.

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 × 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 β parameter, allowing some more timid gangs (with large β) to be more sensitive to the existence of graffiti and some (with small β) 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 β i associated to their graffiti, so that other gangs react more strongly to the graffiti of gangs with a large β i and less strongly to those with a small β i . When gangs have identical membership, this variation results in larger territories for gangs with higher threat level β 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 β 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 β changes; it is highly likely that a phase transition would also occur as λ 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 chemorepellent 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].