Exclusion processes on a growing domain

A Model for Cell Proliferation in a Developing Organism

In mathematical biology, there is a great deal of interest in producing continuum models by scaling discrete agent-based models governed by local stochastic rules. We discuss a particular example of this approach: a model for the proliferation of neural crest cells that can help us understand the development of Hirschprung's disease, a potentially-fatal condition in which the enteric nervous system of a new-born child does not extend all the way through the intestine and colon. Our starting point is a discrete-state, continuous-time Markov chain model proposed by Hywood et al. (2013a) for the location of the neural crest cells that make up the enteric nervous system. Hywood et al. (2013a) scaled their model to derive an approximate second order partial differential equation describing how the limiting expected number of neural crest cells evolve in space and time. In contrast, we exploit the relationship between the above-mentioned Markov chain model and the well-known Yule-Furry process to derive the exact form of the scaled version of the process. Furthermore, we provide expressions for other features of the domain agent occupancy process, such as the variance of the marginal occupancy at a particular site, the distribution of the number of agents that are yet to reach a given site and a stochastic description of the process itself.

Modelling uniaxial non-uniform yeast colony growth: Comparing an agent-based model and continuum approximations

Biological experiments have shown that yeast can be restricted to grow in a uniaxial direction, vertically upwards from an agar plate to form a colony. The growth occurs as a consequence of cell proliferation driven by a nutrient supply at the base of the colony, and the height of the colony has been observed to increase linearly with time. Within the colony the nutrient concentration is non-constant and yeast cells throughout the colony will therefore not have equal access to nutrient, resulting in non-uniform growth. In this work, an agent based model is developed to predict the microscopic spatial distribution of labelled cells within the colony when the probability of cell proliferation can vary in space and time. We also describe a method for determining the average trajectories or pathlines of labelled cells within a colony growing in a uniaxial direction, enabling us to connect the microscopic and macroscopic behaviours of the system. We present results for six cases, which involve different assumptions for the presence or absence of a quiescent region (where no cell proliferation occurs), the size of the proliferative region, and the spatial variation of proliferation rates within the proliferative region. These six cases are designed to provide qualitative insight into likely growth scenarios whilst remaining amenable to analysis. We compare our macroscopic results to experimental observations of uniaxial colony growth for two cases where only a fixed number of cells at the base of the colony can proliferate. The model predicts that the height of the colony will increase linearly with time in both these cases, which is consistent with experimental observations. However, our model shows how different functional forms for the spatial dependence of the proliferation rate can be distinguished by tracking the pathlines of cells at different positions in the colony. More generally, our methodology can be applied to other biological systems exhibiting uniaxial growth, providing a framework for classifying or determining regions of uniform and non-uniform growth.

Detection and characterization of chemotaxis without cell tracking

Swarming has been observed in various biological systems from collective animal movements to immune cells. In the cellular context, swarming is driven by the secretion of chemotactic factors. Despite the critical role of chemotactic swarming, few methods to robustly identify and quantify this phenomenon exist. Here, we present a novel method for the analysis of time series of positional data generated from realizations of agent-based processes. We convert the positional data for each individual time point to a function measuring agent aggregation around a given area of interest, hence generating a functional time series. The functional time series, and a more easily visualized swarming metric of agent aggregation derived from these functions, provide useful information regarding the evolution of the underlying process over time. We extend our method to build upon the modelling of collective motility using drift-diffusion partial differential equations (PDEs). Using a functional linear model, we are able to use the functional time series to estimate the drift and diffusivity terms associated with the underlying PDE. By producing an accurate estimate for the drift coefficient, we can infer the strength and range of attraction or repulsion exerted on agents, as in chemotaxis. Our approach relies solely on using agent positional data. The spatial distribution of diffusing chemokines is not required, nor do individual agents need to be tracked over time. We demonstrate our approach using random walk simulations of chemotaxis and experiments investigating cytotoxic T cells interacting with tumouroids.

Gaussian process approximations for multicolor Pólya urn models

