Interpreting scratch assays using pair density dynamics and approximate Bayesian computation

Modeling adhesion in stochastic and mean-field models of cell migration

Adhesion between cells plays an important role in many biological processes such as tissue morphogenesis and homeostasis, wound healing, and cancer cell metastasis. From a mathematical perspective, adhesion between multiple cell types has been previously analyzed using discrete and continuum models, including the cellular Potts models and partial differential equations (PDEs). While these models can represent certain biological situations well, cellular Potts models can be computationally expensive, and continuum models capture only the macroscopic behavior of a population of cells, ignoring stochasticity and the discrete nature of cell dynamics. Cellular automaton models allow us to address these problems and can be used for a wide variety of biological systems. In this paper we consider a cellular automaton approach and develop an on-lattice agent-based model (ABM) for cell migration and adhesion in a population composed of two cell types. By deriving and comparing the corresponding PDEs to the ABM, we demonstrate that cell aggregation and cell sorting are not possible in the PDE model. Therefore, we propose a set of discrete mean equations which better capture the behavior of the ABM in one and two dimensions.

An off-lattice discrete model to characterise filamentous yeast colony morphology

We combine an off-lattice agent-based mathematical model and experimentation to explore filamentous growth of a yeast colony. Under environmental stress, Saccharomyces cerevisiae yeast cells can transition from a bipolar (sated) to unipolar (pseudohyphal) budding mechanism, where cells elongate and bud end-to-end. This budding asymmetry yields spatially non-uniform growth, where filaments extend away from the colony centre, foraging for food. We use approximate Bayesian computation to quantify how individual cell budding mechanisms give rise to spatial patterns observed in experiments. We apply this method of parameter inference to experimental images of colonies of two strains of S. cerevisiae, in low and high nutrient environments. The colony size at the transition from sated to pseudohyphal growth, and a forking mechanism for pseudohyphal cell proliferation are the key features driving colony morphology. Simulations run with the most likely inferred parameters produce colony morphologies that closely resemble experimental results.

Inference on an interacting diffusion system with application to in vitro glioblastoma migration (publication template)

Glioblastoma multiforme is a highly aggressive form of brain cancer, with a median survival time for diagnosed patients of 15 months. Treatment of this cancer is typically a combination of radiation, chemotherapy and surgical removal of the tumour. However, the highly invasive and diffuse nature of glioblastoma makes surgical intrusions difficult, and the diffusive properties of glioblastoma are poorly understood. In this paper, we introduce a stochastic interacting particle system as a model of in vitro glioblastoma migration, along with a maximum likelihood-algorithm designed for inference using microscopy imaging data. The inference method is evaluated on in silico simulation of cancer cell migration, and then applied to a real data set. We find that the inference method performs with a high degree of accuracy on the in silico data, and achieve promising results given the in vitro data set.

Model-based inference of a dual role for HOPS in regulating guard cell vacuole fusion

Guard cell movements depend, in part, on the remodelling of vacuoles from a highly fragmented state to a fused morphology during stomata opening. Indeed, full opening of plant stomata requires vacuole fusion to occur. Fusion of vacuole membranes is a highly conserved process in eukaryotes, with key roles played by two multi-subunit complexes: HOPS (homotypic fusion and vacuolar protein sorting) and SNARE (soluble NSF attachment protein receptor). HOPS is a vacuole tethering factor that is thought to chaperone SNAREs from apposing vacuole membranes into a fusion-competent complex capable of rearranging membranes. In plants, recruitment of HOPS subunits to the tonoplast has been shown to require the presence of the phosphoinositide phosphatidylinositol 3-phosphate. However, chemically depleting this lipid induces vacuole fusion. To resolve this counter-intuitive observation regarding the role of HOPS in regulating plant vacuole morphology, we defined a quantitative model of vacuole fusion dynamics and used it to generate testable predictions about HOPS-SNARE interactions. We derived our model by using simulation-based inference to integrate prior knowledge about molecular interactions with limited, qualitative observations of emergent vacuole phenotypes. By constraining the model parameters to yield the emergent outcomes observed for stoma opening-as induced by two distinct chemical treatments-we predicted a dual role for HOPS and identified a stalled form of the SNARE complex that differs from phenomena reported in yeast. We predict that HOPS has contradictory actions at different points in the fusion signalling pathway, promoting the formation of SNARE complexes, but limiting their activity.

