Swapping in lattice-based cell migration models

Cell migration is frequently modelled using on-lattice agent-based models (ABMs) that employ the excluded volume interaction. However, cells are also capable of exhibiting more complex cell-cell interactions, such as adhesion, repulsion, pulling, pushing and swapping. Although the first four of these have already been incorporated into mathematical models for cell migration, swapping has not been well studied in this context. In this paper, we develop an ABM for cell movement in which an active agent can 'swap' its position with another agent in its neighbourhood with a given swapping probability. We consider a two-species system for which we derive the corresponding macroscopic model and compare it with the average behaviour of the ABM. We see good agreement between the ABM and the macroscopic density. We also analyse the movement of agents at an individual level in the single-species as well as two-species scenarios to quantify the effects of swapping on an agent's motility.


I. INTRODUCTION
Cell migration is an essential biological process required for the correct development of tissues and organs during embryonic development and their proper maintenance, through wound healing and tissue homeostasis, throughout life [1][2][3][4][5]. During embryonic development, neural crest cells delaminate from the dorsal most aspect of the neural tube and migrate to colonise their target tissues including the gut in the case of enteric ganglia precursors and the skin in the case of melanoblasts the precursors of melanocytes [6].
Diseases of the neural crest are known as neurocristopathies for example failure of the enteric ganglia precursors to colonise the developing gut results in Hirschsprung's disease [7] while failure of melanoblasts to colonise the developing epidermis results in piebaldism [8]. Therefore an in-depth understanding of cell migration is important for identifying the causes of neurocristopathies [9][10][11][12] as well as developing new therapeutic targets to prevent metastasis in cancers [9,13].
Traditionally, many biological problems have been modelled using deterministic methods. However, in cell migration, randomness can play a salient role in determining a cell's trajectory and fate and hence deterministic theory may not be appropriate. Extensive research has gone into modelling cell movement as a stochastic process [14][15][16][17][18].
In one widely used approach, cells are modelled as agents whose positions evolve probabilistically in space and time according to a predefined set of rules. These models are commonly known as agent-based models (ABMs) or individual-based models (IBMs).
The agent-based modelling paradigm can be sub-divided into off-lattice and on-lattice models, both of which have wide applicability to different problems within mathematical biology. Gavagnin and Yates [19] recently reviewed the most commonly used ABMs for cell movement.
In this paper, we only concern ourselves with on-lattice models of cell movement. In a lattice-based approach, the domain is divided into a series of compartments in which the cells reside. Cells take up space, preventing other cells from occupying the same space at the same time. For biological plausibility, it is often desirable that mathematical models of cell migration account for the single occupancy of sites. This realism is incorporated in an ABM via the volume-exclusion principle, which states that a cell attempting to move into a neighbouring site successfully moves only if the neighbouring site is not already occupied at the time of moving. Models with volume exclusion at their core have been used to describe the collective Figure 1. Two Fucci2a labelled NIH 3T3 mouse embryonic fibroblasts swapping places with each other in culture [26]. The colour of the cells represents their cell cycle stage (Red = G1, Green = S/G2/M) in this case making it easy to observe the swap. In (a), we show the initial placement of cells: cell 1 (shown in red) is on the bottom-left of cell 2 (shown in green). In (b) a swap starts to take place and in (c) the swap is complete and now cell 1 is on the top-right of cell 2.
migration of cells for a wide range of biological applications. Mort et al. [8] used an on-lattice ABM with the exclusion principle to model the invasion of the developing epidermis by melanoblasts (the embryonic precursors of melanocytes) and investigate the basis of piebaldism in mice. Lattice-based exclusion models have also been applied to wound healing [20], migration of breast cancer cells [21], developmental processes on growing domains [22], cells' responses to chemotaxis [23] and cells exhibiting pushing [24] and pulling [25] interactions in densely crowded environments.
Biologically, although cells are excluded from the space occupied by other cells, their movement is not completely inhibited by them as typically assumed in volume-exclusion models. For example melanoblasts are able to move freely between keratinocytes in the developing epidermis [8]. Experimental data suggests that cells are often able to move past each other (passing laterally, above or below) exchanging places with one another.
In Figure 1 we show experimental images of two Fucci2a labelled NIH-3T3 fibroblasts exhibiting the swapping behaviour. Swapping has also been observed in blood cells such as leukocytes, erythrocytes and thrombocytes [27] and in pattern formation for maintaining sharp boundaries between different groups of cells as part of a cell sorting mechanism [28]. These examples highlight the importance of incorporating swapping into models of cell migration. To the best of our knowledge, this has not yet been explored thoroughly from a mathematical perspective.
In this paper, we develop a mathematical model to describe and analyse cell-cell swapping in two species setting. By modifying the movement rules of the traditional volumeexclusion process, we show that swapping between agents has an effect on the migration of agents at different spatial resolutions. To investigate how swapping manifests itself in the corresponding population-level model (PLM) we derive a set of partial differential equations (PDEs) describing the macroscopic dynamics of the agents. We compare numerical solutions of the PDEs with the averaged results from the ABM and comment on the agreement or discrepancy between them. We also analyse the movement of agents at an individual level and derive expressions for the individual-level diffusion coefficient.
The remainder of this paper is structured as follows. In Section II A, we develop a model that allows for swapping to take place between pairs of neighbouring agents. In Section II B, we derive the macroscopic PDEs and compare the average behaviour of the ABM to that of the PDEs. In Section III, we analyse the movement of agents at an individual level and derive a relationship between the individual-level diffusion coefficient, swapping probability and background domain density. In Section IV, we give examples to illustrate the applications of swapping and finally in Section V, we conclude the paper with a summary and discussion.

