Quantification of spatial and phenotypic heterogeneity in an agent-based model of tumour-macrophage interactions

We introduce a new spatial statistic, the weighted pair correlation function (wPCF). The wPCF extends the existing pair correlation function (PCF) and cross-PCF to describe spatial relationships between points marked with combinations of discrete and continuous labels. We validate its use through application to a new agent-based model (ABM) which simulates interactions between macrophages and tumour cells. These interactions are influenced by the spatial positions of the cells and by macrophage phenotype, a continuous variable that ranges from anti-tumour to pro-tumour. By varying model parameters that regulate macrophage phenotype, we show that the ABM exhibits behaviours which resemble the ‘three Es of cancer immunoediting’: Equilibrium, Escape, and Elimination. We use the wPCF to analyse synthetic images generated by the ABM. We show that the wPCF generates a ‘human readable’ statistical summary of where macrophages with different phenotypes are located relative to both blood vessels and tumour cells. We also define a distinct ‘PCF signature’ that characterises each of the three Es of immunoediting, by combining wPCF measurements with the cross-PCF describing interactions between vessels and tumour cells. By applying dimension reduction techniques to this signature, we identify its key features and train a support vector machine classifier to distinguish between simulation outputs based on their PCF signature. This proof-of-concept study shows how multiple spatial statistics can be combined to analyse the complex spatial features that the ABM generates, and to partition them into interpretable groups. The intricate spatial features produced by the ABM are similar to those generated by state-of-the-art multiplex imaging techniques which distinguish the spatial distribution and intensity of multiple biomarkers in biological tissue regions. Applying methods such as the wPCF to multiplex imaging data would exploit the continuous variation in biomarker intensities and generate more detailed characterisation of the spatial and phenotypic heterogeneity in tissue samples.

Oxygen (ω) Oxygen is supplied by blood vessels which are represented by static point sources, and is consumed by tumour cells and stromal cells as it diffuses through the domain. We assume that the oxygen concentration at each blood vessel is constant, ω 0 .
Since the timescale of diffusion for an oxygen molecule is faster than the timesteps used in our model to describe cell movement, we assume that the distribution of oxygen can be approximated using a steady state solution [1,2]. For simplicity, we rescale the oxygen concentrations with a factor of ω 0 , so that ω = 1 at each blood vessel. The governing equation is then given by for x ∈ Ω, where x i is the location of stromal or tumour cell i, the parameter D ω is the assumed constant diffusion coefficient of oxygen, κ ω is the oxygen consumption rate and δ is the delta function (δ(x) = 1 when x = 0, δ(x) = 0 otherwise). We impose Neumann boundary conditions ( ∂ω ∂x = 0) on domain boundaries, and assume that initially ω(x, t = 0) = 1 for all x that do not represent a blood vessel.
CXCL12 (C-X-C motif chemokine 12, ξ) CXCL12 is produced by perivascular cancer-associated fibroblasts (CAFs) and acts as a diffusible chemoattractant for M 2 macrophages [6]. We suppose that CAFs localise close to blood vessels, and, hence, make the simplifying assumption that blood vessels act as constant sources of CXCL12. We therefore assume that the concentration of CXCL12, ξ, is maintained at a fixed value ξ = ξ 0 at all blood vessels. We assume that CXCL12 decays naturally at a constant rate, λ ξ . The distribution of CXCL12 in the domain is then given by 0 = D ξ ∇ 2 ξ − λ ξ (2) where D ξ is the diffusion coefficient for CXCL12. We impose Neumann boundary conditions ( ∂ξ ∂x = 0) on domain boundaries, and assume that initially ξ(x, t = 0) = 0 for all x that do not represent a blood vessel.
CSF-1 (Colony stimulating factor-1, c) The macrophage chemoattractant CSF-1 is produced by tumour cells. It acts with EGF as part of a paracrine loop [7,8] and stimulates macrophage extravasation [6]. We assume that CSF-1 is produced at a constant rate, κ c , by each tumour cell, and decays at a constant rate, λ c . The distribution of CSF-1 is therefore described by the equation where D c is the diffusion coefficient of CSF-1 and the sum runs over all tumour cells i. We impose Neumann boundary conditions ( ∂c ∂x = 0) on domain boundaries, and assume that initially c(x, t = 0) = 0 for all x.
TGF-β (Transforming growth factor-β, g) TGF-β causes macrophages to alter their phenotype, and is produced by tumour cells. We suppose that it decays at a constant rate, λ g , and is produced by tumour cells at a constant rate, κ g . Combining these effects, we arrive at the following RDE for g(x, t): where D g is the diffusion coefficient of TGF-β and the sum is over all tumour cells, i. As for CSF-1, we impose Neumann boundary conditions ( ∂g ∂x = 0 on domain boundaries), and assume that g(x, t = 0) = 0 ∀x ∈ Ω.
EGF (Epidermal growth factor, ϵ) EGF is a diffusible chemoattractant for tumour cells that is produced by tumour associated macrophages. We assume that EGF is produced by macrophages at a rate that is linearly dependent on their phenotype p, so that macrophages with an extreme M 1 phenotype (p = 0) produce no EGF and macrophages with an extreme M 2 phenotype (p = 1) produce κ ϵ . We assume that EGF naturally decays at a constant rate, λ ϵ . Combining these effects, we obtain the following RDE for EGF, ϵ(x, t): where D ϵ is the constant diffusion coefficient of EGF, κ ϵ is the maximum rate of production of EGF, the sum is over all macrophages i, and p i ∈ [0, 1] is the phenotype of macrophage i and is defined below. As for CSF-1 and TGF-β, we impose Neumann boundary conditions ( ∂ϵ ∂x = 0 on the boundaries of the domain Ω), and assume that ϵ(x, t = 0) = 0 ∀x ∈ Ω.