Model-based inference of a plant-specific dual role for HOPS in regulating guard cell vacuole fusion

Stomata are the pores on a leaf surface that regulate gas exchange. Each stoma consists of two guard cells whose movements regulate pore opening and thereby control CO2 fixation and water loss. Guard cell movements depend in part on the remodeling of vacuoles, which have been observed to change from a highly fragmented state to a fused morphology during stomata opening. This change in morphology requires a membrane fusion mechanism that responds rapidly to environmental signals, allowing plants to respond to diurnal and stress cues. With guard cell vacuoles being both large and responsive to external signals, stomata represent a unique system in which to delineate mechanisms of membrane fusion. Fusion of vacuole membranes is a highly conserved process in eukaryotes, with key roles played by two multi-subunit complexes: HOPS (homotypic fusion and vacuolar protein sorting) and SNARE (soluble NSF attachment protein receptor). HOPS is a vacuole tethering factor that is thought to chaperone SNAREs from apposing vacuole membranes into a fusion-competent complex capable of rearranging membranes. To resolve a counter-intuitive observation regarding the role of HOPS in regulating plant vacuole morphology, we derived a quantitative model of vacuole fusion dynamics and used it to generate testable predictions about HOPS-SNARE interactions. We derived our model by applying simulation-based inference to integrate prior knowledge about molecular interactions with limited, qualitative observations of emergent vacuole phenotypes. By constraining the model parameters to yield the emergent outcomes observed for stoma opening - as induced by two distinct chemical treatments - we predicted a dual role for HOPS and identified a stalled form of the SNARE complex that differs from phenomena reported in yeast. We predict that HOPS has contradictory actions at different points in the fusion signaling pathway, promoting the formation of SNARE complexes, but limiting their activity.

FitMultiCell: simulating and parameterizing computational models of multi-scale and multi-cellular processes

MOTIVATION

Biological tissues are dynamic and highly organized. Multi-scale models are helpful tools to analyse and understand the processes determining tissue dynamics. These models usually depend on parameters that need to be inferred from experimental data to achieve a quantitative understanding, to predict the response to perturbations, and to evaluate competing hypotheses. However, even advanced inference approaches such as approximate Bayesian computation (ABC) are difficult to apply due to the computational complexity of the simulation of multi-scale models. Thus, there is a need for a scalable pipeline for modeling, simulating, and parameterizing multi-scale models of multi-cellular processes.

RESULTS

Here, we present FitMultiCell, a computationally efficient and user-friendly open-source pipeline that can handle the full workflow of modeling, simulating, and parameterizing for multi-scale models of multi-cellular processes. The pipeline is modular and integrates the modeling and simulation tool Morpheus and the statistical inference tool pyABC. The easy integration of high-performance infrastructure allows to scale to computationally expensive problems. The introduction of a novel standard for the formulation of parameter inference problems for multi-scale models additionally ensures reproducibility and reusability. By applying the pipeline to multiple biological problems, we demonstrate its broad applicability, which will benefit in particular image-based systems biology.

AVAILABILITY AND IMPLEMENTATION

FitMultiCell is available open-source at https://gitlab.com/fitmulticell/fit.

A noble extended stochastic logistic model for cell proliferation with density-dependent parameters

Cell proliferation often experiences a density-dependent intrinsic proliferation rate (IPR) and negative feedback from growth-inhibiting molecules in culture media. The lack of flexible models with explanatory parameters fails to capture such a proliferation mechanism. We propose an extended logistic growth law with the density-dependent IPR and additional negative feedback. The extended parameters of the proposed model can be interpreted as density-dependent cell-cell cooperation and negative feedback on cell proliferation. Moreover, we incorporate further density regulation for flexibility in the model through environmental resistance on cells. The proposed growth law has similarities with the strong Allee model and harvesting phenomenon. We also develop the stochastic analog of the deterministic model by representing possible heterogeneity in growth-inhibiting molecules and environmental perturbation of the culture setup as correlated multiplicative and additive noises. The model provides a conditional maximum sustainable stable cell density (MSSCD) and a new fitness measure for proliferative cells. The proposed model shows superiority to the logistic law after fitting to real cell culture datasets. We illustrate both conditional MSSCD and the new cell fitness for a range of parameters. The cell density distributions reveal the chance of overproliferation, underproliferation, or decay for different parameter sets under the deterministic and stochastic setups.