II. CELL MIGRATION MODEL WITH SWAPPING
We begin by developing an ABM for cell movement with swapping in Section II A and we use this to investigate the effect of swapping on the mobility of agents. In particular, we look at the effect of swapping on mixing of agents in a two-species system. We derive the population-level model and compare this with the average behaviour of the ABM in Section II B. We also analyse the individual-level behaviour in both single and multispecies scenarios in Section III.

A. On-lattice agent-based model
We model cell migration on a two-dimensional lattice. We discretise the domain into compartments (also known as 'sites') such that there are L x compartments in the horizontal direction and L y compartments in the vertical direction. We assume that the compartments are square with side length ∆. Supposing that each compartment can contain no more than one agent, ∆ can be considered a rough proxy for a cell's diameter.
A site (i, j) for i = 1, ..., L x and j = 1, ..., L y can be either occupied by a type-A or a type-B agent or unoccupied. Occupancy of a site (i, j) for a type-A (or type-B) agent is defined as a binary indicator, taking a value of 1 if there is a type-A (or type-B) agent (c) Figure 2. A schematic illustrating the swapping mechanism. Red sites are occupied with agents and black sites are unoccupied. The initial configuration of the lattice before is shown in (a). The agent chosen to move is at site (i, j) and labelled 1. The target site is at position (i + 1, j) and the agent occupying the target site is labelled 2. Agent 1 attempts to move into the target site in (b). The final configuration once the swapping move is complete is shown in (c).
at the site (i, j) or 0 if the site is empty.
We initialise the lattice with species A and species B agents at densities c A and c B , respectively, such that c A + c B = c and 0 ⩽ c ⩽ 1 where c is the overall domain density.
We let the positions of the agents evolve in continuous time according to the Gillespie algorithm [29]. Let r A be the rate of movement of a type-A agent and let r B be the equivalent for a type-B agent. The rates of movement are defined such that r A δt (and equivalently r B δt) is the probability that a type-A (or type-B) agent attempts to move during a finite time interval of duration δt. The agent attempts to move into one of the four sites in its Von Neumann neighbourhood with equal probability. If the chosen neighbouring site is empty, the focal agent successfully moves and its position is updated.
This blocking of the move characterises volume exclusion in our model.
Swapping works by modifying the rules of the exclusion process by allowing an exchange in the positions of a pair of neighbouring agents if the target site is already occupied. We now introduce the swapping parameter ρ to denote the probability of a successful swap between a pair of neighbouring agents conditional on one of the agents attempting to move into the other's position. If ρ = 0 then there are no swaps and we arrive back at the original exclusion process. If ρ > 0 then we can have different levels of swapping based on the value of ρ. For example, for ρ = 1 each scenario in which a move is attempted into an occupied target site will be a successful swap and for ρ = 0.5 half of the attempted moves into occupied sites will be successful.
To implement swapping, we sample a random number u from the uniform distribution over the unit interval (0, 1). If u < ρ, the agent at the site (i, j) swaps with the agent at the target site (Figure 2(a)-(c)), otherwise the move is aborted. Figure 2 shows a successful swap between two agents labelled '1' and '2'.
In this article we present results for the two-species model only unless stated otherwise.
Results for the single-species model can be found in Appendix A.
In Figure 3 we present snapshots of lattice occupancy of a 200 by 20 grid occupied with type-A (red) and type-B (green) agents. The two species move with equal rates, We see that a non-zero swapping probability results in faster dispersion of the agents compared to the ρ = 0 case. We also note that increasing the swapping probability from ρ = 0.5 the ρ = 1 results in faster dispersion of the agents, as expected.
In the next section, we derive the macroscopic PDEs describing the evolution of the mean lattice occupancy. By analysing the PDEs we generate further insight into the behaviours observed in Figure 3.