Motivated by mathematical tissue growth modelling, we consider the problem of approximating the dynamics of multicolor Pólya urn processes that start with large numbers of balls of different colors and run for a long time. Using strong approximation theorems for empirical and quantile processes, we establish Gaussian process approximations for the Pólya urn processes. The approximating processes are sums of a multivariate Brownian motion process and an independent linear drift with a random Gaussian coefficient. The dominating term between the two depends on the ratio of the number of time steps n to the initial number of balls N in the urn. We also establish an upper bound of the form $c(n^{-1/2}+N^{-1/2})$ for the maximum deviation over the class of convex Borel sets of the step-n urn composition distribution from the approximating normal law.

From a discrete model of chemotaxis with volume-filling to a generalized Patlak-Keller-Segel model

We present a discrete model of chemotaxis whereby cells responding to a chemoattractant are seen as individual agents whose movement is described through a set of rules that result in a biased random walk. In order to take into account possible alterations in cellular motility observed at high cell densities (i.e. volume-filling), we let the probabilities of cell movement be modulated by a decaying function of the cell density. We formally show that a general form of the celebrated Patlak-Keller-Segel (PKS) model of chemotaxis can be formally derived as the appropriate continuum limit of this discrete model. The family of steady-state solutions of such a generalized PKS model are characterized and the conditions for the emergence of spatial patterns are studied via linear stability analysis. Moreover, we carry out a systematic quantitative comparison between numerical simulations of the discrete model and numerical solutions of the corresponding PKS model, both in one and in two spatial dimensions. The results obtained indicate that there is excellent quantitative agreement between the spatial patterns produced by the two models. Finally, we numerically show that the outcomes of the two models faithfully replicate those of the classical PKS model in a suitable asymptotic regime.

Agent-based modeling of morphogenetic systems: Advantages and challenges

The complexity of morphogenesis poses a fundamental challenge to understanding the mechanisms governing the formation of biological patterns and structures. Over the past century, numerous processes have been identified as critically contributing to morphogenetic events, but the interplay between the various components and aspects of pattern formation have been much harder to grasp. The combination of traditional biology with mathematical and computational methods has had a profound effect on our current understanding of morphogenesis and led to significant insights and advancements in the field. In particular, the theoretical concepts of reaction–diffusion systems and positional information, proposed by Alan Turing and Lewis Wolpert, respectively, dramatically influenced our general view of morphogenesis, although typically in isolation from one another. In recent years, agent-based modeling has been emerging as a consolidation and implementation of the two theories within a single framework. Agent-based models (ABMs) are unique in their ability to integrate combinations of heterogeneous processes and investigate their respective dynamics, especially in the context of spatial phenomena. In this review, we highlight the benefits and technical challenges associated with ABMs as tools for examining morphogenetic events. These models display unparalleled flexibility for studying various morphogenetic phenomena at multiple levels and have the important advantage of informing future experimental work, including the targeted engineering of tissues and organs.

Understanding interactions between populations: Individual based modelling and quantification using pair correlation functions

Understanding the underlying mechanisms that produce the huge variety of swarming and aggregation patterns in animals and cells is fundamental in ecology, developmental biology, and regenerative medicine, to name but a few examples. Depending upon the nature of the interactions between individuals (cells or animals), a variety of different large-scale spatial patterns can be observed in their distribution; examples include cell aggregates, stripes of different coloured skin cells, etc. For the case where all individuals are of the same type (i.e., all interactions are alike), a considerable literature already exists on how the collective organisation depends on the inter-individual interactions. Here, we focus on the less studied case where there are two different types of individuals present. Whilst a number of continuum models of this scenario exist, it can be difficult to compare these models to experimental data, since real cells and animals are discrete. In order to overcome this problem, we develop an agent-based model to simulate some archetypal mechanisms involving attraction and repulsion. However, with this approach (as with experiments), each realisation of the model is different, due to stochastic effects. In order to make useful comparisons between simulations and experimental data, we need to identify the robust features of the spatial distributions of the two species which persist over many realisations of the model (for example, the size of aggregates, degree of segregation or intermixing of the two species). In some cases, it is possible to do this by simple visual inspection. In others, the features of the pattern are not so clear to the unaided eye. In this paper, we introduce a pair correlation function (PCF), which allows us to analyse multi-species spatial distributions quantitatively. We show how the differing strengths of inter-individual attraction and repulsion between species give rise to different spatial patterns, and how the PCF can be used to quantify these differences, even when it might be impossible to recognise them visually.

