Linking discrete and continuous models of cell birth and migration

Headnote

Self-organization of individuals within large collectives occurs throughout biology. Mathematical models can help elucidate the individual-level mechanisms behind these dynamics, but analytical tractability often comes at the cost of biological intuition. Discrete models provide straightforward interpretations by tracking each individual yet can be computationally expensive. Alternatively, continuous models supply a large-scale perspective by representing the 'effective' dynamics of infinite agents, but their results are often difficult to translate into experimentally relevant insights. We address this challenge by quantitatively linking spatio-temporal dynamics of continuous models and individual-based data in settings with biologically realistic, time-varying cell numbers. Specifically, we introduce and fit scaling parameters in continuous models to account for discrepancies that can arise from low cell numbers and localized interactions. We illustrate our approach on an example motivated by zebrafish-skin pattern formation, in which we create a continuous framework describing the movement and proliferation of a single cell population by upscaling rules from a discrete model. Our resulting continuous models accurately depict ensemble average agent-based solutions when migration or proliferation act alone. Interestingly, the same parameters are not optimal when both processes act simultaneously, highlighting a rich difference in how combining migration and proliferation affects discrete and continuous dynamics.

Received: 22 December 2023 Accepted: 10 May 2024

Subject Category:

Mathematics

Subject Areas:

applied mathematics/mathematical modelling/ theoretical biology

Keywords:

non-local interactions, zebrafish, self-organization, aggregation equations

1. Introduction

Self-organization of individual agents is a key feature of life. It occurs ubiquitously throughout the natural world, from the macroscopic example of bird flocking [1-4] to the microscopic phenomenon of cell sorting during development [5-9]. The degree to which members of a group coordinate their movement, proliferation and competition accounts for pattern diversity across biological scales. Alongside experimental approaches, mathematical models can help identify the underlying behaviours that give rise to specific collective dynamics. However, a trade-off often exists between tractability and detail when building models of pattern formation, due in part to the multiscale nature of biological systems. Consequently, better quantitative characterization of the relationship between analytically tractable models and more biologically representative approaches will improve our understanding of self-organization throughout nature.

Here, we help address this open challenge using pigment cell dynamics in zebrafish patterns as a motivation. The zebrafish (Danio rerio) is a popular model organism for studying pattern formation, as dark stripes and gold interstripes emerge in its skin during development [10-13]. As we show in figure 1, these stripes result from the coordination of interactions among several types of cells, including black melanophores and gold (dense) xanthophores [10,16-20]. Experiments that perturb stripes-e.g. by laser ablation [16,21]-demonstrate how cell-cell signalling and external cues contribute to the creation of alternative motifs such as spots or labyrinths. A rich diversity of mutant patterns, including widened or curvy stripes, also emerge when cell interactions are altered due to genetic mutations [19,22].

Data-driven mathematical models can help uncover the drivers of zebrafish pattern formation and other biological phenomena exhibiting self-organization by identifying important phase transitions, isolating the effects of specific processes such as cell division, and providing hypotheses that can guide the design of in vivo experiments [23-26]. Different modelling frameworks yield insight at the population or individual level, depending on how they represent members of a group. One modelling approach involves tracking how the position of each individual changes in time. These so-called 'discrete' systems include centre-based models [27,28], cellular automata [29,30], cellular Potts models [31,32] and vertex models [33,34]. Within the setting of zebrafish patterning, agent-based models (ABMs) have been developed that restrict cells to occupy certain locations 'on-lattice' [35-37] or allow them to roam freely, 'off-lattice', in the domain [14,38-40]. Due to their ability to work on the same length scales as empirical data, ABMs provide an intuitive connection to experiments and allow for detailed predictions about how interactions between agents drive group behaviours. However, ABMs can be prohibitive to simulate when the number of individuals is large, and understanding their long-time behaviour under alternative rules and parameters relies on extensive computation [41].

A second modelling approach uses continuous functions to represent the 'average' density of agents in a collective, with their dynamics governed by a partial differential equation (PDE) in space and time. Continuous models, including reaction-diffusion equations, Boltzmann-like kinetic equations and integro-differential equations (IDEs), typically cannot resolve individuals and, instead, track the ensemble average (EA) behaviour of a population. However, these models are more amenable to mathematical analysis and more readily provide insight into long-term behaviour than discrete frameworks do [42,43]. For example, changes in patterning may arise because of Turing-like instabilities [44-46] or due to alterations in physically based interactions such as cell-cell adhesion [7,38,47-50]. In the case of zebrafish patterns, researchers have applied a wide swath of continuous models-including reaction-diffusion equations [16,21,35,51,52] and non-local PDEs [7,53-55]-to better understand cell dynamics.

Despite the differences between discrete and continuous approaches, it is possible to establish a mathematical link between these representations in the limit of infinite individuals. This procedure, known as 'coarse-graining', derives differential equations from a given discrete model and yields information about its EA behaviour [56-62]. For example, the authors in [63-65] derive logistic IDEs from stochastic processes that describe the birth and death of individuals undergoing Darwinian evolution in the limit of large numbers. Coarse-grained descriptions become inaccurate when relatively few individuals are present, however, as is the case in many biological contexts such as pattern formation in zebrafish. Many approaches also neglect potentially important spatial correlations between cells-caused, for instance, by division or competition-that may play a critical role in pattern dynamics [52,66-71]. While it is possible to go beyond this 'mean-field' setting by deriving continuous models that respect higher-order correlations [72-76] and hard-core interactions [77], these methods still rely on simplifying assumptions that introduce errors between the discrete and continuous frameworks. Controlling these errors in a biological setting is an important objective.

We help tackle this problem by developing a pipeline to minimize spatio-temporal discrepancies between continuous models and individual-based data in settings with biologically relevant, dynamic cell numbers. The method relies on fitting parameters that effectively dilate the time variable of continuous-model solutions. Our approach can be used to describe biological self-organization in systems whose macroscopic description can be derived or inferred. We apply our approach to a case study motivated by stripe formation in zebrafish skin, as this allows us to illustrate its utility and interpretablity for experimentally measurable quantities with biologically meaningful spatial and temporal units. In §2, we describe the ABM that we simulate to generate synthetic individual-based data associated with self-organizing phenomena. In its full form, the ABM motivating our work [14] admits pattern formation via non-local rules for cell birth, death and movement that are inspired by the underlying biology of zebrafish-skin patterns (figure 1) [14,40]. To focus on the presentation of our pipeline, however, we simplify some biological complexity by reducing this model [14] to focus on a single cell type-melanophores or dense xanthophores. (We plan to extend this pipeline to multiple cell types in future work.) We then detail the corresponding continuous descriptions and our method for matching their solutions to EA ABM data even in scenarios with finite and time-dynamic cell numbers. We present our results for black melanophores in §3, and-as a means of demonstrating the generality of our methodology-apply the same approach to dense xanthophores in electronic supplementary material, S6.

2. Mathematical models and methods

In §2.1, we develop our ABM for cell migration and derive its continuous counterpart. Subsequently, in §2.2, we introduce our discrete model for cell birth and develop a corresponding continuous IDE model. We present our full models of migration and proliferation in §2.3. Lastly, we present our approach to estimating scaling parameters in our continuous models from EA ABM data in §2.4. Following previous ABMs [14,40] of pattern formation in zebrafish, we assume (i) migration is governed by conservative forces between pigment cells and (ii) non-local interactions inform cell birth in a two-dimensional (2D) plane [14,39,40]. Throughout this paper, we refer to:with the exception of figure 7 where we consider a one-dimensional (ID) domain; there M(x, t) is the number density of melanophores in cells mm-1. Because it appears several times, we define the indicator function l(Concirtion} nere as

2.1. Models of migration