B. Continuum model
Let A k ij (t) be the occupancy of site (i, j) at time t on the kth repeat of the ABM such that A k ij (t) = 1 if the site is occupied by a type-A agent and 0 otherwise. Let B k ij (t) be the same for a type-B agent. Then the average density of type-A and type-B agents after K repeats is given by, In what follows we will typically drop the notation for time dependence of our species densities, i.e. ⟨A ij (t)⟩ = ⟨A ij ⟩ and ⟨B ij (t)⟩ = ⟨B ij ⟩, for conciseness. By considering all the possible ways in which the site (i, j) can gain or lose occupancy of either type-A or type-B agents during the time step δt, we can write down the corresponding occupancy master equations at time t + δt: For illustration, we describe the terms in Equation (2). The terms in Equation (3) carry similar interpretation. A site (i, j) can gain occupancy of type A in one of the following three ways: 1. The site (i, j) is unoccupied and a type-A agent moves in from a neighbouring site (line 1 in Equation (2)).
2. The site (i, j) is occupied by a type-B agent, which initiates and completes a swap with a type-A agent at a neighbouring site (line 3 in Equation (2)).
3. The site (i, j) is occupied by a type-B agent and a type-A agent at a neighbouring site initiates a swap to exchange positions with the type-B agent in the site (i, j) (also line 3 in Equation (2)).
In all three cases, a type-A agent moves into the site (i, j). Similarly, there are three ways for a type-A agent to move out of the site (i, j) leading to a loss in the corresponding occupancy: 1. The site is occupied by a type-A agent which jumps out to an unoccupied neighbouring site, leaving the site (i, j) empty (line 2 in Equation (2)).
2. The site (i, j) is occupied by a type-A agent which initiates a swap with a type-B agent in its neighbourhood.
3. The site (i, j) is occupied by a type-A agent and a type-B agent in the neighbouring site initiates a swap to exchange positions with the agent in site (i, j) (line 4 in Equation (2)).
In all three cases, a type-A agent moves out of the site (i, j).
To obtain the continuum model, we Taylor expand the appropriate terms on the RHS of Equations (2) and (3) around the site (i, j) keeping terms of up to second order. By letting ∆ → 0 and δt → 0 such that ∆ 2 /δt is held constant, we arrive at the coupled PDEs, where, Here, are the macroscopic diffusion coefficients corresponding to species A and B, respectively.
Setting ρ = 0 in equations (4) and (5) leads to, which are the macroscopic equations for the two-species volume-exclusion process [31].
In Figure 4, we compare the column-averaged density of the ABM given by, to the numerical solution of the one-dimensional analogue of Equations (4) and (5) with reflective boundary conditions by averaging the PDEs over the y direction [32]. In the ρ = 0 case ( Figure 4, first column), the PDE solutions and the ABM do not agree well as evidenced by the disparity between the two profiles. This discrepancy can also be seen in Simpson et al. [31] where the authors devise an ABM for multi-species exclusion processes (ρ = 0 case in our model) and compare their ABM with the corresponding continuum model. We remark that one reason for the discrepancy is that in crowded environments where the movement of agents is frequently inhibited by other agents, lattice occupancies cannot be considered independent of each other and spatial correlation are not dissipated efficiently [33][34][35]. Independence of lattice sites is a key assumption that is typically made when deriving the continuum models such as the one above [19,21,24,25,30,31,36].
For non-zero swapping probabilities, the agreement is significantly improved between the deterministic and stochastic model ( Figure 4, second and third columns). Swapping helps to break down the spatial correlations improving the agreement between the PLM and the ABM. We also see that in the zero swapping case, crowding of the green agents behind the red agents leads to profiles for which the maximum density at non-zero time (shown t = 100 and t = 1000 in Figure 4(a)) is higher than the initial maximum density (t = 0).
The reason for this is that for a multispecies migration process with cross-diffusion terms there is no maximum principle for the individual species [31,37,38]. We note that the enhanced diffusion which swapping engenders eliminates this effect. However, this does not necessarily mean that a maximum principle now holds for the systems under consideration. Investigating this further is beyond the scope of this article.