How domain growth is implemented determines the long-term behavior of a cell population through its effect on spatial correlations

Domain growth plays an important role in many biological systems, and so the inclusion of domain growth in models of these biological systems is important to understanding how these systems function. In this work we present methods to include the effects of domain growth on the evolution of spatial correlations in a continuum approximation of a lattice-based model of cell motility and proliferation. We show that, depending on the way in which domain growth is implemented, different steady-state densities are predicted for an agent population. Furthermore, we demonstrate that the way in which domain growth is implemented can result in the evolution of the agent density depending on the size of the domain. Continuum approximations that ignore spatial correlations cannot capture these behaviors, while those that account for spatial correlations do. These results will be of interest to researchers in developmental biology, as they suggest that the nature of domain growth can determine the characteristics of cell populations.

A mechanistic study on tumour spheroid formation in thermosensitive hydrogels: experiments and mathematical modelling

Thermo-reversible microgels to culture and harvest uniform-sized tumour spheroids with a narrow size-distribution.

Survival probability for a diffusive process on a growing domain

We consider the motion of a diffusive population on a growing domain, 0<x<L(t), which is motivated by various applications in developmental biology. Individuals in the diffusing population, which could represent molecules or cells in a developmental scenario, undergo two different kinds of motion: (i) undirected movement, characterized by a diffusion coefficient, D, and (ii) directed movement, associated with the underlying domain growth. For a general class of problems with a reflecting boundary at x=0, and an absorbing boundary at x=L(t), we provide an exact solution to the partial differential equation describing the evolution of the population density function, C(x,t). Using this solution, we derive an exact expression for the survival probability, S(t), and an accurate approximation for the long-time limit, S=lim(t→∞)S(t). Unlike traditional analyses on a nongrowing domain, where S≡0, we show that domain growth leads to a very different situation where S can be positive. The theoretical tools developed and validated in this study allow us to distinguish between situations where the diffusive population reaches the moving boundary at x=L(t) from other situations where the diffusive population never reaches the moving boundary at x=L(t). Making this distinction is relevant to certain applications in developmental biology, such as the development of the enteric nervous system (ENS). All theoretical predictions are verified by implementing a discrete stochastic model.

Exact solutions of linear reaction-diffusion processes on a uniformly growing domain: criteria for successful colonization

Many processes during embryonic development involve transport and reaction of molecules, or transport and proliferation of cells, within growing tissues. Mathematical models of such processes usually take the form of a reaction-diffusion partial differential equation (PDE) on a growing domain. Previous analyses of such models have mainly involved solving the PDEs numerically. Here, we present a framework for calculating the exact solution of a linear reaction-diffusion PDE on a growing domain. We derive an exact solution for a general class of one-dimensional linear reaction—diffusion process on 0<x<L(t), where L(t) is the length of the growing domain. Comparing our exact solutions with numerical approximations confirms the veracity of the method. Furthermore, our examples illustrate a delicate interplay between: (i) the rate at which the domain elongates, (ii) the diffusivity associated with the spreading density profile, (iii) the reaction rate, and (iv) the initial condition. Altering the balance between these four features leads to different outcomes in terms of whether an initial profile, located near x = 0, eventually overcomes the domain growth and colonizes the entire length of the domain by reaching the boundary where x = L(t).

Spatial and temporal dynamics of cell generations within an invasion wave: a link to cell lineage tracing

Mathematical models of a cell invasion wave have included both continuum partial differential equation (PDE) approaches and discrete agent-based cellular automata (CA) approaches. Here we are interested in modelling the spatial and temporal dynamics of the number of divisions (generation number) that cells have undergone by any time point within an invasion wave. In the CA framework this is performed from agent lineage tracings, while in the PDE approach a multi-species generalized Fisher equation is derived for the cell density within each generation. Both paradigms exhibit qualitatively similar cell generation densities that are spatially organized, with agents of low generation number rapidly attaining a steady state (with average generation number increasing linearly with distance) behind the moving wave and with evolving high generation number at the wavefront. This regularity in the generation spatial distributions is in contrast to the highly stochastic nature of the underlying lineage dynamics of the population. In addition, we construct a method for determining the lineage tracings of all agents without labelling and tracking the agents, but through either a knowledge of the spatial distribution of the generations or the number of agents in each generation. This involves determining generation-dependent proliferation probabilities and using these to define a generation-dependent Galton-Watson (GDGW) process. Monte-Carlo simulations of the GDGW process are used to determine the individual lineage tracings. The lineages of the GDGW process are analyzed using Lorenz curves and found to be similar to outcomes generated by direct lineage tracing in CA realizations. This analysis provides the basis for a potentially useful technique for deducing cell lineage data when imaging every cell is not feasible.