Agents
Our model distinguishes four types of agents: macrophages, tumour cells, stromal cells and necrotic cells. A fixed number of static blood vessels are also distributed throughout the domain. The location of each agent is represented by its cell centre.
We use Newton's second law to determine the equations of motion for macrophages, tumour cells, stromal cells and necrotic cells. Neglecting inertial effects in the overdamped limit, the force balance for cell i can be written as: where ν is the drag coefficient (assuming that the drag on a cell is proportional to its velocity), and F i denotes the net force acting on cell i. The forces that act on each cell type are indicated in Fig 2 in the main text. All cells are subject to mechanical forces due to cell-cell interactions (incorporating repulsion due to volume exclusion and attraction due to intercellular adhesion). Here we use the overlapping spheres approach [9,10], assuming that two cells interact if the distance between their cell centres is less than a fixed distance R int . Specifically, for cells at locations x i and x j , if |x i − x j | < R int then the interaction force between cells i and j is parallel to the vector (x i − x j ) connecting their centres. The resting spring length between the cell centres, s i,j , is the sum of the equilibrium spring lengths associated with each cell (s i,j = s i + s j ). For most cells i, the resting spring length is equal to the approximate radius of a cell (s i = R Cell ≡ 0.5 in non-dimensional units). For newly divided and necrotic cells, s i changes over time to account for cell growth and shrinkage (see below).
The net mechanical force acting on cell i is the sum of all interaction forces due to cells within its interaction radius: where F m i,j is the mechanical force between cells i and j. This force always points in the direction of the vector connecting the cell centres and has magnitude: where x = |x i − x j | − s i,j is the overlap between cells i and j, µ is a parameter describing the spring stiffness and α determines the strength of intercellular adhesion between neighbouring cells. We normalise lengths with the lengthscale R Cell , assuming that 1 cell diameter = 2R Cell = 20µm. Following cell division, the radius of the daughter cells is set to s i = R Cell 2 and increases linearly over one hour until s i = R Cell . For necrotic cells, s i decreases linearly to zero overτ i hours (see below), and then the cell is removed from the simulation. The associated spring stiffness µ of springs attached to a necrotic cell is reduced linearly at the same rate.