III. INDIVIDUAL-LEVEL ANALYSIS
In this section, we analyse the movement of agents at the individual level in the singlespecies and two-species case to assess how swapping affects the movement of agents.

A. Single-species individual-level analysis
The single-species swapping discrete and continuum models are given in Appendix A.
Here, we present the individual-level analysis for the single-species model. Our aim is to quantify the movement of individually tagged agents by their individual-level diffusion coefficient. For the analysis that we present next, we neglect any long range temporal correlations in the agents' movement. We first derive the individual-level time-uncorrelated diffusion coefficient analytically and then we compare it with its ABM approximation.
Let P ij (t) = P ij denote the probability that a focal agent is at position (i, j) at time t.
Defining δt as an infinitesimally small change in time, we can write down the probability of the agent being at position (i, j) at time t + δt as, A lattice site is occupied with probability c and empty with probability 1 − c. The first four terms in Equation (10)  with probability (1 − rδt) and the probability that the agent attempts to swap with an agent at an occupied neighbouring site but fails to complete the swap with probability rc(1 − ρ)δt and lastly the probability that a neighbouring agent successfully swaps with the agent at site (i, j) with probability rcρδt.
By rearranging, dividing both sides of Equation (10) by δt and taking the limit as δt → 0 leads to the system of ODEs given by, which describe the time evolution of the probability of finding an agent at position (i, j) at time t. From Equation (11) it can be shown that the expected net displacement of an agent is zero. This gives us no information about the statistical fluctuations in the movement of an agent and therefore we use the variance to quantify the net displacement [30,39]. The equations describing the time-evolution of the variances ⟨i(t) 2 ⟩ and ⟨j(t) 2 ⟩ are given by, and, Under the initial condition that at time t = 0, ⟨i 2 ⟩ = ⟨j 2 ⟩ = 0, Equations (12) and (13) solve to give, The time-uncorrelated individual-level diffusion coefficient can be retrieved as, As briefly discussed above, the master equation neglects temporal correlations in an agent's position. Consequently, the diffusion coefficient derived from the master equation will necessarily be inaccurate and fail to represent the true dynamics of the agents in the system. As a result, we refer to the expression given in Equation (15)  We analyse individual movement of agents using the sum of squared displacement (SSD) [30,36], and, By only taking the difference between successive positions, the SSD neglects temporal correlations in an agent's position which is consistent with our master equation in Equation (10). For each value of c and ρ we average the SSD over an ensemble of tracks and fit a linear model of the formŜ x t = a x t andŜ y t = a y t in each orthogonal direction x and y, respectively. The ABM approximation of the time-uncorrelated individual-level diffusion coefficient can be extracted from the gradient of these linear equations, where d = 2 is the dimension dimension of the lattice. We put a 'hat' symbol over D ⋆ to differentiate it from the exact time-uncorrelated diffusion coefficient D ⋆ in Equation (15).
In Figure 5, is as if the focal agent were on an unoccupied domain irrespective of density. For c = 1 and ρ = 1 we note that D ⋆ = 0.5 (i.e. twice as large compared to an agent moving on a domain with zero background density) as every attempted move by the focal agent is executed successfully and the focal agent is also moved equally often by neighbouring agents swapping into its position.