Mathematical Modeling Can Advance Wound Healing Research

Significance: For over 30 years, there has been sustained interest in the development of mathematical models for investigating the complex mechanisms underlying each stage of the wound healing process. Despite the immense associated challenges, such models have helped usher in a paradigm shift in wound healing research. Recent Advances: In this article, we review contributions in the field that span epidermal, dermal, and corneal wound healing, and treatments of nonhealing wounds. The recent influence of mathematical models on biological experiments is detailed, with a focus on wound healing assays and fibroblast-populated collagen lattices. Critical Issues: We provide an overview of the field of mathematical modeling of wound healing, highlighting key advances made in recent decades, and discuss how such models have contributed to the development of improved treatment strategies and/or an enhanced understanding of the tightly regulated steps that comprise the healing process. Future Directions: We detail some of the open problems in the field that could be addressed through a combination of theoretical and/or experimental approaches. To move the field forward, we need to have a common language between scientists to facilitate cross-collaboration, which we hope this review can support by highlighting progress to date.

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.

Discrete Manhattan and Chebyshev pair correlation functions in k dimensions

Pair correlation functions provide a summary statistic which quantifies the amount of spatial correlation between objects in a spatial domain. While pair correlation functions are commonly used to quantify continuous-space point processes, the on-lattice discrete case is less studied. Recent work has brought attention to the discrete case, wherein on-lattice pair correlation functions are formed by normalizing empirical pair distances against the probability distribution of random pair distances in a lattice with Manhattan and Chebyshev metrics. These distance distributions are typically derived on an ad hoc basis as required for specific applications. Here we present a generalized approach to deriving the probability distributions of pair distances in a lattice with discrete Manhattan and Chebyshev metrics, extending the Manhattan and Chebyshev pair correlation functions to lattices in k dimensions. We also quantify the variability of the Manhattan and Chebyshev pair correlation functions, which is important to understanding the reliability and confidence of the statistic.

Predicting population extinction in lattice-based birth-death-movement models

The question of whether a population will persist or go extinct is of key interest throughout ecology and biology. Various mathematical techniques allow us to generate knowledge regarding individual behaviour, which can be analysed to obtain predictions about the ultimate survival or extinction of the population. A common model employed to describe population dynamics is the lattice-based random walk model with crowding (exclusion). This model can incorporate behaviour such as birth, death and movement, while including natural phenomena such as finite size effects. Performing sufficiently many realizations of the random walk model to extract representative population behaviour is computationally intensive. Therefore, continuum approximations of random walk models are routinely employed. However, standard continuum approximations are notoriously incapable of making accurate predictions about population extinction. Here, we develop a new continuum approximation, the state-space diffusion approximation, which explicitly accounts for population extinction. Predictions from our approximation faithfully capture the behaviour in the random walk model, and provides additional information compared to standard approximations. We examine the influence of the number of lattice sites and initial number of individuals on the long-term population behaviour, and demonstrate the reduction in computation time between the random walk model and our approximation.

Identifying density-dependent interactions in collective cell behaviour

Scratch assays are routinely used to study collective cell behaviour in vitro. Typical experimental protocols do not vary the initial density of cells, and typical mathematical modelling approaches describe cell motility and proliferation based on assumptions of linear diffusion and logistic growth. Jin et al. (Jin et al. 2016 J. Theor. Biol. 390, 136-145 (doi:10.1016/j.jtbi.2015.10.040)) find that the behaviour of cells in scratch assays is density-dependent, and show that standard modelling approaches cannot simultaneously describe data initiated across a range of initial densities. To address this limitation, we calibrate an individual-based model to scratch assay data across a large range of initial densities. Our model allows proliferation, motility, and a direction bias to depend on interactions between neighbouring cells. By considering a hierarchy of models where we systematically and sequentially remove interactions, we perform model selection analysis to identify the minimum interactions required for the model to simultaneously describe data across all initial densities. The calibrated model is able to match the experimental data across all densities using a single parameter distribution, and captures details about the spatial structure of cells. Our results provide strong evidence to suggest that motility is density-dependent in these experiments. On the other hand, we do not see the effect of crowding on proliferation in these experiments. These results are significant as they are precisely the opposite of the assumptions in standard continuum models, such as the Fisher-Kolmogorov equation and its generalizations.