Approximating spatially exclusive invasion processes

A number of biological processes, such as invasive plant species and cell migration, are composed of two key mechanisms: motility and reproduction. Due to the spatially exclusive interacting behavior of these processes a cellular automata (CA) model is specified to simulate a one-dimensional invasion process. Three (independence, Poisson, and 2D-Markov chain) approximations are considered that attempt to capture the average behavior of the CA. We show that our 2D-Markov chain approximation accurately predicts the state of the CA for a wide range of motility and reproduction rates.

Distinguishing between mechanisms of cell aggregation using pair-correlation functions

Many cell types form clumps or aggregates when cultured in vitro through a variety of mechanisms including rapid cell proliferation, chemotaxis, or direct cell-to-cell contact. In this paper we develop an agent-based model to explore the formation of aggregates in cultures where cells are initially distributed uniformly, at random, on a two-dimensional substrate. Our model includes unbiased random cell motion, together with two mechanisms which can produce cell aggregates: (i) rapid cell proliferation and (ii) a biased cell motility mechanism where cells can sense other cells within a finite range, and will tend to move towards areas with higher numbers of cells. We then introduce a pair-correlation function which allows us to quantify aspects of the spatial patterns produced by our agent-based model. In particular, these pair-correlation functions are able to detect differences between domains populated uniformly at random (i.e. at the exclusion complete spatial randomness (ECSR) state) and those where the proliferation and biased motion rules have been employed - even when such differences are not obvious to the naked eye. The pair-correlation function can also detect the emergence of a characteristic inter-aggregate distance which occurs when the biased motion mechanism is dominant, and is not observed when cell proliferation is the main mechanism of aggregate formation. This suggests that applying the pair-correlation function to experimental images of cell aggregates may provide information about the mechanism associated with observed aggregates. As a proof of concept, we perform such analysis for images of cancer cell aggregates, which are known to be associated with rapid proliferation. The results of our analysis are consistent with the predictions of the proliferation-based simulations, which supports the potential usefulness of pair correlation functions for providing insight into the mechanisms of aggregate formation.

Discrete and continuous models for tissue growth and shrinkage

The incorporation of domain growth into stochastic models of biological processes is of increasing interest to mathematical modellers and biologists alike. In many situations, especially in developmental biology, the growth of the underlying tissue domain plays an important role in the redistribution of particles (be they cells or molecules) which may move and react atop the domain. Although such processes have largely been modelled using deterministic, continuum models there is an increasing appetite for individual-based stochastic models which can capture the fine details of the biological movement processes which are being elucidated by modern experimental techniques, and also incorporate the inherent stochasticity of such systems. In this work we study a simple stochastic model of domain growth. From a basic version of this model, Hywood et al. (2013) were able to derive a Fokker-Plank equation (FPE) (in this case an advection-diffusion partial differential equation on a growing domain) which describes the evolution of the probability density of some tracer particles on the domain. We extend their work so that a variety of different domain growth mechanisms can be incorporated and demonstrate a good agreement between the mean tracer density and the solution of the FPE in each case. In addition we incorporate domain shrinkage (via element death) into our individual-level model and demonstrate that we are able to derive coefficients for the FPE in this case as well. For situations in which the drift and diffusion coefficients are not readily available we introduce a numerical coefficient estimation approach and demonstrate the accuracy of this approach by comparing it with situations in which an analytical solution is obtainable.

Cell lineage tracing in the developing enteric nervous system: superstars revealed by experiment and simulation