B. Two-species individual-level analysis
In this section, we perform the individual-level analysis for a two-species system as set out for the single-species case in Section III A.
Let P A ij (t) = P A ij be the probability that a focal agent of type-A occupies the position (i, j) and let P B ij (t) = P B ij be the equivalent for a type-B focal agent. Recalling δt as a small change in time we can write down the master equations for species A and B at time A full derivation of Equations (19) and (20) can be found in Appendix B accompanying this article.
The individual-level time-uncorrelated diffusion coefficients D ⋆ A and D ⋆ B for the species A and B, respectively, are given by, In order to investigate the effect of swapping on a two-species system and to verify our theoretical results, we simulate the two-species model, track a set of tagged agents and analyse their movement using the SSD as done in Section III A.
In Figure 6  the movement of the agents in the multi-species setting compared to the pure volumeexcluded scenario (ρ = 0). Furthermore, since species B has a higher movement rate than species A, species B diffuses faster than species A, apart from the c = 1 and ρ = 0 trivial case in which D ⋆ A and D ⋆ B are both 0 since the agents have nowhere to go.
In Figure 7 we present a similar comparison as in Figure 6 but this time the two species are not present in equal proportions. In this particular case, we consider a scenario in which type-A agents make up 25% of the total population and type-B agents 75%. Again, we see good agreement between the ABM and theoretical values and that swapping speeds up the movement of agents.
Comparing Figure 6 to

IV. ILLUSTRATIVE EXAMPLES
In this section, we show examples of the situations in which swapping has important applications. In Section IV A, we build the swapping mechanism into a cell migration model with proliferation and in Section IV B, we show how the swapping mechanism in conjunction with cell-cell adhesion can facilitate spontaneous pattern formation in densely crowded environments.

A. Swapping model with cell proliferation
We look at the role of swapping in cell migration with proliferation. For this example, we concern ourselves with the two-species cell migration model. The movement kinetics of the agents are the same as the swapping model described in Section II A but in addi-tion to migrating, agents can attempt to proliferate, placing a daughter at a randomly chosen neighbouring site if the site is empty, otherwise the division event is aborted. The proliferation rates per unit time for the two species A and B are denoted by r A p and r B p , respectively.
We initialise the domain with L x = 100 sites in the horizontal direction and L y = 20 sites in the vertical direction. We fill all the sites in the range 41 ⩽ x ⩽ 60 with agents of type A and all the remaining sites with type-B agents at a density of 0.5. The movement rates of agents are set to r A = r B = 1 and the proliferation rates as r A p = 0.01 and r B p = 0, i.e. only the agents of type A divide and the number of type-B agents are held constant. We let the system evolve according to the specified ABM.
We see, at the same time points, that cells are unsurprisingly more well-mixed in the case of non-zero swapping (second and third columns in Figure 8) compared to the zeroswapping situation (first column in Figure 8). We also see faster colonisation of the domain overall in the non-zero swapping cases than without swapping. This is because swapping allows the proliferating red agents to disperse more quickly into less dense regions, which in turn increases the probability of a successful division events for these agents. Without swapping, it takes longer for proliferative red agents to find the space to proliferate into. This trend of decreasing colonisation time with increasing swapping probability is reinforced in Figure 9 where we see that the time to reach the domain's carrying capacity is a decreasing function of the swapping probability, ρ.
The mean-field PDEs describing the approximate population-level dynamics of the agent are given by, where D 1 , D 2 , D 3 and D 4 are as defined previously. Notice that proliferation of type-A agents gives rise to an additive source term in Equation (23). The derivation of this source is standard and can be found in Plank and Simpson [40], Simpson et al. [41], for example.  In Figure 10 we compare the average column density of the ABM with the numerical solution of the one-dimensional analogue of the mean-field PDEs obtained by averaging over the y direction with r A p = 0.01. We chose this value for the proliferation rate as we wanted to keep the ratio r A /r A p ≪ 1 [41,42]. There are two reasons for this: firstly, averaged over 100 repeats of the ABM described above. We also plot the corresponding meanfield PDE solutions with r A p = 0.01 in red for species A and in green for species B.
a modelling choice to prevent agents clustering into proliferation-induced patches and secondly, it is biological realism that given the parameters of the model we could expect real biological cells will attempt proliferation events less frequently than movement events. We see good agreement between the two profiles for non-zero swapping probability. However, when the swapping probability is set to 0, we see discrepancies arising that are amplified as time increases. Recall that the disparity between the PDE and ABM profiles can be also observed in Figure 3 due to spatial correlations that are not accounted for by the mean-field PDEs. The addition of proliferation into the model increases the spatial correlations between site occupancies, leading to greater disparity. One way to partially rectify this problem is by modelling the higher order moments in the PDE description [33][34][35] however for the purposes of this work, we simply note that allowing swapping breaks up the correlations more effectively than in its absence leading to better agreement between the ABM and the population-level densities.