Macrophages
Extravasation Macrophages extravasate from blood vessels in response to local levels of CSF-1. At each timestep, the probability that a macrophage extravasating in a one hour time period is P ex , where: where the parameter P ⋆ controls the maximum possible probability of macrophage extravasation per hour from each vessel, and c 1/2 is the concentration of CSF-1 at which this probability is half-maximal.
Phenotype The diffusible species included in our model bind to receptors on the outer membranes of macrophages and tumour cells and regulate behaviours including chemotaxis and the production of other species. While we do not explicitly model surface receptors, the receptor CXCR4 (C-X-C motif chemokine receptor type 4) plays an important role in this system: when exposed to TGF-β, macrophages increase expression of CXCR4 and become sensitive to gradients of CXCL12 [6]. We model this by associating each macrophage with a phenotype, p ∈ [0, 1], which determines the extent to which it exhibits M 1 -or M 2 behaviour. When a macrophage first enters the spatial domain, it has phenotype p = 0. Exposure to TGF-β causes its phenotype to increase irreversibly. The rate of change of the phenotype of macrophage i at location , the non-negative parameter ∆p determines the rate at which phenotype changes, and g crit is a critical TGF-β threshold above which macrophage phenotype increases. When p = 1 the macrophage has a fully M 2 phenotype, and p no longer changes. Phenotype affects the following macrophage behaviours: • the rate at which they produce EGF; • the rate at which they kill tumour cells; • their chemotactic sensitivity to spatial gradients of CSF-1 and CXCL12.
Equation (5) describes how the rate of EGF production depends on p; we now describe how cell killing and phenotype-dependent chemotaxis are implemented.
Cell killing Macrophage phenotype determines the probability that a macrophage kills a tumour cell on a given timestep. The probability that macrophage i of phenotype p i kills a tumour cell in one hour is where P ⋆ φ is the maximum probability that a macrophage kills a tumour cell in one hour. t φ is the time since the macrophage last killed a tumour cell, and the parameter t cool defines a 'cooldown' period during which the macrophage cannot kill another tumour cell. On each timestep the probability of each macrophage attempting cell killing is evaluated. If a macrophage attempts killing, any cells within distance 1 of the macrophage are identified and, if there is at least one tumour cell, one of the tumour cells is selected at random to be killed. The tumour cell is labelled as necrotic, and t φ is reset to 0 for that macrophage.
Force laws for macrophages In addition to mechanical forces, macrophages are subject to forces which account for random movement and chemotactic forces due to spatial gradients of CXCL12 and CSF-1. We account for random movement via the force where the coefficient D describes the strength of the random force and n = (n x , n y ), where n x and n y are random variables drawn from a standard normal distribution. The net force due to chemotaxis acting on macrophage i depends on its phenotype, p i , as follows: where the parameters χ m c and χ m ξ control macrophage sensitivity to spatial gradients of CSF-1 and CXCL12 respectively. We suppose that macrophage sensitivity to CSF-1 and CXCL12 gradients scales linearly with phenotype.
The net force acting on macrophage i in Equation (6) is: Tumour cells Cell cycle and proliferation Each tumour cell has an internal cell cycle which advances at a rate that depends on the local oxygen concentration and two oxygen thresholds, ω tum H < 1 and ω tum N ≤ ω tum H . When ω tum N < ω ≤ ω tum H , the cell becomes hypoxic and pauses its cell cycle until ω increases above this threshold. If 0 ≤ ω ≤ ω tum N , then the tumour cell dies and becomes a necrotic cell. We also account for contact inhibition in our model. If a cell's area falls below a proportion A H i tum of its target area, then its cell cycle pauses until space is available for proliferation. Each tumour cell i has a subcellular variable denoted T i which tracks its progress through the cell cycle cell as follows: where A i is the area of the cell calculated as A i = πr 2 i and r i is the estimated cell radius calculated via the average separation between cells within the interaction radius of i.
Each cell has a target cell cycle duration τ i , drawn at random from a uniform distribution of U (0.75τ tum , 1.25τ tum ) with average cell cycle duration τ tum . When T i = τ i , cell division occurs and a new cell is placed at a distance half a cell diameter away, in a randomly chosen direction. Both cells are assigned new cell cycle durations. T i is set to 0 for each cell, and then evolves according to Eq (15). The equilibrium spring length associated with the new cells is set to s i = 0.5R Cell to account for their reduced size, and increases linearly over the course of 1 hour until it reaches s i = R Cell .
Force laws for tumour cells Tumour cell movement is governed by the force balance described in Equation (6). The net force incorporates intercellular interactions as described above, and a chemotactic force due to spatial gradients in EGF. The chemotactic force, denoted F χϵ i , has the form where the parameter χ T ϵ determines tumour cell sensitivity to EGF gradients. The net force used in Equation (6) is therefore