Simulation and inference algorithms for stochastic biochemical reaction networks: from basic concepts to state-of-the-art

Stochasticity is a key characteristic of intracellular processes such as gene regulation and chemical signalling. Therefore, characterizing stochastic effects in biochemical systems is essential to understand the complex dynamics of living things. Mathematical idealizations of biochemically reacting systems must be able to capture stochastic phenomena. While robust theory exists to describe such stochastic models, the computational challenges in exploring these models can be a significant burden in practice since realistic models are analytically intractable. Determining the expected behaviour and variability of a stochastic biochemical reaction network requires many probabilistic simulations of its evolution. Using a biochemical reaction network model to assist in the interpretation of time-course data from a biological experiment is an even greater challenge due to the intractability of the likelihood function for determining observation probabilities. These computational challenges have been subjects of active research for over four decades. In this review, we present an accessible discussion of the major historical developments and state-of-the-art computational techniques relevant to simulation and inference problems for stochastic biochemical reaction network models. Detailed algorithms for particularly important methods are described and complemented with Matlab^® implementations. As a result, this review provides a practical and accessible introduction to computational methods for stochastic models within the life sciences community.

Pulling in models of cell migration

There are numerous biological scenarios in which populations of cells migrate in crowded environments. Typical examples include wound healing, cancer growth, and embryo development. In these crowded environments cells are able to interact with each other in a variety of ways. These include excluded-volume interactions, adhesion, repulsion, cell signaling, pushing, and pulling. One popular way to understand the behavior of a group of interacting cells is through an agent-based mathematical model. A typical aim of modellers using such representations is to elucidate how the microscopic interactions at the cell-level impact on the macroscopic behavior of the population. At the very least, such models typically incorporate volume-exclusion. The more complex cell-cell interactions listed above have also been incorporated into such models; all apart from cell-cell pulling. In this paper we consider this under-represented cell-cell interaction, in which an active cell is able to "pull" a nearby neighbor as it moves. We incorporate a variety of potential cell-cell pulling mechanisms into on- and off-lattice agent-based volume exclusion models of cell movement. For each of these agent-based models we derive a continuum partial differential equation which describes the evolution of the cells at a population level. We study the agreement between the agent-based models and the continuum, population-based models and compare and contrast a range of agent-based models (accounting for the different pulling mechanisms) with each other. We find generally good agreement between the agent-based models and the corresponding continuum models that worsens as the agent-based models become more complex. Interestingly, we observe that the partial differential equations that we derive differ significantly, depending on whether they were derived from on- or off-lattice agent-based models of pulling. This hints that it is important to employ the appropriate agent-based model when representing pulling cell-cell interactions.

Corrected pair correlation functions for environments with obstacles

Environments with immobile obstacles or void regions that inhibit and alter the motion of individuals within that environment are ubiquitous. Correlation in the location of individuals within such environments arises as a combination of the mechanisms governing individual behavior and the heterogeneous structure of the environment. Measures of spatial structure and correlation have been successfully implemented to elucidate the roles of the mechanisms underpinning the behavior of individuals. In particular, the pair correlation function has been used across biology, ecology, and physics to obtain quantitative insight into a variety of processes. However, naively applying standard pair correlation functions in the presence of obstacles may fail to detect correlation, or suggest false correlations, due to a reliance on a distance metric that does not account for obstacles. To overcome this problem, here we present an analytic expression for calculating a corrected pair correlation function for lattice-based domains containing obstacles. We demonstrate that this obstacle pair correlation function is necessary for isolating the correlation associated with the behavior of individuals, rather than the structure of the environment. Using simulations that mimic cell migration and proliferation we demonstrate that the obstacle pair correlation function recovers the short-range correlation known to be present in this process, independent of the heterogeneous structure of the environment. Further, we show that the analytic calculation of the obstacle pair correlation function derived here is significantly faster to implement than the corresponding numerical approach.