B. Swapping model with cell-cell adhesion
Another interesting application of swapping is the formation of patterns in densely crowded environments. In this section, we use a cell-cell adhesion model with swapping to investigate how biologically plausible patterns can form starting from a randomly seeded domain. Our model is based on the similar cell-cell adhesion model studied previously by [20,21,23] who consider adhesion between identical agents. Here, we extend the model to incorporate two types of agents with swapping to facilitate the movement events.
For the purpose of this paper, we assume adhesion between two species, A and B, on a fully with its neighbours is given by, However, if the focal agent is a type-B agent then, Here, A z are binary taking a value of unity if the site with position z ∈ Z ij is occupied by a type-A agent or 0 otherwise. Therefore, z∈Z A z is the sum of occupancies of the sites in Z ij that are occupied by type-A agents. Likewise, z∈Z B z is the sum of occupancies of the sites in Z ij that are occupied by type-B agents. Similarly, for a type-A agent occupying the target site the probability that it breaks existing connections with its neighbours is given by, However, if the agent occupying the target site is a type-B agent then, Here, the set Y z contains the positions of sites in the neighbourhood of the target site, z.
Therefore, y∈Yz A y is the sum of occupancies of all the sites in Y z that are occupied by a type-A agent and y∈Yz B y denotes the sum of occupancies of all the sites in Y z that are occupied by a type-B agent.
The probability of a successful swap given a movement event has been attempted is therefore a product of the swapping probability, the probability of the focal agent at site (i, j) breaking links with its neighbours in order to move out and the probability of the agent at the target site breaking links with its neighbours in order to move in, i.e., Here we are considering that the link between the two swapping agents has to be broken in (c) Figure 11. A schematic illustrating swapping in the cell-cell adhesion model. A red site represents a type-A agent and a green site represent a type-B agent. Agent 1 at site (i, j) attempts to swap with agent 2 at site (i + 1, j). We have coloured in black the neighbouring sites that are unimportant in this context (i.e. occupancy of these sites does not affect the probability of a successful swap happening between agent 1 and agent 2). The initial configuration of the lattice is shown in (a). Agent 1 attempts to swap with agent 2 (b) and the state of the lattice after the successful swap is shown in (c).
order for the swap to take place. As an example, consider the case in which the focal agent at site (i, j) attempts to jump into the neighbouring site (i + 1, j) ( Figure 11). The agent occupying the focal site is a type-A agent, shown in red and labelled as 1 and the target site (i + 1, j) is occupied by a type-B agent, shown in green and labelled as 2. Since the focal agent is red, it is adhesive to other red agents in its Von Neumann neighbourhood.
There is only one site (with position (i, j − 1)) neighbouring agent 1 that is occupied by a red agent and therefore, z∈Z ij A z = 1. Since agent 2 (target site) is a green agent, it is adhesive to the other green agents in its Von Neumann neighbourhood. There is only one green agent in the neighbourhood (position (i + 2, j)) and hence, y∈Y i+1,j B y = 1.
For an arbitrary ρ, p and q, the probability of swap in this situation given a neighbouring site to move into has already been chosen can be written as In Figure 12, we present some results that demonstrate the importance and impact of swapping in a model of adhesion-mediated pattern formation. For this, we consider a square lattice with L = 20 sites in both the horizontal and vertical direction. As before, one compartment can accommodate no more than a single agent at a time. We impose periodic boundary conditions. We seed the lattice with type-A and type-B agents at a density of 0.5 each, where the initial positions of the agents on the lattice are assigned uniformly at random (Figure 12(a)). We let positions of the agents evolve according to the kinetics described above using the movement rate r A = r B = 1 and ρ = 1.
Snapshots of the evolving lattice occupancy at t = 15, 000 for two different sets of p and q values are shown in Figure 12 even amongst some densely packed configurations. We considered a two-species system and allowed it to evolve according to the dynamics of our model. We found that swapping enhances the movement of agents by allowing agents to mix more compared to the pure volume-exclusion model. Comparing the ABM to the population-level model we found excellent agreement between the two as long as the swapping probability was sufficiently large.
To understand how swapping affects agent movement at an individual level, we analysed simulated agent tracks to determine the time-uncorrelated individual-level diffusion coefficient. We found that swapping enhanced the movement of agents in all cases compared to the volume-exclusion model. Using the probability master equation, we were able to analytically derive an expression for the diffusion coefficient that confirms the relationship obtained via the simulated tracks.
In Section IV, we demonstrated the importance of swapping via a couple of examples. In the first application, we considered a cell migration model with proliferation. We found that swapping accelerates the proliferation process by allowing the agents to disperse more and breaking spatial correlations. We found that the time to reach the domain's carrying capacity varies inversely with the swapping probability. Deriving the PLM and comparing it to the average density of the ABM, we showed that there is a good agreement between the two profiles for sufficiently large swapping probability. For the ρ = 0 case, we note discrepancies that are caused by the build up of spatial correlations between sites.
For the second example, by incorporating swapping into a cell-cell adhesion model, we showed that agents can spontaneously rearrange themselves into clusters to form patterns even on densely populated domains. We stress that without swapping these patterns would not be realised. Biological patterns involving cell-cell adhesion have been observed amply in experimental work. Honda et al. [43] reported chequered patterns in Japanese quail oviduct epithelium which were later modelled mathematically by Glazier and Graner [14] using a cellular Potts model. Armstrong [44] observed sorting of neural and pigmented retinal epithelial cells in chick embryo where neural cells completely engulfed the pigmented cells from an initially disordered arrangement. Engulfment was modelled mathematically using a cellular Potts model [14,15] and by Armstrong et al. [45] using a continuous approach.
Models incorporating swapping may prove useful when analysing multi-species systems in which cells need to migrate in a crowded epithelium or tissue parenchyma. These systems are often interrogated experimentally using the generation of chimaeras or mosaics, analysing the resulting patterns generated by cell mixing [46][47][48][49][50]. For example the generation of corneal epithelial stripes in X-inactivation mosaic female mice hemizygous for an X-linked copy of the LacZ gene, expressing the enzyme β-galactosidase [51]. Here, a mixture of β-gal positive and β-gal negative limbal epithelial cells are specicifed as limbal stem cells whose progeny then migrate centripetally to generate β-gal positive and β-gal negative corneal epithelial stripes [52]. Experiments suggest that the degree of mixing of the limbal progenitors will generate stripes of different width composed of multiple clones of the same colour [51]. Swapping would be useful in understanding how cell mixing at the limbus impacts the final stripe pattern generated.
In conclusion, motivating our study by real-life examples, we have developed a cell migration model that incorporates swapping as a viable movement process. As well as adding biological realism, our model has the added benefit of better agreement between the corresponding continuum description and the ABM compared to the classical volumeexclusion process. We also saw that the ABM with swapping and cell-cell adhesion, when applied to cells in densely crowded environments, leads to pattern formation. We once again stress that the patterns would be unattainable under the traditional volumeexclusion model. The results in this paper hint that swapping is an important and overlooked mechanism in the context of modelling real biological scenarios and merits further explorations in conjunction with experimental data.