Cell lineage tracing is a powerful tool for understanding how proliferation and differentiation of individual cells contribute to population behaviour. In the developing enteric nervous system (ENS), enteric neural crest (ENC) cells move and undergo massive population expansion by cell division within self-growing mesenchymal tissue. We show that single ENC cells labelled to follow clonality in the intestine reveal extraordinary and unpredictable variation in number and position of descendant cells, even though ENS development is highly predictable at the population level. We use an agent-based model to simulate ENC colonization and obtain agent lineage tracing data, which we analyse using econometric data analysis tools. In all realizations, a small proportion of identical initial agents accounts for a substantial proportion of the total final agent population. We term these individuals superstars. Their existence is consistent across individual realizations and is robust to changes in model parameters. This inequality of outcome is amplified at elevated proliferation rate. The experiments and model suggest that stochastic competition for resources is an important concept when understanding biological processes which feature high levels of cell proliferation. The results have implications for cell-fate processes in the ENS.

Modeling biological tissue growth: discrete to continuum representations

There is much interest in building deterministic continuum models from discrete agent-based models governed by local stochastic rules where an agent represents a biological cell. In developmental biology, cells are able to move and undergo cell division on and within growing tissues. A growing tissue is itself made up of cells which undergo cell division, thereby providing a significant transport mechanism for other cells within it. We develop a discrete agent-based model where domain agents represent tissue cells. Each agent has the ability to undergo a proliferation event whereby an additional domain agent is incorporated into the lattice. If a probability distribution describes the waiting times between proliferation events for an individual agent, then the total length of the domain is a random variable. The average behavior of these stochastically proliferating agents defining the growing lattice is determined in terms of a Fokker-Planck equation, with an advection and diffusion term. The diffusion term differs from the one obtained Landman and Binder [J. Theor. Biol. 259, 541 (2009)] when the rate of growth of the domain is specified, but the choice of agents is random. This discrepancy is reconciled by determining a discrete-time master equation for this process and an associated asymmetric nonexclusion random walk, together with consideration of synchronous and asynchronous updating schemes. All theoretical results are confirmed with numerical simulations. This study furthers our understanding of the relationship between agent-based rules, their implementation, and their associated partial differential equations. Since tissue growth is a significant cellular transport mechanism during embryonic growth, it is important to use the correct partial differential equation description when combining with other cellular functions.

Quantifying spatial structure in experimental observations and agent-based simulations using pair-correlation functions

We define a pair-correlation function that can be used to characterize spatiotemporal patterning in experimental images and snapshots from discrete simulations. Unlike previous pair-correlation functions, the pair-correlation functions developed here depend on the location and size of objects. The pair-correlation function can be used to indicate complete spatial randomness, aggregation, or segregation over a range of length scales, and quantifies spatial structures such as the shape, size, and distribution of clusters. Comparing pair-correlation data for various experimental and simulation images illustrates their potential use as a summary statistic for calibrating discrete models of various physical processes.

Macroscopic limits of individual-based models for motile cell populations with volume exclusion

Partial differential equation models are ubiquitous in studies of motile cell populations, giving a phenomenological description of events which can be analyzed and simulated using a wide range of existing tools. However, these models are seldom derived from individual cell behaviors and so it is difficult to accurately include biological hypotheses on this spatial scale. Moreover, studies which do attempt to link individual- and population-level behavior generally employ lattice-based frameworks in which the artifacts of lattice choice at the population level are unclear. In this work we derive limiting population-level descriptions of a motile cell population from an off-lattice, individual-based model (IBM) and investigate the effects of volume exclusion on the population-level dynamics. While motility with excluded volume in on-lattice IBMs can be accurately described by Fickian diffusion, we demonstrate that this is not the case off lattice. We show that the balance between two key parameters in the IBM (the distance moved in one step and the radius of an individual) determines whether volume exclusion results in enhanced or slowed diffusion. The magnitude of this effect is shown to increase with the number of cells and the rate of their movement. The method we describe is extendable to higher-dimensional and more complex systems and thereby provides a framework for deriving biologically realistic, continuum descriptions of motile populations.

Going from microscopic to macroscopic on nonuniform growing domains