A Bayesian Sequential Learning Framework to Parameterise Continuum Models of Melanoma Invasion into Human Skin

We present a novel framework to parameterise a mathematical model of cell invasion that describes how a population of melanoma cells invades into human skin tissue. Using simple experimental data extracted from complex experimental images, we estimate three model parameters: (i) the melanoma cell proliferation rate, [Formula: see text]; (ii) the melanoma cell diffusivity, D; and (iii) [Formula: see text], a constant that determines the rate that melanoma cells degrade the skin tissue. The Bayesian sequential learning framework involves a sequence of increasingly sophisticated experimental data from: (i) a spatially uniform cell proliferation assay; (ii) a two-dimensional circular barrier assay; and (iii) a three-dimensional invasion assay. The Bayesian sequential learning approach leads to well-defined parameter estimates. In contrast, taking a naive approach that attempts to estimate all parameters from a single set of images from the same experiment fails to produce meaningful results. Overall, our approach to inference is simple-to-implement, computationally efficient, and well suited for many cell biology phenomena that can be described by low-dimensional continuum models using ordinary differential equations and partial differential equations. We anticipate that this Bayesian sequential learning framework will be relevant in other biological contexts where it is challenging to extract detailed, quantitative biological measurements from experimental images and so we must rely on using relatively simple measurements from complex images.

The impact of experimental design choices on parameter inference for models of growing cell colonies

To better understand development, repair and disease progression, it is useful to quantify the behaviour of proliferative and motile cell populations as they grow and expand to fill their local environment. Inferring parameters associated with mechanistic models of cell colony growth using quantitative data collected from carefully designed experiments provides a natural means to elucidate the relative contributions of various processes to the growth of the colony. In this work, we explore how experimental design impacts our ability to infer parameters for simple models of the growth of proliferative and motile cell populations. We adopt a Bayesian approach, which allows us to characterize the uncertainty associated with estimates of the model parameters. Our results suggest that experimental designs that incorporate initial spatial heterogeneities in cell positions facilitate parameter inference without the requirement of cell tracking, while designs that involve uniform initial placement of cells require cell tracking for accurate parameter inference. As cell tracking is an experimental bottleneck in many studies of this type, our recommendations for experimental design provide for significant potential time and cost savings in the analysis of cell colony growth.

Pair correlation functions for identifying spatial correlation in discrete domains

Identifying and quantifying spatial correlation are important aspects of studying the collective behavior of multiagent systems. Pair correlation functions (PCFs) are powerful statistical tools that can provide qualitative and quantitative information about correlation between pairs of agents. Despite the numerous PCFs defined for off-lattice domains, only a few recent studies have considered a PCF for discrete domains. Our work extends the study of spatial correlation in discrete domains by defining a new set of PCFs using two natural and intuitive definitions of distance for a square lattice: the taxicab and uniform metric. We show how these PCFs improve upon previous attempts and compare between the quantitative data acquired. We also extend our definitions of the PCF to other types of regular tessellation that have not been studied before, including hexagonal, triangular, and cuboidal. Finally, we provide a comprehensive PCF for any tessellation and metric, allowing investigation of spatial correlation in irregular lattices for which recognizing correlation is less intuitive.

Modeling persistence of motion in a crowded environment: The diffusive limit of excluding velocity-jump processes

Persistence of motion is the tendency of an object to maintain motion in a direction for short time scales without necessarily being biased in any direction in the long term. One of the most appropriate mathematical tools to study this behavior is an agent-based velocity-jump process. In the absence of agent-agent interaction, the mean-field continuum limit of the agent-based model (ABM) gives rise to the well known hyperbolic telegraph equation. When agent-agent interaction is included in the ABM, a strictly advective system of partial differential equations (PDEs) can be derived at the population level. However, no diffusive limit of the ABM has been obtained from such a model. Connecting the microscopic behavior of the ABM to a diffusive macroscopic description is desirable, since it allows the exploration of a wider range of scenarios and establishes a direct connection with commonly used statistical tools of movement analysis. In order to connect the ABM at the population level to a diffusive PDE at the population level, we consider a generalization of the agent-based velocity-jump process on a two-dimensional lattice with three forms of agent interaction. This generalization allows us to take a diffusive limit and obtain a faithful population-level description. We investigate the properties of the model at both the individual and population levels and we elucidate some of the models' key characteristic features. In particular, we show an intrinsic anisotropy inherent to the models and we find evidence of a spontaneous form of aggregation at both the micro- and macroscales.