Our ABM for cell movement tracks the position, M,(f), of each cell, indexed by ie {1,..., NM(f)}, at time t >0. The movement of each melanophore depends on its interactions with surrounding melanophores, following an overdamped version of Newton's second law. The forces are assumed to be conservative, i.e. they may be written as the gradient of a potential. This leads to the following system:

Here, Fv , is the net force arising from all cell-cell interactions according to the potential Wmm (r)/ where r is the inter-particle displacement. The potential can encode cell-cell repulsion and adhesion, depending on the sign of its gradient along the direction between two cell centres. Many choices for Wmm(i) are possible, including harmonic, power-law, Morse and Lennard-Jones potentials, among others. Here, we use exponential potentials given by in accordance with prior ABMs of pattern formation in zebrafish [14,40]; see figure \c and table 1 for parameter values and their biological interpretations. To model cells communicating through cellular extensions or dendrites [19,20], secreted signals [80] or cell-cell contact [81], we assume forces on cells are zero beyond some cut-off distance dmax. We represent this using the notation VWJjjjr) - VWmmW l{|r|<dm"}(r)/ and set NM(t) =NM when there is no cell birth. We remark that the model given by equation (2.2) is deterministic in the sense that it does not include Brownian motion.

The associated continuous model describes the melanophore density, M(x, t). Integrating M(x, t) over a bounded region yields the total number of melanophores within that area at time t. Following the coarse-graining procedure in [60,82,83], an outline of which is found in electronic supplementary material, §1.1, we obtain the PDE below:where the force VW^M is the same as in equation (2.2) and * is the convolution operator [14,40]. The parameter aMM in equation (2.4) is not inherent to the coarse-graining procedure; rather, it accounts for possible differences between the discrete and continuous models. Indeed, simulating equation (2.4) with the parameter values listed in table 1 and aMM - 1 (the value expected from the mean-field approximation) does not always capture the ABM dynamics; see figure 2.

The individual and EA ABM results demonstrate that cells disperse until they are about 55-115 |im apart at t = 150 days. The PDE with %m - 1/ however, predicts that cells travel about 250 |im further in the same time period. Additionally, the PDE cell density is lower than the EA ABM density near the centre of the domain, implying that cells are more separated there. The continuous solution at earlier times more closely resembles the EA ABM result at t = 150 days, however, which suggests that PDE solutions evolve at a faster time scale than that of the discrete model. The parameter aMM effectively dilates the time variable, such that solutions travel "mm times more quickly. Thus, a non-unitary value of "mm is likely to produce a better match between the discrete and continuous solutions. To our knowledge, the value of aMM cannot be derived a priori. Instead, we develop an approach for estimating its value based on ABM data in §2.4, and will pursue an analytic derivation for this parameter in future work.

2.2. Models of cell birth

Our ABM for cell birth consists of stochastic, discrete-time rules which we adapt and simplify from [14] (motivation for these rules may be found in electronic supplementary material, SI .2). Specifically, at each time step (i.e. day) in a simulation, we select Nuv £= ^ locations uniformly at random from Q and evaluate them synchronously for possible cell birth. Each selected location, z, represents the position of a precursor cell that may differentiate into a melanophore based on the signals that it receives. The conditions for melanophore birth in the ABM [14] depend on both neighbouring melanophores and dense xanthophores, as we show in electronic supplementary material, figure SI. Since we restrict to one population in this paper, we simplify the rules from [14]; see electronic supplementary material, SI .2, for details. In particular, a new melanophore emerges at position z according to the rule

According to equation (2.5), new cells appear near existing melanophores until the maximum number of cells-namely, c^BM-in Hfoc, the interaction region between cells, is reached; see table 1 for parameter values. While equation (2.5) is deterministic, stochasticity enters our ABM through our Nbir randomly selected positions {z}. Similar stochastic rules can also be used to model cell death, as in [14,40], although we do not consider them here.

We do not know of existing methods for rigorously deriving continuous models of cell birth from off-lattice ABMs with this noise structure. We therefore adopt a phenomenological modelling approach, in which we create a continuous model whose governing equations mimic the stochastic interaction rules. We reason that the number density of cells in this setting must increase at a constant rate (proportional to Nbir) when continuous versions of the overcrowding and short-range activation restrictions are met, since this occurs at the individual level in the ABM. Furthermore, we represent the density restrictions with an indicator function using the integral of the number density over £2ioc as an argument, since the latter quantity yields the total number of cells within that region. This leads to the following model:

where c+ is the continuous equivalent of the density-limiting parameter c^BM in equation (2.5); Nbir has the same value as in our corresponding ABM; and yE U.+ is a parameter that effectively dilates the time variable in a similar way as "mm in the cell movement model. The units of y must be inversely proportional to those of the domain size in order to make the dimensions of equation (2.7) consistent. Its value is unknown; however, we can employ a phenomenological argument to determine an expected value by integrating equation (2.7) over the whole domain. This yields an upper bound on the number of cells born per unit time of yNbir \i2\, hence one expects /" IQI _1 to maintain a maximum rate of Nbir cells born per day as in our ABM. We note, however, that this argument does not take into account possible clustering or other spatial correlations that can occur in the discrete setting, which may change the values of y and c+ from their expected values. While we could address this by allowing both parameters to depend on the proportion of the domain in which the birth conditions are fulfilled, we leave this extension for a future study and simply estimate uniform values for y and c+ by fitting to EA ABM data, as we do for "mm in the movement-only model. We overview our approach for estimating the values of c+ and y in §2.4.

2.3. Full models of cell movement and birth

We combine our descriptions of cell movement and proliferation to form our full discrete and continuous models. For our full ABM, we move cells according to equation (2.2) and then introduce new agents based on equation (2.5) at each simulated day; see electronic supplementary material, S2, for details. For our continuous model, we combine the terms related to movement and birth, such that the cell density evolves according to where the parameters "mm, J, Nbir and c+ have the same interpretations as in §§2.1 and 2.2. Importantly, by assuming that these parameters have the same interpretations, we are assuming that migration and proliferation are additive, so that combining them has no extra influence. Our fitting approach for these parameters, discussed below, allows us to evaluate this choice and better understand the interplay of these two mechanisms in discrete and continuous settings.

2.4. Parameter estimation procedure

We identify the values of the parameters-"mm in equation (2.4) and {y, c+} in equation (2.7)-by minimizing the sum of squared differences (hereafter referred to as the 'L2 error') between the continuous and EA discrete solutions over time and space. Because we are able to model the isolated processes of cell birth and movement separately, or consider them acting simultaneously, there are two ways of estimating parameters: by fitting all three parameters simultaneously to data from the combined model, or by fitting them in a modular fashion by considering cell movement and birth in isolation from each other. For the remainder of this paper, we adopt a modular approach because it allows us to probe the particular effects of cell movement and birth in detail (see figure 3 for an overview), and we present a study of simultaneous estimation in electronic supplementary material, S3. In particular, by using our modular parameter values in the combined PDE model, we can investigate their interplay and better understand the additive effects of individual-level mechanisms on the accuracy of continuous models. As we discuss in electronic supplementary material, S3, our modular approach may also supply additional information that can improve parameter estimation.

For example, we show in electronic supplementary material, S3, that hmm and y are not uniquely identifiable if they are fit only to the combined EA ABM data, whereas figure 4b and figure 7/ suggest we can uniquely identify them with a modular approach. We overview our method for parameter estimation below; for parameter values, see table 1. We refer to electronic supplementary material, S2-S3, for further details about our implementation of the pipeline in addition to alternative choices that could be taken in parameter estimation (such as fitting to earlier times, less refined spatial data etc.).

We consider biologically meaningful time scales (i.e. days), length scales (i.e. mm), and cell densities and stress empirical units throughout our results. This choice supports future studies that may treat pattern formation with multiple cell types. Throughout our simulations, we consider a domain of size 3x3 mm2 (with one ID exception in figure 7). We implement four initial conditions to extract common features of cell interactions from different geometric scenarios. The first involves a square region of melanophores in the centre of the domain (Box), which mimics the symmetry of the domain. For the second initial condition, we place a single stripe of melanophores (Stripe), which is motivated by the typical patterning observed in wild-type zebrafish. For the third initial condition, we consider two rectangular regions of melanophores (Offset rectangles), which take into account non-standard geometries and the meeting of two disjoint melanophore populations. Finally, the last initial condition we consider involves two melanophore stripes (Two stripes), which explores the interactions between two disjoint melanophore populations with biologically realistic sharp fronts. (See electronic supplementary material, figure S6, for a summary of these initial conditions.) We initialize individual ABM simulations by sampling cell positions uniformly in these regions for each respective initial condition, and initialize our continuous models by setting the cell density uniformly equal to the estimated biological density of 400 cells mm-2 [14].

In Step 2 of our pipeline in figure 3, we solve our discrete models with an explicit approach. Specifically, we solve equation (2.2) with an explicit forward Euler scheme. To model differentiation from uniformly distributed precursor cells [14], we solve the birth-only ABM by selecting Nbir sites in the domain uniformly at random at a fixed time step (here, 1 day) and placing a new cell at each position that meets the conditions given by equation (2.5). (Following the approach in [14,39,40], we evaluate all Nbir locations for potential cell proliferation at the same time. This synchronous evaluation means that it is possible, though uncommon, for more than c^BM cells to be present in a local neighbourhood, and the choice of parameters in our model, based on the ABMs [14,39,40], accounts for this possibility.) We solve our ABM combining migration and birth by simulating equation (2.2) and then implementing cell birth as above, with migration evaluated using a shorter time step than the time step for birth events.

To compare ABM results directly with the cell density from our continuous models, we obtain an EA distribution by simulating many ABM realizations (here, between 103 and 104 simulations), sorting all the cell locations into a histogram of Nhist x Nhist voxels (or Nhist x 1 voxels in ID), and normalizing by the number of simulations and the voxel area for each day simulated; see Step 3 of our pipeline in figure 3. Other ways of relating ABM and PDE results are also possible, for example by introducing Gaussian kernels at each cell location [84]. However, as far as the construction of the histogram is concerned, we expect comparable results for particles and localized Gaussian kernels. Furthermore, we show in electronic supplementary material, S3, that histogram voxel size used to bin EA ABM data and the final time used to fit the continuous equations only play a minor role in affecting the parameter values we obtain, at least for the cell movement model.

As part of Step 4 of our pipeline in figure 3, we need to solve our continuous models, and we do so with explicit approaches. Specifically, we apply a first-order finite volume scheme for the migration model (equation (2.4)), a forward Euler method for the continuous cell birth model (equation (2.7)), and a combined finite volume/forward Euler scheme for the full continuous framework (equation (2.8)). (More details about the particular time steps used to simulate the discrete and continuous models are in electronic supplementary material, §2.) We simulate the continuous models on an Nbin x Nbin mesh and, to match with EA ABM solutions, record the average cell density at each day on a (possibly coarser) grid of Nhist x Nhist voxels.

We compute continuous model parameters by minimizing the L2 error between the continuous and EA ABM results across time. Notably, this nonlinear least-squares problem is equivalent to maximum likelihood parameter estimation when the densities produced from the ABM simulations are independent, identically distributed normal random variables with constant variance and mean equal to the continuous solution. The L2 error that we minimize is given by where Afrecord denotes the time steps at which data are collected; Ax and Ay are the spatial step sizes of the histogram used to compare the EA ABM and PDE data; M,"sl- ¦ is the continuous-model solution at time tn and position (x;, i/y); and M^gMl- is the corresponding EA ABM result. For the birth-only model, we consider fitting to either the L error as above or simply the difference in the total cell count of the two datasets (we verify in electronic supplementary material, S2, table S4, that fitting to the L2 error produces similar parameter estimates). When we consider cell birth, we simulate our models with different values of Nbir and estimate parameters by minimizing the sum of the errors across these Nbir values. We fit parameters related to cell proliferation sequentially-that is, we determine the optimal value for c+ before y. We verify in ID that sequential and simultaneous estimation does not lead to significant difference in parameter values; see electronic supplementary material, S3.

3. Results

We now present our results linking discrete and continuous models of cell migration (§2.1), birth (§2.2), and migration and birth (§2.3). We first isolate each interaction process, separately identifying the values of "mm in equation (2.4) and {y, c+} in equation (2.7). As we note in §2.4, this choice allows us to extract the distinct effects of each mechanism. We then determine how this simplification affects the ability of the full continuous framework, given by equation (2.8), to approximate EA ABM solutions. By considering different initial conditions (we discuss the details and motivations for these in electronic supplementary material, S4), we demonstrate the robustness of our fitting procedure. Our results show how the time scales of proliferation and movement in our continuous model may depend on numerical implementation and the frequency of stochastic cell birth controlled by Nbir. Moreover, our modular fitting approach highlights important considerations to account for in more general systems where agents are moving and changing in number.

3.1. Cell migration

We estimate "mm/ the scaling parameter that controls the dynamics of melanophore movement. Figure 4" presents the values of "mm that minimize the L2 error between the continuous solution of equation (2.4) and EA ABM results for our four initial conditions (see §2.4 and electronic supplementary material, S2S4). In each case, the optimal value of "mm is positively correlated with our PDE mesh resolution, i.e. greater values of "mm are associated with larger Nb;n - N^t values. This unitless parameter appears to converge to around 0.60-0.66 as the mesh resolution increases. There is at most a 2.5% relative difference between the values of "mm that we find when Nbm - 240 versus when Nb;n - 480 for our Box initial condition. These results suggest that "mm is independent of the mesh resolution when the latter contains at least 240 x 240 voxels, corresponding to a mesh spacing of 12.5 |im. As we show in figure 2c, melanophores tend to separate by between 60 and 100 |xm in our ABM results, so this mesh spacing is less than one quarter of the typical distance between agents.

At each mesh resolution in figure 4a, the estimated optimal value of "mm does not appear to depend greatly on the initial condition. For example, in the case of a mesh with Nb;n - N^st - 240, the maximum relative difference between the four parameter values is at most 6.5%. This similarity suggests that there is an inherent time scale at which migratory melanophore-melanophore interactions occur. Figure 4b, which presents the log L2 error for Nb;n - 240 as a function of "mm, further supports this conclusion. Although the errors associated with different initial conditions can vary by an order of magnitude, the minimum value of each (roughly convex) curve appears nearly identical and is located near the values shown in figure 4a.

Figure 5 presents snapshots of the EA ABM results across 104 realizations of equation (2.2) and the optimized PDE solution associated with the Box initial condition. The first row shows the expansion in time of the EA ABM support, i.e. the area occupied by the cells, due to melanophore-melanophore repulsion. For more intuition, we superimpose the cell positions from one ABM realization on our number-density results in this figure and throughout the paper. In all cases, we crop out approximately the upper half of cell positions. Visual inspection of cell positions in figure 5 suggests that melanophore-melanophore distances increase near the edge of the collective. Similarly, the speed at which the support expands appears to slow down for the EA ABM result, consistent with melanophores experiencing weaker forces from comparatively distant cells in this region.

We also observe in figure 5a-d that a band of high cell density emerges around the edge of the support which surrounds a ring-like region of low density. These bands may result from the combined effects of cell-cell repulsion and the fine mesh resolution that we use to sort agent positions in the EA solution. Repulsion causes cells at the edge of the collective to travel towards empty regions, while more centrally located agents move more slowly due to the balance of forces from their neighbours. When repulsion separates cells by distances greater than the mesh resolution, we expect regions of low density within the solution support to appear. These oscillatory bands should become less evident when the repulsive potentials in figure Id exhibit shallower gradients, as this permits cells to cluster more closely, or when coarser histograms with fewer bins are used to visualize the EA ABM data. As we discuss in electronic supplementary material, S6, the forces acting on xanthophores are about an order of magnitude smaller than those for melanophores, and we indeed observe less pronounced bands there. Notably, fitting to EA ABM data on coarser histograms leads to similar parameter estimates; see electronic supplementary material, S3, for details.

We present snapshots of the continuous model, equation (2.4), under our estimated value of %m in figure 5e-h. This PDE solution captures the dynamics of our example ABM realization significantly better than the case in figure 2b, when hmm - 1 o However, unlike the E A ABM result, the PDE does not exhibit bands of high and low cell density. This discrepancy can be further appreciated in figure 5i-k, which presents snapshots of the pointwise difference between the PDE and EA ABM solutions. Here, positive values indicate that the discrete solution is larger than the continuous one. The lack of bands in the PDE setting is likely because the mean-field assumption used to derive the continuous system is invalid where density is low. We do not expect this discrepancy to be as pronounced in models that include cell birth, as this mechanism increases density; see §3.3. Moreover, the PDE support expands more quickly than that of the ABM. This result is likely due to our choice of error function to fit "mm- Specifically, this parameter is biased towards values that produce accurate approximations in the bulk as these regions have a larger contribution to the L2 norm. Since we have already determined that the assumptions underlying the continuous model break down in low density regions, however, we choose to fit to the bulk of the cell density and focus on the L2 difference.

To demonstrate that our observations for the Box case are consistent across initial conditions, we compare the EA ABM and PDE dynamics for the Stripe and Two stripes initial conditions in figure 6; electronic supplementary material, figure S7, presents results for the Offset rectangles initial condition. In figure 6a-d, the column-averaged PDE solution, i.e. the solution average over the x variable, has a larger support than that of the EA ABM and does not exhibit oscillatory bands. (Comparing column averages is justified because both results are nearly uniform along the x-axis.) Nevertheless, the continuous solution closely approximates the EA ABM density, particularly in regions where the latter is high. For example, we find that the average pointwise error (over voxels where at least one of the EA ABM or PDE solution is non-zero) is equal to about 40 cells mm-2 at f = 0days, 24 cells mm-2 at t = 50days, 23 cells mm-2 at t = 100days and 23 cells mm-2 at t = 150 days (these correspond to roughly 10%, 6%, 6% and 6% of the maximum cell density of 400 cells mm-2, respectively). Both solutions invade empty space in time, and the speed of this travelling wavefront appears to slow as cells become more diffuse. For the Two stripes initial condition in figure 6e-h, the ABM and PDE predict that cells move into the initially empty space between stripes to approach a characteristic profile also observed in the one-stripe case. The EA ABM model does not appear to form oscillatory bands in the interstripe region, corroborating our hypothesis that these bands are more likely to arise near the edge of the solution support. In this case, the average pointwise error between the continuous and discrete data is roughly 18%, 11%, 11% and 11% of the maximum cell density at t = 0, 50, 100 and 150 days, respectively.

3.2. Cell birth

We identify the density-limiting parameter c+ and growth rate y in our IDE model, equation (2.7), by comparing with agent-based data from equation (2.5). Importantly, the dynamics of discrete-model proliferation, unlike cell migration, involve stochasticity beyond the initial condition. To gain intuition, we thus start with ID simulations: for each value of Nbire{l, 2,..., 10}, we compute the EA of 103 ABM realizations from an initial condition in which a single melanophore is placed at the origin in a ID domain. In figure 7a, we show the EA result for Nbir = l and the corresponding IDE model solution with the optimal values of y and c+ in figure 7b. The continuous solution appears to have a smaller radius of support than the EA ABM result at every time point; see figure 7c. This result holds across all Nbir values in figure 7A. While the IDE predicts a piecewise linear growth of the total number of cells, the corresponding EA ABM result increases linearly before slowly saturating as the domain fills, as we depict in figure 7e. This behaviour likely arises from our overcrowding condition that prevents cell densities from exceeding c+. As the domain fills with cells, it becomes less likely to select a location z that satisfies the overcrowding condition in the ABM. This reduces the population growth rate at later times. By contrast, the IDE model specifies that the support increases by the same amount at each time step until it reaches the domain boundaries. As we discuss in §4, capturing discrete model behaviour more accurately at higher cell numbers may require replacing y in our IDE with a density-dependent function.

Our ID simulations provide a baseline case to test our estimation process. As we note in §2.4, we employ a sequential procedure, first fitting c+ with y= I Q. I _1 and then estimating y with c+ fixed. In the ID case, this leads to optimal values c+ = 7.592 cells and 7=0.2822. If we instead estimate both parameters simultaneously, we find c+ = 7.430 cells and y= 0.2902. This is a difference of about 2.1% in c+ and 2.8% in y, suggesting that sequential estimation reduces computational complexity without strongly affecting parameter values. To understand if a coarser discrepancy measure based only on cell numbers at each time is sufficient, we also fit c+ and y by minimizing the squared difference in the total cell numbers over time; see electronic supplementary material, S2 and table S4, for the resulting parameter values. The corresponding parameter estimates differ from the density-based case by approximately 1.2% for and c+ and y, suggesting both error measures are reasonable. Both approaches also appear to exhibit similar sensitivity as parameters are varied (compare figure 7; electronic supplementary material, figure S2).

Figure 8 and electronic supplementary material, figure S8, respectively, show that proliferation in 2D broadens the solution support from the Box and Offset rectangles initial conditions over time, and the IDE model accurately captures the total cell mass of the ABM system for all Nbir values considered. Our estimated optimal values of c+ and y for these two initial conditions differ by about 2.5% and 0.31%, respectively, suggesting that our estimation procedure is robust to the initial condition. We also highlight that a region of higher density forms at the edge of the initial condition's support for both the ABM and IDE in figure 8a-h. Indeed, if z is near the support boundary, fi^oc covers only a fraction of the occupied domain, thereby meeting both conditions for birth. Conversely, the cell density at the centre of the domain is comparatively low throughout time because the total number of cells contained within discs of size |fiioc is already close to the threshold c+. Interestingly, as in the ID case with only proliferation, the ADM hA support is larger than that or the iDfc solution, the reverse or the behaviour that we observed for cell migration in figure 5.

3.3. Cell movement and proliferation

To obtain a full continuous model, we may substitute our estimated values of the migration scaling parameter %M; density-limiting parameter c+ and birth-rate scaling parameter y into equation (2.8). However, comparing this model to the dynamics of our full ABM shows that migration and proliferation have interwoven effects. To illustrate this phenomenon, we present a PDE solution with our optimal values of %M; c+ and y from §§3.1 and 3.2 at t = 70 days in figure 9a. We observe that this PDE model produces a significantly higher cell density than its discrete counterpart in figure 9b. This discrepancy occurs regardless of the value of Nbir, which influences the speed of cell birth. Related to this, we note that the long-time cell density in our ABM results is much lower when both mechanisms operate simultaneously than it is when only birth occurs; compare figure 8d and figure 9b. On the other hand, the inclusion of movement does not influence the long-time density of the continuous model solution; see the colourbar in figure 8h in comparison to the one in figure 9. Although we do not furnish these observations with an analytical explanation here, they demonstrate an interesting difference in how 'adding' mechanisms or terms impact PDE and ABM dynamics.

One approach to addressing these discrepancies is to refit all three scaling parameters ("mm, y and c+) simultaneously, and we present the results of this approach in electronic supplementary material, S3. (Indeed, we show there that the errors produced with a simultaneous estimation approach can be relatively small, although the model parameters may not all be identifiable.) Because we are interested in understanding the interplay of individual-based mechanisms of proliferation and movement in continuous models, however, we instead take a simpler theoretical approach. Namely, we note that the parameter c+ is largely responsible for controlling the maximum cell density over long time periods. (We determine this by integrating equation (2.8) over space and identifying the steady-state dynamics; this analysis reveals that equilibrium is reached when the density within any neighbourhood fi*oc is below c+.) In order to limit the maximum density to our estimated empirical value of 400 cells mm"2 [14,78], we let c+ = 4001 Q.\oc I " 7.0686 cells. As we show in figure 9c, using this value of c+, alongside our previously fitted values of %m and y, produces PDE densities that are much closer to the corresponding ABM results. We thus fix c+ = 7.0686 cells for the remainder of this paper, which allows us to highlight the time dynamics of our full PDE model in comparison to the EA ABM result with Box and Offset rectangles initial conditions in figures 10 and 11, respectively.

Figure 10j, which depicts the time evolution of the estimated radius of support for the PDE and EA ABM results, shows that a reasonably accurate continuous description of the combined model can be obtained by using the scaling parameters obtained from a modular approach. The supports of the ABM and PDE solutions both increase at roughly the same rate, although the degree to which the solutions agree can be affected by Nbir: when this parameter is small, the PDE solution travels at a faster rate than the ABM solution, whereas the opposite occurs when this value is large (greater than Nbir = 50 positions day-1). At intermediate values of Nbir (i.e. Nbir = 50 positions day-1), however, the ABM and PDE solution curves are almost identical. Figure lOi, which presents the number of cells over time, yields similar observations: the ABM and PDE solutions exhibit similar dynamics over the time period investigated here and there are certain values of Nbir for which the solution curves are nearly identical. Figure 11 further demonstrates that these observations do not depend on the specific choice of initial condition. Figures 10 and 11 demonstrate that combining movement with proliferation also dissipates the oscillatory bands that we observed for movement alone in figure 5. This is likely because the stochastic addition of cells in the birth model disrupts the regular cell spacing created by the movement model. Furthermore, the EA ABM and PDE solutions exhibit similar characteristic profiles without regions of high cell density around the edge of the initial condition support, in contrast to the birth-only model (figure 8).

4. Discussion

We presented a procedure for constructing experimentally interpretable continuous models of cell migration and birth in biologically relevant settings of low numbers of individuals and localized interactions which may lie outside the validity of the mean-field regime. Specifically we introduced and estimated scaling parameters in continuous models to account for realistic-i.e. relatively small and changing-numbers of cells with localized interactions. We applied this methodology to an illustrative, simplified example inspired by zebrafish pattern formation, in which we used a reduced ABM to generate individual-level data with biologically meaningful spatial and temporal units. Nonlocal rules for cell birth and migration, based on the ABM [14], informed our discrete and continuous descriptions and allowed us to transfer biological length scales and units to the macroscopic setting. Throughout our work, we stressed matching the spatio-temporal behaviour of our continuous and discrete models. We adopted a modular approach by estimating parameters in cases with either movement or birth before considering both mechanisms simultaneously. This allowed us to examine the specific contributions of each mechanism to self-organization and provided insight into their interplay in discrete and continuous settings.

We observed that the solutions of our continuous models expand at a different rate than EA ABM results and feature smoother profiles. Indeed, inaccuracies in mean-field descriptions for 'intermediate' numbers of individuals appear to be common in other biological phenomena described by simpler dynamics such as Fisher-KPP-type equations, cf. [85,86]. In fact, both references analytically derive corrections to the wave speed, a procedure we cannot adopt due to our use of off-lattice models. However, this substantiates our introduction of scaling parameters to handle the discrepancy. By introducing and estimating parameters that rescale the time variable, we produced more accurate descriptions of agent-based movement or birth. However, when we used the same parameter values in a continuous model of both cell migration and birth in §3.3, the PDE did not produce close estimates of the full ABM. Specifically, our full continuous model yielded larger long-time densities than the EA ABM results, motivating us to re-estimate the threshold value c+ with a theoretical approach. This generated a more faithful continuous description and highlighted that the effects of movement and proliferation are not simply additive. We thus stress that parameters must be fitted to data in which all mechanisms of interest act simultaneously, in order to capture their interplay. This is particularly crucial for contexts such as cancer biology, where cell migration, proliferation and death are known to play critical roles in tumour progression and immune response [87].

Our results highlight how choices in numerical implementation affect parameter estimates and suggest several directions for future work that may improve our approach. For example, the optimal value of our parameter controlling the time scale of cell migration ("mm) appears to be independent of the initial condition and the mesh resolution that we used to construct PDE solutions, provided the latter is sufficiently refined. One drawback of our current approach, however, is that we may need to estimate hmm and y for each new choice of discrete rules governing migration and cell birth, respectively, because these rules perturb the short-range interactions between relatively small numbers of cells. This naturally leads to the question of whether an analytic expression can be derived for these parameters. Several coarse-graining techniques that take into account higher-order correlations between cells in on- or off-lattice models may produce scaling factors similar to those introduced in this paper, but these may only apply to certain classes of ABMs [88]. Alternatively, it may be possible to estimate scaling parameters of continuous models by investigating the convergence of EA ABM results to features of their solutions such as the speed of solution propagation, as in [85,86] for on-lattice models; by accounting for the dynamics of the two-particle distribution, as in [73,74]; or by exploring scaling relationships as in [89]. Adapting these approaches to our setting is an interesting avenue for future work.

Additionally, our continuous models more accurately represent ABM results within the bulk of the solution support because the L2 norm more strongly penalizes discrepancies there. In the future, other norms, such as the L00 error, could be used to match the solution supports given by our discrete and continuous models. Replacing our birth-rate scaling parameter y with a density-dependent function- through either rigorous derivation or an equation-learning approach [90]-is another exciting future direction. In particular, because cell proliferation in the ABM involves selecting positions uniformly at random from the domain each day, the chance that we select a location z that permits birth appears to depend in a nonlinear way on the solution support. More generally, our computational study does not provide theoretical explanations for our parameter values, and we plan to build on the intuition that we established here to develop these arguments in the future.

To simplify our initial study, we considered the dynamics of one cell population (i.e. melanophores in the main text and xanthophores in electronic supplementary material, S6), but pattern formation in zebrafish skin involves multiple cell types and longer-range interactions, as we highlight in electronic supplementary material, figure SI. Future work may extend our pipeline to construct more realistic continuous models with multiple cell types and interaction neighbourhoods. Related to this, the initial conditions that we designed allowed us to make one-to-one comparisons between discrete- and continuous-model densities, but this may not always be possible. More realistic zebrafish models (i.e. [14,37,39,40]) produce patterns that are more complicated than our box and stripe motifs. This means that ensemble-averaging stochastic ABM realizations may not retain information about the length scales inherent in patterns. For such cases, fitting parameters based on summary statistics (e.g. pair-correlation functions [91], pattern-simplicity scores [92] or persistent-homology approaches [93]) may be more useful, and we plan to address this in future work. These and other directions move us toward constructing interpretable, analytically tractable continuous models of self-organization, increasing our understanding of biological pattern formation more broadly. Ethics. This work did not require ethical approval from a human subject or animal welfare committee. Data accessibility. Data and relevant code for this research work are stored in GitHub: https://github.com/ wdmartmson/Self-Organization-One-Species/tree/firstversion and have been archived within the Zenodo repository: https://doi.org/10.5281/zenodo.11160930 [94]. Electronic supplementary material, figures and videos are available on Figshare at https://figshare.com/projects/Linking_discrete_and_continuous_models_of_cell_birth_ and_migration_in_one_population/171234. Raw data used to generate figures in the main text and electronic supplementary material have been archived in the Dryad repository https://doi.org/10.5061/dryad.s4mw6m9cb [95]. Supplementary material is available online [96]. Declaration of Al use. We have not used Al-assisted technologies in creating this article. Authors' contributions. W.D.M.: conceptualization, data curation, formal analysis, investigation, methodology, software, validation, visualization, writing-original draft, writing-review and editing; A.V.: conceptualization, data curation, investigation, methodology, software, validation, visualization, writing-original draft, writing-review and editing; VI.M.: conceptualization, data curation, investigation, methodology, software, validation, writing-original draft, writing-review and editing; C.V.: conceptualization, data curation, investigation, methodology, software, validation, visualization, writing-original draft, writing-review and editing; J.A.C.: conceptualization, funding acquisition, methodology, project administration, supervision, writing-review and editing. All authors gave final approval for publication and agreed to be held accountable for the work performed therein. Conflict of interest declaration. We declare we have no competing interests. Funding. J.A.C. and W.D.M. were supported by the Advanced Grant Non-local-CPD (Non-local PDEs for Complex Particle Dynamics: Phase Transitions, Patterns and Synchronization) of the European Research Council Executive Agency (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 883363). J.A.C. is also partially supported by the 'Maria de Maeztu' Excellence Unit IMAG, reference CEX2020-001105-M, funded by MCIN/AEI/10.13039/501100011033/. W.D.M., J.A.C. and A.V. would also like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme Vlathematics of Movement, where work on this paper was undertaken. This work was supported by EPSRC grant no. EP/R014604/1. C.V. acknowledges support from the Dr Perry James (Jim) Browne Research Centre on Vlathematics and its Applications (University of Sussex). Acknowledgements. We are grateful to Shigeru Kondo for helpful discussions during an early stage of this research.

References

References

1. Mogilner A, Edelstein-Keshet L. 1999 A non-local model for a swarm. J. Math. Biol. 38, 534-570. (doi:10.1007/s002850050158)

2. D'Orsogna MR, Chuang YL, Bertozzi AL, Chayes LS. 2006 Self-propelled particles with soft-core interactions: patterns, stability, and collapse. Phys. Rev. Lett. 96,104302. (doi:10.1103/PhysRevLett.96.104302)

3. Cucker F, Smale S. 2007 Emergent behavior in flocks. SEE Trans. Mom. Control 52, 852-862. (doi:10.1109/TAC.2007.895842)

4. Carrillo JA, Fornasier M, Toscani G, Vecil F. 2010 Particle, kinetic, and hydrodynamic models of swarming. In Mathematical modeling of collective behavior in socio-economic and life sciences, pp. 297-336. New York, NY: Springer.

5. Amack JD, Manning ML. 2012 Knowing the boundaries: extending the differential adhesion hypothesis in embryonic cell sorting. Science 338, 212-215. (doi:10.1126/ science.1223953)

6. Burger M, Francesco MD, Fagioli S, Stevens A. 2018 Sorting phenomena in a mathematical model for two mutually attracting/repelling species. 51AM J. Math. Anal. 50, 3210-3250. (doi:10.1137/17M1125716)

7. Carrillo JA, Murakawa H, Sato M, Togashi H, Trush 0. 2019 A population dynamics model of cell-cell adhesion incorporating population pressure and density saturation. / Theor. Biol. 474, 14-24. (doi:10.1016/j.jtbi.2019.04.023)

8. Buttenschon A, Edelstein-Keshet L. 2020 Bridging from single to collective cell migration: a review of models and links to experiments. PLo5 Comput. Biol. 16, e1008411. (doi:10.1371/journal.pcbi.1008411)

9. Tsai TYC, Garner RM, Megason SG. 2022 Adhesion-based self-organization in tissue patterning. Anna. Rev. Cell. Dev. Biol. 38,349-374. (doi:10.1146/annurev-cellbio-120420-100215)

10. Frohnhofer HG, Krauss J, Maischein HM, Niisslein-Volhard C. 2013 Iridophores and their interactions with other chromatophores are required for stripe formation in zebrafish. Development 140, 2997-3007.

11. Irion U, Niisslein-Volhard C. 2019 the identification of genes involved in the evolution of color patterns in fish. Curr. Opin. Genet. Dev. 57,31-38. (doi:10.1016/j.gde.2019. 07.002)

12. Parichy DM. 2021 Evolution of pigment cells and patterns: recent insights from teleost fishes. Curr. Opin. Genet. Dev. 69, 88-96. (doi:10.1016/j.gde.2021.02.006)

13. Kondo S, Watanabe M, Miyazawa S. 2021 Studies of luring pattern formation in zebrafish skin. Phil. Trans. R. Sot A 379, 20200274. (doi:10.1098/rsta.2020.0274)

14. Volkening A, Sandstede B. 2015 Modelling stripe formation in zebrafish: an agent-based approach. J. ft Soc. Interface 12, 20150812. (doi:10.1098/rsif.2015.0812)

15. Fadeev A, Krauss J, Frohnhofer HG, Irion U, Niisslein-Volhard C. 2015 Tight Junction Protein la regulates pigment cell organisation during zebrafish colour patterning. eLife 4, e06545. (doi:10.7554/eLife.06545)

16. Nakamasu A, Takahashi G, Kanbe A, KondoS. 2009 Interactions between zebrafish pigment cells responsible for the generation of Turing patterns. Proc. Natl Acad. Sci. USA 106, 8429-8434. (doi:10.1073/pnas.0808622106)

17. Singh AP, Schach U, Niisslein-Volhard C. 2014 Proliferation, dispersal and patterned aggregation of iridophores in the skin prefigure striped colouration of zebrafish. Nat. Cell Biol. 16, 604-611. (doi:10.1038/ncb2955)

18. Gur D et al. 2020 In situ differentiation of iridophore crystallotypes underlies zebrafish stripe patterning. Nat. Commun. 11, 6391. (doi:10.1038/s41467-020-20088-1)

19. Hamada H, Watanabe M, Lau HE, Nishida T, Hasegawa T, Parichy DM, Kondo S. 2014 Involvement of Delta/Notch signaling in zebrafish adult pigment stripe patterning. Development 141,318-324. (doi:10.1242/dev.099804)

20. Inaba M, Yamanaka H, Kondo S. 2012 Pigment pattern formation by contact-dependent depolarization. Science 335, 677-677. (doi:10.1126/science.1212821)

21. Yamaguchi M, Yoshimoto E, Kondo S. 2007 Pattern regulation in the stripe of zebrafish suggests an underlying dynamic and autonomous mechanism. Proc. Natl Acad. Sci. USA 104, 4790-4793. (doi:10.1073/pnas.0607790104)

22. Kroll F ef al. 2021 A simple and effective F0 knockout method for rapid screening of behaviour and other complex phenotypes. eLife 10, e59683. (doi:10.7554/eLife.59683)

23. Byrne HM. 2010 Dissecting cancer through mathematics: from the cell to the animal model. Nat. Rev. Cancer 10, 221-230. (doi:10.1038/nrc2808)

24. Gompper G ef al. 2020 The 2020 motile active matter roadmap. J. Phys. Condens. Matter 32, 193001. (doi:10.1088/1361-648X/ab6348)

25. Stillman NR, Mayor R. 2023 Generative models of morphogenesis in developmental biology. Semin. Cell Dev. Biol. 147, 83-90. (doi:10.1016/j.semcdb.2023.02.001)

26. Volkening A. 2020 Linking genotype, cell behavior, and phenotype: multidisciplinary perspectives with a basis in zebrafish patterns. Curr. Opin. Genet. Dev. 63, 78-85. (doi:10.1016/j.gde.2020.05.010)

27. Alert R, Trepat X. 2020 Physical models of collective cell migration. Annu. Rev. Condens. Matter Phys. 11, 77-101. (doi:10.1146/annurev-conmatphys-031218-013516)

28. Metzcar J, Wang Y, Heiland R, Macklin P. 2019 A review of cell-based computational modeling in cancer biology. ICO Clin. Cancer Inform. 2,1-13. (doi:10.1200/CCI.18. 00069)

29. Deutsch A. 2015 Cellular automaton models for collective cell behaviour. In Cellular automata and discrete complex systems (ed. J Kari), pp. 1-10. Berlin, Germany: Springer. (doi:10.1007/978-3-662-47221-7_1)

30. Deutsch A, Dormann S. 2017 Cellular automaton modeling of biological pattern formation, characterization, examples, and analysis. New York, NY: Springer.

31. Hirashima T, Rens EG, Merks RMH. 2017 Cellular Potts modeling of complex multicellular behaviors in tissue morphogenesis. Dev. Growth Differ. 59,329-339. (doi:10. HH/dgd.12358)

32. Rens EG, Edelstein-Keshet L. 2019 From energy to cellular forces in the cellular Potts model: an algorithmic approach. PLoS Comput. Biol. 15, e1007459. (doi:10.1371/ journal.pcbi.1007459)

33. Alt S, Ganguly P, Salbreux G. 2017 Vertex models: from cell mechanics to tissue morphogenesis. Phil. Trans, ft Soc. B 372, 20150520. (doi:10.1098/rstb.2015.0520)

34. Fletcher A, Osterfield M, Baker R, Shvartsman S. 2014 Vertex models of epithelial morphogenesis. Biophys. J. 106, 2291-2304. (doi:10.1016/j.bpj.2013.11.4498)

35. Bullara D, De Decker Y 2015 Pigment cell movement is not required for generation of Turing patterns in zebrafish skin. Wof. Commun. 6,6971. (doi:10.1038/ncomms7971)

36. Moreira J, Deutsch A. 2005 Pigment pattern formation in zebrafish during late larval stages: a model based on local interactions. Dev. Dyn. 232,33-42. (doi:10.1002/dvdy. 20199)

37. Owen JP, Kelsh RN, Yates CA. 2020 A quantitative modelling approach to zebrafish pigment pattern formation. eLife 9, e52998. (doi:10.7554/eLife.52998)

38. Caicedo-Carvajal CE, Shinbrot T. 2008 In silica zebrafish pattern formation. Dev. Biol. 315,397-403. (doi:10.1016/j.ydbio.2007.12.036)

39. Volkening A, Abbott MR, Chandra N, Dubois B, Lim F, Sexton D, Sandstede B. 2020 Modeling stripe formation on growing zebrafish tailfins. Bull. Math. Biol. 82,56. (doi:10.1007/s11538-020-00731-0)

40. Volkening A, Sandstede B. 2018 Iridophores as a source of robustness in zebrafish stripes and variability in Danio patterns. Nat. Commun. 9,3231. (doi:10.1038/s41467-018-05629-z)

41. Osborne JM, Fletcher AG, Pitt-Francis JM, Maini PK, Gavaghan DJ. 2017 Comparing individual-based approaches to modelling the self-organization of multicellular tissues. PLoS Comput. Biol. 13, e1005387. (doi:10.1371/journal.pcbi.100S387)

42. Murray JD. 2002 Mathematical biology: 1. An introduction. New York, NY: Springer.

43. Perthame B. 2006 Transport eguations in biology. Berlin, Germany: Springer Science & Business Media.

44. Turing AM. 1952 The chemical basis of morphogenesis. Phil. Trans, ft Soc. Land. B 237, 37-72. (doi:10.1098/rstb.l952.0012)

45. Maini P, Painter K, Chau HP. 1997 Spatial pattern formation in chemical and biological systems. J. Chem. Soc. 93, 3601-3610. (doi:10.1039/a702602a)

46. Marciniak-Czochra A, Karch G, Suzuki K. 2017 Instability of Turing patterns in reaction-diffusion-ODE systems.J. Math. Biol. 74,583-618. (doi:10.1007/s00285-016-1035-z)

47. Armstrong NJ, Painter KJ, Sherratt JA. 2006 A continuum approach to modelling cell-cell adhesion. J. Theor. Biol. 243, 98-113. (doi:10.1016/j.jtbi.2006.05.030)

48. Murakawa H, Togashi H. 2015 Continuous models for cell-cell adhesion. J. Theor. Biol. 374, 1-12. (doi:10.1016/j.jtbi.2015.03.002)

49. Carrillo JA, Huang Y, Schmidtchen M. 2018 Zoology of a nonlocal cross-diffusion model for two species. SIAM J. Appl. Math. 78,1078-1104. (doi:10.1137/17M1128782)

50. Carrillo JA, Colombi A, Scianna M. 2018 Adhesion and volume constraints via nonlocal interactions determine cell organisation and migration profiles. J. Theor. Biol. 445, 75-91. (doi:10.1016/j.jtbi.2018.02.022)

51. Gaffney EA, Seirin Lee S. 2015 The sensitivity of Turing self-organization to biological feedback delays: 2D models of fish pigmentation. Math. Med. Biol. 32,57-79. (doi:10.1093/imammb/dqt017)

52. Konow C, Li Z, Shepherd S, Bullara D, Epstein IR. 2021 Influence of survival, promotion, and growth on pattern formation in zebrafish skin. Sci. Rep. 11, 9864. (doi:10. 1038A41598-021-89116-4)

53. Bloomfield JM, Painter KJ, Sherratt JA. 2011 How does cellular contact affect differentiation mediated pattern formation? Bull. Math. Biol. 73,1529-1558. (doi:10.1007/ si 1538-010-9578-4)

54. Kondo S. 2017 An updated kernel-based Turing model for studying the mechanisms of biological pattern formation. J. Theor. Biol. 414,120-127. (doi:10.1016/j.jtbi.2016. 11.003)

55. Painter KJ, Bloomfield JM, Sherratt JA, Gerisch A. 2015 A nonlocal model for contact attraction and repulsion in heterogeneous cell populations. Bull. Math. Biol. 77, 1132-1165. (doi:10.1007/s11538-015-0080-x)

56. Giacomin G, Lebowitz JL. 1997 Phase segregation dynamics in particle systems with long range interactions 1. Macroscopic limits. J. Stat. Phys. 87,37-61. (doi:10.1007/ BF02181479)

57. Painter KJ, Hillen T. 2002 Volume-filing and quorum-sensing in models for chemosensitive movement. Can. Appl. Math. Quart. 10,501-543.

58. Penington 0, Hughes BD, Landman KA. 2011 Building macroscale models from microscale probabilistic models: a general probabilistic approach for nonlinear diffusion and multispecies phenomena. Phys. Rev. E 84, 041120. (doi:10.1103/PhysRevE.84.041120)

59. Hillen T, Painter KJ. 2013 Transport and anisotropic diffusion models for movement in oriented habitats. In Dispersal, individual movement and spatial ecology (eds M Lewis, P Maini, S Petrovskii), pp. 177-222. Berlin, Germany: Springer. (doi:10.1007/978-3-642-35497-7_7)

60. Carrillo JA, Choi YP, Hauray M. 2014 The derivation of swarming models: mean-feld limit and Wasserstein distances. In Collective dynamic from bacteria to crowds: an excursion through modeling, analysis and simulation (eds A Muntean, F Toschi), pp. 1-46. Vienna, Austria: Springer. (doi:10.1007/978-3-7091 -1785-9_1)

61. Chen Y, Kolokolnikov T. 2014 A minimal model of predator-swarm interactions. J. R Soc. Interface 11, 20131208. (doi:10.1098/rsif.2013.1208)

62. Bruna M, Burger M, Pietschmann JF, Wolfram MT. 2022 Active crowds. In Active particles, vol. 3 (eds N Bellomo, JA Carrillo, E Tadmor), pp. 35-73. Cham, Switzerland: Birkhauser. (doi:10.1007/978-3-030-93302-9_2)

63. Champagnat N, Ferriere R, Meleard S. 2006 Unifying evolutionary dynamics: from individual stochastic processes to macroscopic models. Theor. Popul. Biol. 69, 297-321. (doi :10.1016/j.tpb.2005.10.004)

64. Champagnat N, Ferriere R, Meleard S. 2008 Individual-based probabilistic models of adaptive evolution and various scaling approximations. Prog. Probab. 59, 75. (doi:10. 1007/978-3-7643-8458-6_6)

65. Champagnat N, Ferriere R, Meleard S. 2008 From individual stochastic processes to macroscopic models in adaptive evolution. Stock Models 24, 2-44. (doi:10.1080/ 15326340802437710)

66. Hansen JP, McDonald IR. 2006 Theory of simple liguids. London, UK: Academic Press.

67. Lushnikov PM, Chen N, Alber M. 2008 Macroscopic dynamics of biological cells interacting via chemotaxis and direct contact. Phys. Rev. E 78, 061904. (doi:10.1103/ PhysRevE.78.061904)

68. Simpson MJ, Landman KA, Hughes BD. 2009 Multi-species simple exclusion processes. Phys. A: Stat. Mech. 388,399-406. (doi:10.1016/j.physa.2008.10.038)

69. Simpson MJ, Baker RE. 2011 Corrected mean-feld models for spatially dependent advection-diffusion-reaction phenomena. Phys. Rev. £83,51922. (doi:10.1103/PhysRevE. 83.051922)

70. Markham DC, Simpson MJ, Baker RE. 2013 Simplified method for including spatial correlations in mean-feld approximations. Phys. Rev. E 87, 62702. (doi:10.1103/ PhysRevE.87.062702)

71. Wieczorek R. 2023 Hydrodynamic limit of a stochastic model of proliferating cells with chemotaxis. Kinet. Relat. Models 16,373-393. (doi:10.3934/krm.2022032)

72. Binny RN, Plank MJ, James A. 2015 Spatial moment dynamics for collective cell movement incorporating a neighbour-dependent directional bias. J. R. Soc. Interface 12, 20150228. (doi:10.1098/rsif.2015.0228)

73. Binny RN, James A, Plank MJ. 2016 Collective cell behaviour with neighbour-dependent proliferation, death and directional bias. Bull. Math. Biol. 78, 2277-2301. (doi:10. 1007/sl 1538-016-0222-9)

74. Middleton AM, Fleck C, Grima R. 2014 A continuum approximation to an off-lattice individual-cell based model of cell migration and adhesion. J. Theor. Biol. 359, 220-232. (doi:10.1016/j.jtbi.2014.06.011)

75. Bruna M, Chapman SJ, Robinson M. 2017 Diffusion of particles with short-range interactions. SIAM J. Appl. Math. 77, 2294-2316. (doi:10.1137/17M1118543)

76. Johnston ST, Baker RE, Simpson MJ. 2017 A new and accurate continuum description of moving fronts. New J. Phys. 19, 033010. (doi:10.1088/1367-2630/aa5bf5)

77. Bruna M, Chapman SJ. 2012 Excluded-volume effects in the diffusion of hard spheres. Phys. Rev. E 85, 011103. (doi:10.1103/PhysRevE.85.011103)

78. Takahashi G, Kondo S. 2008 Melanophores in the stripes of adult zebrafish do not have the nature to gather, but disperse when they have the space to move. Pigment Cell Melanoma Res. 21, 677-686. (doi:10.1111/j.l755-148X.2008.00504.X)

79. Parichy DM, Turner JM. 2003 Zebrafish puma mutant decouples pigment pattern and somatic metamorphosis. Dev. Biol. 256, 242-257. (doi:10.1016/S0012-1606(03)00015-0)

80. Patterson LB, Bain EJ, Parichy DM. 2014 Pigment cell interactions and differential xanthophore recruitment underlying zebrafish stripe reiteration and Danio pattern evolution. Nat. Commun. 5,5299. (doi:10.1038/ncomms6299)

81. Walderich B, Singh AP, Mahalwar P, Niisslein-Volhard C. 2016 Homotypic cell competition regulates proliferation and tiling of zebrafish pigment cells during colour pattern formation. Nat. Commun. 7, 11462. (doi:10.1038/ncomms11462)

82. Golse F. 2003 The mean-feld limit for the dynamics of large particle systems. Journees Equations aux derivees Partielles.

83. Di Francesco M, Fagioli S. 2013 Measure solutions for non-local interaction PDEs with two species. Nonlinearity 26, 2777. (doi:10.1088/0951 -7715/26/10/2777)

84. Chiari G, Delitala ME, Morselli D, Scianna M. 2022 A hybrid modeling environment to describe aggregates of cells heterogeneous for genotype and behavior with possible phenotypic transitions. Int. J. Non-Linear Mech. 144, 104063. (doi:10.1016/j.ijnonlinmec.2022.104063)

85. Brunet E, Derrida B. 1997 Shift in the velocity of a front due to a cutoff. Phys. Rev. E 56, 2597-2604. (doi:10.1103/PhysRevE.56.2597)

86. Sontag A, Rogers T, Yates CA. 2023 Stochastic drift in discrete waves of nonlocally interacting particles. Phys. Rev. E 107, 014128. (doi:10.1103/PhysRevE.107.014128)

87. Bull JA, Byrne HM. 2022 The hallmarks of mathematical oncology. Proc. IEEE 110, 523-540. (doi:10.1109/JPROC2021.3136715)

88. Bruna M, Chapman SJ, Robinson M. 2017 Diffusion of particles with short-range interactions. SIAM J. Appl. Math. 77, 2294-2316. (doi:10.1137/17M1118543)

89. Franke F, Aland S, Bohme HJ, Voss-Bohme A, Lange S. 2022 Is cell segregation like oil and water: asymptotic versus transitory regime. PLoS Comput. Biol. 18, e1010460. (doi:10.1371/journal.pcbi.1010460)

90. Nardini JT, Baker RE, Simpson MJ, Flores KB. 2021 Learning differential equation models from stochastic agent-based model simulations. J. R. Soc Interface 18, 20200987. (doi:10.1098/rsif.2020.0987)

91. Bull JA, Byrne HM. 2023 Quantification of spatial and phenotypic heterogeneity in an agent-based model of tumour-macrophage interactions. PLoS Comput. Biol. 19, e1010994. (doi:10.1371/journal.pcbi.1010994)

92. Miyazawa S. 2020 Pattern blending enriches the diversity of animal colorations. Sci. Adv. 6, eabb9107. (doi:10.1126/sciadv.abb9107)

93. McGuirl MR, Volkening A, Sandstede B. 2020 Topological data analysis of zebrafsh patterns. Proc. Natl Acad. Sci. USA 117,5113-5124. (doi:10.1073/pnas.1917763117)

94. Martinson WD. 2024 wdmartinson/Self-Organization-One-Species: Self-Organization-One-Species-Code-for-Article. zenodo. (doi:10.5281/zenodo.11160930)

95. William M, Alexandria V, Markus S, Chandrasekhar V, Jose C. 2024 Data from: Linking continuous and discrete models of cell birth and migration. Dryad Digital Repository. (doi:10.5061/dryad.s4mw6m9cb)

96. Martinson WD, Volkening A, Schmidtchen M, Venkataraman C, Carrillo JA. 2024 Linking discrete and continuous models of cell birth and migration. Figshare. (doi:10.6084/ m9.fgshare.c.7303179)