APPENDICES Appendix A: Single-species PDE vs. ABM results
Below we summarise the single-species swapping model and present results comparing the discrete and continuum models.
The single-speices model is a simplified version of the two-species model in the sense that there is only species on the domain of interest at density c and hence only a single movement rate. We summarise the the model here.
We let r be the rate of movement of an agent such that rδt is the probability that the agent attempts to move during a finite time interval of duration δt. The agent attempts to move into one of its four neighbouring sites with equal probability. If the chosen neighbouring site is empty, the focal agent successfully moves and its position is updated. If another agent already occupies the site, the agent at site (i, j) attempts to swap positions with the neighbouring agent with probability ρ. If the swap is successful, the two agents exchange positions with each other. Otherwise, the move is aborted.
In Figure 13 we present snapshots of the lattice occupancy for the single-species model at t = 0, 100, 1000 and for swapping probabilities ρ = 0, 0.5, 1 with r = 1 with reflective boundary conditions. We see that swapping seems to have no effect on the dispersion of the agents at a macroscopic scale as the density profiles are indistinguishable regardless of the swapping probability. If two agents are identical to each other and unlabelled then exchanging positions by swapping is equivalent to an aborted movement attempt in the volume-exclusion process and produces no change in the state of the system.
In the next section, we derive the macroscopic PDE describing the evolution of the mean lattice occupancy in the single-species model.