Throughout development, chemical cues are employed to guide the functional specification of underlying tissues while the spatiotemporal distributions of such chemicals can be influenced by the growth of the tissue itself. These chemicals, termed morphogens, are often modeled using partial differential equations (PDEs). The connection between discrete stochastic and deterministic continuum models of particle migration on growing domains was elucidated by Baker, Yates, and Erban [Bull. Math. Biol. 72, 719 (2010)] in which the migration of individual particles was modeled as an on-lattice position-jump process. We build on this work by incorporating a more physically reasonable description of domain growth. Instead of allowing underlying lattice elements to instantaneously double in size and divide, we allow incremental element growth and splitting upon reaching a predefined threshold size. Such a description of domain growth necessitates a nonuniform partition of the domain. We first demonstrate that an individual-based stochastic model for particle diffusion on such a nonuniform domain partition is equivalent to a PDE model of the same phenomenon on a nongrowing domain, providing the transition rates (which we derive) are chosen correctly and we partition the domain in the correct manner. We extend this analysis to the case where the domain is allowed to change in size, altering the transition rates as necessary. Through application of the master equation formalism we derive a PDE for particle density on this growing domain and corroborate our findings with numerical simulations.

Generalized index for spatial data sets as a measure of complete spatial randomness

Spatial data sets, generated from a wide range of physical systems can be analyzed by counting the number of objects in a set of bins. Previous work has been limited to equal-sized bins, which are inappropriate for some domains (e.g., circular). We consider a nonequal size bin configuration whereby overlapping or nonoverlapping bins cover the domain. A generalized index, defined in terms of a variance between bin counts, is developed to indicate whether or not a spatial data set, generated from exclusion or nonexclusion processes, is at the complete spatial randomness (CSR) state. Limiting values of the index are determined. Using examples, we investigate trends in the generalized index as a function of density and compare the results with those using equal size bins. The smallest bin size must be much larger than the mean size of the objects. We can determine whether a spatial data set is at the CSR state or not by comparing the values of a generalized index for different bin configurations-the values will be approximately the same if the data is at the CSR state, while the values will differ if the data set is not at the CSR state. In general, the generalized index is lower than the limiting value of the index, since objects do not have access to the entire region due to blocking by other objects. These methods are applied to two applications: (i) spatial data sets generated from a cellular automata model of cell aggregation in the enteric nervous system and (ii) a known plant data distribution.

Spatial analysis of multi-species exclusion processes: application to neural crest cell migration in the embryonic gut

Hindbrain (vagal) neural crest cells become relatively uniformly distributed along the embryonic intestine during the rostral to caudal colonization wave which forms the enteric nervous system (ENS). When vagal neural crest cells are labeled before migration in avian embryos by in ovo electroporation, the distribution of labeled neural crest cells in the ENS varies vastly. In some cases, the labeled neural crest cells appear evenly distributed and interspersed with unlabeled neural crest cells along the entire intestine. However, in most specimens, labeled cells occur in relatively discrete patches of varying position, area, and cell number. To determine reasons for these differences, we use a discrete cellular automata (CA) model incorporating the underlying cellular processes of neural crest cell movement and proliferation on a growing domain, representing the elongation of the intestine during development. We use multi-species CA agents corresponding to labeled and unlabeled neural crest cells. The spatial distributions of the CA agents are quantified in terms of an index. This investigation suggests that (i) the percentage of the initial neural crest cell population that is labeled and (ii) the ratio of cell proliferation to motility are the two key parameters producing the extreme differences in spatial distributions observed in avian embryos.

Quantifying evenly distributed states in exclusion and nonexclusion processes

Spatial-point data sets, generated from a wide range of physical systems and mathematical models, can be analyzed by counting the number of objects in equally sized bins. We find that the bin counts are related to the Pólya distribution. New measures are developed which indicate whether or not a spatial data set, generated from an exclusion process, is at its most evenly distributed state, the complete spatial randomness (CSR) state. To this end, we define an index in terms of the variance between the bin counts. Limiting values of the index are determined when objects have access to the entire domain and when there are subregions of the domain that are inaccessible to objects. Using three case studies (Lagrangian fluid particles in chaotic laminar flows, cellular automata agents in discrete models, and biological cells within colonies), we calculate the indexes and verify that our theoretical CSR limit accurately predicts the state of the system. These measures should prove useful in many biological applications.