Bayesian inference of agent-based models: a tool for studying kidney branching morphogenesis

The adult mammalian kidney has a complex, highly-branched collecting duct epithelium that arises as a ureteric bud sidebranch from an epithelial tube known as the nephric duct. Subsequent branching of the ureteric bud to form the collecting duct tree is regulated by subcellular interactions between the epithelium and a population of mesenchymal cells that surround the tips of outgrowing branches. The mesenchymal cells produce glial cell-line derived neurotrophic factor (GDNF), that binds with RET receptors on the surface of the epithelial cells to stimulate several subcellular pathways in the epithelium. Such interactions are known to be a prerequisite for normal branching development, although competing theories exist for their role in morphogenesis. Here we introduce the first agent-based model of ex vivo kidney uretic branching. Through comparison with experimental data, we show that growth factor-regulated growth mechanisms can explain early epithelial cell branching, but only if epithelial cell division depends in a switch-like way on the local growth factor concentration; cell division occurring only if the driving growth factor level exceeds a threshold. We also show how a recently-developed method, "Approximate Approximate Bayesian Computation", can be used to infer key model parameters, and reveal the dependency between the parameters controlling a growth factor-dependent growth switch. These results are consistent with a requirement for signals controlling proliferation and chemotaxis, both of which are previously identified roles for GDNF.

Quantifying the dominant growth mechanisms of dimorphic yeast using a lattice-based model

A mathematical model is presented for the growth of yeast that incorporates both dimorphic behaviour and nutrient diffusion. The budding patterns observed in the standard and pseudohyphal growth modes are represented by a bias in the direction of cell proliferation. A set of spatial indices is developed to quantify the morphology and compare the relative importance of the directional bias to nutrient concentration and diffusivity on colony shape. It is found that there are three different growth modes: uniform growth, diffusion-limited growth (DLG) and an intermediate region in which the bias determines the morphology. The dimorphic transition due to nutrient limitation is investigated by relating the directional bias to the nutrient concentration, and this is shown to replicate the behaviour observed in vivo Comparisons are made with experimental data, from which it is found that the model captures many of the observed features. Both DLG and pseudohyphal growth are found to be capable of generating observed experimental morphologies.

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.

Using approximate Bayesian computation to quantify cell-cell adhesion parameters in a cell migratory process

In this work, we implement approximate Bayesian computational methods to improve the design of a wound-healing assay used to quantify cell–cell interactions. This is important as cell–cell interactions, such as adhesion and repulsion, have been shown to play a role in cell migration. Initially, we demonstrate with a model of an unrealistic experiment that we are able to identify model parameters that describe agent motility and adhesion, given we choose appropriate summary statistics for our model data. Following this, we replace our model of an unrealistic experiment with a model representative of a practically realisable experiment. We demonstrate that, given the current (and commonly used) experimental set-up, our model parameters cannot be accurately identified using approximate Bayesian computation methods. We compare new experimental designs through simulation, and show more accurate identification of model parameters is possible by expanding the size of the domain upon which the experiment is performed, as opposed to increasing the number of experimental replicates. The results presented in this work, therefore, describe time and cost-saving alterations for a commonly performed experiment for identifying cell motility parameters. Moreover, this work will be of interest to those concerned with performing experiments that allow for the accurate identification of parameters governing cell migratory processes, especially cell migratory processes in which cell–cell adhesion or repulsion are known to play a significant role.

Cell motility is a central process in wound healing and relies on complex cell-cell interactions. A team of mathematicians led by Ruth Baker and Kit Yates at the University of Oxford utilised computer simulations to re-design wound-healing assays that efficiently identify cell motility parameters. New experimental designs through computer simulation can more accurately identify cell motility parameters by expanding the size of the domain upon which the experiment is performed, as opposed to increasing the number of experimental replicates. The results describe time and cost-saving alterations for an experimental method for evaluate complex cell-cell interactions.