Single-species continuum model
Let C k ij (t) be the occupancy of site (i, j) on the kth repeat at time t, where C k ij (t) = 1 if the site (i, j) is occupied and C k ij (t) = 0 otherwise. The average occupancy of site (i, j) at time t after K runs is then given by, Let C ij (t) = C ij for conciseness. By considering the possible movement events of the agent at the site (i, j) during the small time step δt [24,25], we can write down the master equation for the occupancy of the site at time t + δt, The Here, is the diffusion coefficient given that ∆ 2 /δt is held constant in the diffusive limit. Equation (A3) describes the evolution of the lattice occupancy over time. It is well-known that for the simple exclusion process, the occupancy is described by the diffusion equation [19,31]. It makes sense therefore that in Figure 13 swapping made no difference to the overall occupancy of the lattice. In Figure 14 we compare the average column density of the ABM which is given by, to the solution of the one-dimensional analogue of Equation (A3) with reflective boundary conditions by averaging the PDEs over the y direction.
As expected, we see an excellent agreement between the two density profiles. We also see that there is perfect agreement even in the zero-swapping case which we did not see in the two-species example in Section II B. This is because in a single-species system there is no way of distinguishing between a successful swapping of the position of a pair of neighbouring agents and an aborted movement event due to the volume-exclusion principle. This leads to identical profiles for different ρ regardless of its magnitude.
In contrast to the two-species case, here there is perfect agreement between the ABM and the PDE profiles since in the single-species case the agents are indistinguishable from each other. In the two-species scenario, because the agents are competing for space, in the no-swapping situation one species affects the occupancy of the other, leading to spatial correlation and discrepancy between the ABM and PDE descriptions. Swapping serves to break up the correlation and improve agreement between the discrete and continuum models.

Appendix B: Two-species individual-level diffusion coefficients
Let P A ij (t) = P ij be the probability that a type-A agent occupies the site (i, j) and let P B ij (t) = P B ij be the equivalent for species B. We can write down the master equation for species A and species B at time t + δt where δt is a small change in time, We explain the meaning of the terms in the discrete-time master Equation (B1). Equation (B2) can be interpreted similarly. The equations describe the evolution of the probability that a focal agent of type A is occupying site (i, j) by considering possible movement events of the focal agent at site (i, j) and agents at neighbouring sites. Firstly, the terms which correspond movements that increase the probability that the focal agent sits at 5. Finally, the fifth term in the square brackets corresponds to a type-B agent at a neighbouring site successfully swapping positions with the type-A focal agent occupying site (i, j).
Subtracting P A ij from both side and dividing by δt and taking the limit as δt → 0 gives, After using the definitions in Equation (12) and (13), and suitably transforming the indices i and j, after simplifying it can be shown that the equations describing the evolution of the second moment of the position i and j of species A and B are given by, Given that ⟨i(0)⟩ = ⟨j(0)⟩ = 0, Thus, are the individual-level time-uncorrelated diffusion coefficients of species A and B, respectively.