Stromal cells
Cell cycle and proliferation The same cell cycle model is used for stromal and tumour cells, but it is parameterised differently to account for the increased ability of tumour cells to survive in adverse conditions such as lower oxygen environments or increased mechanical pressure from neighbouring cells. Like tumour cells, stromal cells possess a subcellular age variable T i which evolves according to the equation where the parameter ω str H determines the oxygen threshold below which stromal cells become hypoxic and the parameter A H i str defines the proportion of a stromal cell's target area below which the cell cycle stops.
As with tumour cells, we define a second oxygen threshold ω str N , below which stromal cells become necrotic. We suppose that ω str H < ω tum H and A H i str < A H i tum to account for the ability of tumour cells to proliferate in more adverse environments than stromal cells.
Force laws Stromal cells are subject to the same intercellular forces as macrophages and tumour cells, and defined by Equations (7)- (8). The force balance for stromal cells used in Equation (6) is

Necrotic cells
When a stromal or tumour cell is marked for cell death, as a result of oxygen starvation or being killed by a macrophage, it irreversibly becomes necrotic. A necrotic cell i occupies space forτ i hours, whereτ i is drawn from a uniform distribution U (0.75τ , 1.25τ ) andτ is the average duration of necrosis. Over this time period, the necrotic cell shrinks in size, by reducing its equilibrium spring constant, s i , at a constant rate until s i = 0. The cell is then removed from the simulation. While s i is being reduced, the spring constant µ associated with the cell is also reduced to 0 at a constant rate, to account for weakening of the intercellular forces between degrading necrotic cells and other cells.

Initial and boundary conditions for cells
RDEs describing the diffusible species are solved numerically on a regular triangular mesh with edge length 1 cell diameter. We initialise the simulation by selecting lattice sites to act as point vessels, which do not occupy space or interact directly with cells in our model. All lattice sites more than R B cell diameters from the centre of the domain are possible locations for blood vessels, and N B of these sites are chosen, at random, to be blood vessels. This ensures that the centre of the domain is at least R B cell diameters from the nearest blood vessel and, hence, that there is sufficient space to observe macrophage movement between blood vessels and the tumour, while also ensuring that cells near the domain boundaries are well-oxygenated. The domain is initialised with stromal cells filling the domain in rows that are 0.75 cell diameters apart, with alternating rows offset by 0.375 cell diameters (i.e., forming a hexagonal lattice). We place four tumour cells in a cluster at the centre of the domain approximately 0.5 cell diameters apart. Stromal and tumour cells are assigned cell cycle durations τ i from the relevant distributions, and cell cycle progression times T i selected at random from a uniform distribution U (0, τ i ). All cells are constrained to remain within the domain by imposing reflective boundary conditions.

Schematic
The schematic in Fig S0 shows

Table of ABM parameters
In Table 1 we list the model parameters, their default dimensional and dimensionless values or ranges, and supporting references where these are available. Values of some parameters, indicated with * , have been estimated based on model behaviour.
Our model is implemented such that typical scales are given by: • Length: 1 cell diameter (taken as 20µm) is 1 unit of length.
• Time: 1 hour is 1 unit of time.
• Concentration: the boundary concentration of oxygen is 1.