Parallelization and High-Performance Computing Enables Automated Statistical Inference of Multi-scale Models

Mechanistic understanding of multi-scale biological processes, such as cell proliferation in a changing biological tissue, is readily facilitated by computational models. While tools exist to construct and simulate multi-scale models, the statistical inference of the unknown model parameters remains an open problem. Here, we present and benchmark a parallel approximate Bayesian computation sequential Monte Carlo (pABC SMC) algorithm, tailored for high-performance computing clusters. pABC SMC is fully automated and returns reliable parameter estimates and confidence intervals. By running the pABC SMC algorithm for ∼10⁶ hr, we parameterize multi-scale models that accurately describe quantitative growth curves and histological data obtained in vivo from individual tumor spheroid growth in media droplets. The models capture the hybrid deterministic-stochastic behaviors of 10⁵-10⁶ of cells growing in a 3D dynamically changing nutrient environment. The pABC SMC algorithm reliably converges to a consistent set of parameters. Our study demonstrates a proof of principle for robust, data-driven modeling of multi-scale biological systems and the feasibility of multi-scale model parameterization through statistical inference.

Identifying the necrotic zone boundary in tumour spheroids with pair-correlation functions

Automatic identification of the necrotic zone boundary is important in the assessment of treatments on in vitro tumour spheroids. This has been difficult especially when the difference in cell density between the necrotic and viable zones of a tumour spheroid is small. To help overcome this problem, we developed novel one-dimensional pair-correlation functions (PCFs) to provide quantitative estimates of the radial distance of the necrotic zone boundary from the centre of a tumour spheroid. We validate our approach on synthetic tumour spheroids in which the position of the necrotic zone boundary is known a priori It is then applied to nine real tumour spheroids imaged with light sheet-based fluorescence microscopy. PCF estimates of the necrotic zone boundary are compared with those of a human expert and an existing standard computational method.

Filling the gaps: A robust description of adhesive birth-death-movement processes

Existing continuum descriptions of discrete adhesive birth-death-movement processes provide accurate predictions of the average discrete behavior for limited parameter regimes. Here we present an alternative continuum description in terms of the dynamics of groups of contiguous occupied and vacant lattice sites. Our method provides more accurate predictions, is valid in parameter regimes that could not be described by previous continuum descriptions, and provides information about the spatial clustering of occupied sites. Furthermore, we present a simple analytic approximation of the spatial clustering of occupied sites at late time, when the system reaches its steady-state configuration.

Quantifying the effect of experimental design choices for in vitro scratch assays

Scratch assays are often used to investigate potential drug treatments for chronic wounds and cancer. Interpreting these experiments with a mathematical model allows us to estimate the cell diffusivity, D, and the cell proliferation rate, λ. However, the influence of the experimental design on the estimates of D and λ is unclear. Here we apply an approximate Bayesian computation (ABC) parameter inference method, which produces a posterior distribution of D and λ, to new sets of synthetic data, generated from an idealised mathematical model, and experimental data for a non-adhesive mesenchymal population of fibroblast cells. The posterior distribution allows us to quantify the amount of information obtained about D and λ. We investigate two types of scratch assay, as well as varying the number and timing of the experimental observations captured. Our results show that a scrape assay, involving one cell front, provides more precise estimates of D and λ, and is more computationally efficient to interpret than a wound assay, with two opposingly directed cell fronts. We find that recording two observations, after making the initial observation, is sufficient to estimate D and λ, and that the final observation time should correspond to the time taken for the cell front to move across the field of view. These results provide guidance for estimating D and λ, while simultaneously minimising the time and cost associated with performing and interpreting the experiment.

Spatial moment dynamics for collective cell movement incorporating a neighbour-dependent directional bias

The ability of cells to undergo collective movement plays a fundamental role in tissue repair, development and cancer. Interactions occurring at the level of individual cells may lead to the development of spatial structure which will affect the dynamics of migrating cells at a population level. Models that try to predict population-level behaviour often take a mean-field approach, which assumes that individuals interact with one another in proportion to their average density and ignores the presence of any small-scale spatial structure. In this work, we develop a lattice-free individual-based model (IBM) that uses random walk theory to model the stochastic interactions occurring at the scale of individual migrating cells. We incorporate a mechanism for local directional bias such that an individual's direction of movement is dependent on the degree of cell crowding in its neighbourhood. As an alternative to the mean-field approach, we also employ spatial moment theory to develop a population-level model which accounts for spatial structure and predicts how these individual-level interactions propagate to the scale of the whole population. The IBM is used to derive an equation for dynamics of the second spatial moment (the average density of pairs of cells) which incorporates the neighbour-dependent directional bias, and we solve this numerically for a spatially homogeneous case.

Estimating cell diffusivity and cell proliferation rate by interpreting IncuCyte ZOOM™ assay data using the Fisher-Kolmogorov model

Background

Standard methods for quantifying IncuCyte ZOOM™ assays involve measurements that quantify how rapidly the initially-vacant area becomes re-colonised with cells as a function of time. Unfortunately, these measurements give no insight into the details of the cellular-level mechanisms acting to close the initially-vacant area. We provide an alternative method enabling us to quantify the role of cell motility and cell proliferation separately. To achieve this we calibrate standard data available from IncuCyte ZOOM™ images to the solution of the Fisher-Kolmogorov model.

Results

The Fisher-Kolmogorov model is a reaction-diffusion equation that has been used to describe collective cell spreading driven by cell migration, characterised by a cell diffusivity, D, and carrying capacity limited proliferation with proliferation rate, λ, and carrying capacity density, K. By analysing temporal changes in cell density in several subregions located well-behind the initial position of the leading edge we estimate λ and K. Given these estimates, we then apply automatic leading edge detection algorithms to the images produced by the IncuCyte ZOOM™ assay and match this data with a numerical solution of the Fisher-Kolmogorov equation to provide an estimate of D. We demonstrate this method by applying it to interpret a suite of IncuCyte ZOOM™ assays using PC-3 prostate cancer cells and obtain estimates of D, λ and K. Comparing estimates of D, λ and K for a control assay with estimates of D, λ and K for assays where epidermal growth factor (EGF) is applied in varying concentrations confirms that EGF enhances the rate of scratch closure and that this stimulation is driven by an increase in D and λ, whereas K is relatively unaffected by EGF.

Conclusions

Our approach for estimating D, λ and K from an IncuCyte ZOOM™ assay provides more detail about cellular-level behaviour than standard methods for analysing these assays. In particular, our approach can be used to quantify the balance of cell migration and cell proliferation and, as we demonstrate, allow us to quantify how the addition of growth factors affects these processes individually.

Spectral analysis of pair-correlation bandwidth: application to cell biology images

Images from cell biology experiments often indicate the presence of cell clustering, which can provide insight into the mechanisms driving the collective cell behaviour. Pair-correlation functions provide quantitative information about the presence, or absence, of clustering in a spatial distribution of cells. This is because the pair-correlation function describes the ratio of the abundance of pairs of cells, separated by a particular distance, relative to a randomly distributed reference population. Pair-correlation functions are often presented as a kernel density estimate where the frequency of pairs of objects are grouped using a particular bandwidth (or bin width), Δ>0. The choice of bandwidth has a dramatic impact: choosing Δ too large produces a pair-correlation function that contains insufficient information, whereas choosing Δ too small produces a pair-correlation signal dominated by fluctuations. Presently, there is little guidance available regarding how to make an objective choice of Δ. We present a new technique to choose Δ by analysing the power spectrum of the discrete Fourier transform of the pair-correlation function. Using synthetic simulation data, we confirm that our approach allows us to objectively choose Δ such that the appropriately binned pair-correlation function captures known features in uniform and clustered synthetic images. We also apply our technique to images from two different cell biology assays. The first assay corresponds to an approximately uniform distribution of cells, while the second assay involves a time series of images of a cell population which forms aggregates over time. The appropriately binned pair-correlation function allows us to make quantitative inferences about the average aggregate size, as well as quantifying how the average aggregate size changes with time.