Addressing Interdisciplinary Difficulties in Developmental Biology/Mathematical Collaborations: A Neural Crest Example.

Mathematical modeling can allow insight into the biological processes that can be difficult to access by conventional biological means alone. Such projects are becoming increasingly attractive with the appearance of faster and more powerful quantitative techniques in both biological data acquisition and data storage, manipulation, and presentation. However, as is frequent in interdisciplinary research, the main hurdles are not within the mindset and techniques of each discipline but are usually encountered in attempting to meld the different disciplines together. Based upon our experience in applying mathematical methods to investigate how neural crest cells interact to form the enteric nervous system, we present our views on how to pursue biomathematical modeling projects, what difficulties to expect, and how to overcome, or at least survive, these hurdles. The main advice being: persevere.

Distinguishing cell shoving mechanisms.

Motivated by in vitro time-lapse images of ovarian cancer spheroids inducing mesothelial cell clearance, the traditional agent-based model of cell migration, based on simple volume exclusion, was extended to include the possibility that a cell seeking to move into an occupied location may push the resident cell, and any cells neighbouring it, out of the way to occupy that location. In traditional discrete models of motile cells with volume exclusion such a move would be aborted. We introduce a new shoving mechanism which allows cells to choose the direction to shove cells that expends the least amount of shoving effort (to account for the likely resistance of cells to being pushed). We call this motility rule 'smart shoving'. We examine whether agent-based simulations of different shoving mechanisms can be distinguished on the basis of single realisations and averages over many realisations. We emphasise the difficulty in distinguishing cell mechanisms from cellular automata simulations based on snap-shots of cell distributions, site-occupancy averages and the evolution of the number of cells of each species averaged over many realisations. This difficulty suggests the need for higher resolution cell tracking.

Stochastic clonal expansion of "superstars" enhances the reserve capacity of enteric nervous system precursor cells.

We quantified cell population increase in the quail embryo enteric nervous system (ENS) from E2.5 (about 1500 cells) to E12 (about 8 million cells). We then probed ENS proliferative capacity by grafting to the chorio-allantoic membrane large (600 cells) and small (40 cells) populations of enteric neural crest (ENC) cells with aneural gut. This demonstrated that ENC cells show an extremely high capacity to regulate their proliferation while forming the ENS. Previous mathematical models and clonal label experiments revealed that a few dominant ENS "superstar" cell clones emerge but most clones are small. The model implied that "superstars" arise stochastically, but the same outcome could arise if "superstars" were pre-determined. We investigated these two modes mathematically and by grafting experiments with large and small numbers of ENCs, each including one EGFP-labelled ENC cell. The stochastic model predicts that the frequency of "superstar" detection increases as the ENC population decreases, the pre-determined model does not. Experimentally, as predicted by the stochastic model, the frequency of "superstar" detection increased with small ENC cell number. We conclude that ENS "superstar" clones achieve this status stochastically. Clonal dominance implies that clonal diversity is greatly reduced and in this case, somatic mutations may affect the phenotype. We suggest that somatic mutations coupled with loss of clonal diversity may contribute to variable penetrance and expressivity in individuals with genetically identical ENS pathologies.

A Multicellular Model of Intestinal Crypt Buckling and Fission.

Crypt fission is an in vivo tissue deformation process that is involved in both intestinal homeostasis and colorectal tumourigenesis. Despite its importance, the mechanics underlying crypt fission are currently poorly understood. Recent experimental development of organoids, organ-like buds cultured from crypt stem cells in vitro, has shown promise in shedding light on crypt fission. Drawing inspiration from observations of organoid growth and fission in vivo, we develop a computational model of a deformable epithelial tissue layer. Results from in silico experiments show the stiffness of cells and the proportions of cell subpopulations affect the nature of deformation in the epithelial layer. In particular, we find that increasing the proportion of stiffer cells in the layer increases the likelihood of crypt fission occurring. This is in agreement with and helps explain recent experimental work.

What is the optimal distribution of myelin along a single axon?

The myelin sheath that insulates some axons in the central nervous system allows for faster signal conduction. Previously, axons were thought to be either unmyelinated or fully myelinated. Recent experimental work has discovered a new pattern of myelination (intermittent myelination) along axons in the mouse brain, in which long unmyelinated axon segments are followed by myelinated segments of comparable length. We use a computational model to explore how myelin distribution (in particular intermittent myelination) affects conduction velocity. We find that although fully myelinated axons minimize conduction velocity, varying the spatial distribution of a fixed amount of myelin along a partially myelinated axon leads to considerable variation in the conduction velocity for action potentials. Whether sodium ion channel number or sodium ion channel density is held constant as the area of the unmyelinated segments increases has a strong influence on the optimal pattern of myelin and the conduction velocity.

Differential Clonal Expansion in an Invading Cell Population: Clonal Advantage or Dumb Luck?

In neoplastic cell growth, clones and subclones are variable both in size and mutational spectrum. The largest of these clones are believed to represent those cells with mutations that make them the most "fit," in a Darwinian sense, for expansion in their microenvironment. Thus, the degree of quantitative clonal expansion is regarded as being determined by innate qualitative differences between the cells that originate each clone. Here, using a combination of mathematical modelling and clonal labelling experiments applied to the developmental model system of the forming enteric nervous system, we describe how cells which are qualitatively identical may consistently produce clones of dramatically different sizes: most clones are very small while a few clones we term "superstars" contribute most of the cells to the final population. The basis of this is minor stochastic variations ("luck") in the timing and direction of movement and proliferation of individual cells, which builds a local advantage for daughter cells that is cumulative. This has potentially important consequences. In cancers, especially before strongly selective cytotoxic therapy, the assumption that the largest clones must be the cells with deterministic proliferative ability may not always hold true. In development, the gradual loss of clonal diversity as "superstars" take over the population may erode the resilience of the system to somatic mutations, which may have occurred early in clonal growth.

Evidence for Cooperative Selection of Axons for Myelination by Adjacent Oligodendrocytes in the Optic Nerve.

The cellular mechanisms that regulate the topographic arrangement of myelin internodes along axons remain largely uncharacterized. Recent clonal analysis of oligodendrocyte morphologies in the mouse optic nerve revealed that adjacent oligodendrocytes frequently formed adjacent internodes on one or more axons in common, whereas oligodendrocytes in the optic nerve were never observed to myelinate the same axon more than once. By modelling the process of axonal selection at the single cell level, we demonstrate that internode length and primary process length constrain the capacity of oligodendrocytes to myelinate the same axon more than once. On the other hand, probabilistic analysis reveals that the observed juxtaposition of myelin internodes among common sets of axons by adjacent oligodendrocytes is highly unlikely to occur by chance. Our analysis may reveal a hitherto unknown level of communication between adjacent oligodendrocytes in the selection of axons for myelination. Together, our analyses provide novel insights into the mechanisms that define the spatial organization of myelin internodes within white matter at the single cell level.

Is cell migration or proliferation dominant in the formation of linear arrays of oligodendrocytes?

Oligodendrocytes are the myelin-producing cells of the central nervous system that are responsible for electrically insulating axons to speed the propagation of electrical impulses. A striking feature of oligodendrocyte development within white matter is that the cell bodies of many oligodendrocyte progenitor cells become organised into discrete linear arrays of three or more cells before they differentiate into myelin-producing oligodendrocytes. These linear arrays align parallel to the direction of the axons within white matter tracts and are believed to play an important role in the co-ordination of myelination. Guided by experimental data on the abundance and composition of linear arrays in the corpus callosum of the postnatal mouse brain, we construct discrete and continuous models of linear array generation to specifically investigate the relative influence of cell migration, proliferation, differentiation and death of oligodendroglia upon the genesis of linear arrays during early postnatal development. We demonstrate that only models that incorporate significant cell migration can replicate all of the experimental observations on number of arrays, number of cells in arrays and total cell count of oligodendroglia within a given area of the corpus callosum. These models are also necessary to accurately reflect experimental data on the abundance of linear arrays composed of oligodendrocytes that derive from progenitors of different clonal origins.

Incomplete penetrance: The role of stochasticity in developmental cell colonization.

Cell colonization during embryonic development involves cells migrating and proliferating over growing tissues. Unsuccessful colonization, resulting from genetic causes, can result in various birth defects. However not all individuals with the same mutation show the disease. This is termed incomplete penetrance, and it even extends to discordancy in monozygotic (identical) twins. A one-dimensional agent-based model of cell migration and proliferation within a growing tissue is presented, where the position of every cell is recorded at any time. We develop a new model that approximates this agent-based process - rather than requiring the precise configuration of cells within the tissue, the new model records the total number of cells, the position of the most advanced cell, and then invokes an approximation for how the cells are distributed. The probability mass function (PMF) for the most advanced cell is obtained for both the agent-based model and its approximation. The two PMFs compare extremely well, but using the approximation is computationally faster. Success or failure of colonization is probabilistic. For example for sufficiently high proliferation rate the colonization is assured. However, if the proliferation rate is sufficiently low, there will be a lower, say 50%, chance of success. These results provide insights into the puzzle of incomplete penetrance of a disease phenotype, especially in monozygotic twins. Indeed, stochastic cell behavior (amplified by disease-causing mutations) within the colonization process may play a key role in incomplete penetrance, rather than differences in genes, their expression or environmental conditions.

An in silico agent-based model demonstrates Reelin function in directing lamination of neurons during cortical development.

The characteristic six-layered appearance of the neocortex arises from the correct positioning of pyramidal neurons during development and alterations in this process can cause intellectual disabilities and developmental delay. Malformations in cortical development arise when neurons either fail to migrate properly from the germinal zones or fail to cease migration in the correct laminar position within the cortical plate. The Reelin signalling pathway is vital for correct neuronal positioning as loss of Reelin leads to a partially inverted cortex. The precise biological function of Reelin remains controversial and debate surrounds its role as a chemoattractant or stop signal for migrating neurons. To investigate this further we developed an in silico agent-based model of cortical layer formation. Using this model we tested four biologically plausible hypotheses for neuron motility and four biologically plausible hypotheses for the loss of neuron motility (conversion from migration). A matrix of 16 combinations of motility and conversion rules was applied against the known structure of mouse cortical layers in the wild-type cortex, the Reelin-null mutant, the Dab1-null mutant and a conditional Dab1 mutant. Using this approach, many combinations of motility and conversion mechanisms can be rejected. For example, the model does not support Reelin acting as a repelling or as a stopping signal. In contrast, the study lends very strong support to the notion that the glycoprotein Reelin acts as a chemoattractant for neurons. Furthermore, the most viable proposition for the conversion mechanism is one in which conversion is affected by a motile neuron sensing in the near vicinity neurons that have already converted. Therefore, this model helps elucidate the function of Reelin during neuronal migration and cortical development.

Erosion by a one-dimensional random walk.

We consider a model introduced by Baker et al. [Phys. Rev. E 88, 042113 (2013)] of a single lattice random walker moving on a domain of allowed sites, surrounded by blocked sites. The walker enlarges the allowed domain by eroding the boundary at its random encounters with blocked boundary sites: attempts to step onto blocked sites succeed with a given probability and convert these sites to allowed sites. The model interpolates continuously between the Pólya random walker on the one-dimensional lattice and a "blind" walker who attempts freely, but always aborts, moves to blocked sites. We obtain some exact results about the walker's location and the rate of erosion.

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.

A dynamical model of tumour immunotherapy.

A coupled ordinary differential equation model of tumour-immune dynamics is presented and analysed. The model accounts for biological and clinical factors which regulate the interaction rates of cytotoxic T lymphocytes on the surface of the tumour mass. A phase plane analysis demonstrates that competition between tumour cells and lymphocytes can result in tumour eradication, perpetual oscillations, or unbounded solutions. To investigate the dependence of the dynamic behaviour on model parameters, the equations are solved analytically and conditions for unbounded versus bounded solutions are discussed. An analytic characterisation of the basin of attraction for oscillatory orbits is given. It is also shown that the tumour shape, characterised by a surface area to volume scaling factor, influences the size of the basin, with significant consequences for therapy design. The findings reveal that the tumour volume must surpass a threshold size that depends on lymphocyte parameters for the cancer to be completely eliminated. A semi-analytic procedure to calculate oscillation periods and determine their sensitivity to model parameters is also presented. Numerical results show that the period of oscillations exhibits notable nonlinear dependence on biologically relevant conditions.

Interacting motile agents: taking a mean-field approach beyond monomers and nearest-neighbor steps.

We consider a discrete agent-based model on a one-dimensional lattice, where each agent occupies L sites and attempts movements over a distance of d lattice sites. Agents obey a strict simple exclusion rule. A discrete-time master equation is derived using a mean-field approximation and careful probability arguments. In the continuum limit, nonlinear diffusion equations that describe the average agent occupancy are obtained. Averaged discrete simulation data are generated and shown to compare very well with the solution to the derived nonlinear diffusion equations. This framework allows us to approach a lattice-free result using all the advantages of lattice methods. Since different cell types have different shapes and speeds of movement, this work offers insight into population-level behavior of collective cellular motion.

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.

Analytic study of three-dimensional single cell migration with and without proteolytic enzymes.

Cell motility is a fundamental physiological process that regulates cellular fate in healthy and diseased systems. Cells cultured in 3D environments often exhibit biphasic dependence of migration speed with cell adhesion. Much is not understood about this very common behavior. A phenomenological model for 3D single-cell migration that exhibits biphasic behavior and highlights the important role of steric hindrance is developed and studied analytically. Changes in the biphasic behavior in the presence of proteolytic enzymes are investigated. Our methods produce a framework to determine analytic formulae for the mean cell speed, allowing general statements in terms of parameters to be explored, which will be useful when interpreting future experimental results. Our formula for mean cell speed as a function of ligand concentration generalizes and extends previous computational models that have shown good agreement with in vitro experiments.

Random walks in nonuniform environments with local dynamic interactions.

We consider a class of lattice random walk models in which the random walker is initially confined to a finite connected set of allowed sites but has the opportunity to enlarge this set by colliding with its boundaries, each such collision having a given probability of breaking through. The model is motivated by an analogy to cell motility in tissue, where motile cells have the ability to remodel extracellular matrix, but is presented here as a generic model for stochastic erosion. For the one-dimensional case, we report some exact analytic results, some mean-field type analytic approximate results and simulations. We compute exactly the mean and variance of the time taken to enlarge the interval from a single site to a given size. The problem of determining the statistics of the interval length and the walker's position at a given time is more difficult and we report several interesting observations from simulations. Our simulations include the case in which the initial interval length is random and the case in which the initial state of the lattice is a random mixture of allowed and forbidden sites, with the walker placed at random on an allowed site. To illustrate the extension of these ideas to higher-dimensional systems, we consider the erosion of the simple cubic lattice commencing from a single site and report simulations of measures of cluster size and shape and the mean-square displacement of the walker.

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.

Mathematical modelling of gubernaculum during inguino-scrotal migration shows limb bud characteristics.

PURPOSE: The gubernaculum is postulated to grow like an embryonic limb bud during inguinoscrotal descent in rodents. Recently, modelling of limb bud growth suggests the undifferentiated, distal "progress zone" provides molecular morphogenic signals, rather than cell division, as previously thought. We aimed to develop a mathematical gubernacular growth model, hypothesising that it would mimic limb buds through evolutionary conservation. METHODS: Histology was done on Sprague-Dawley rats (day 2, 8; n=7/group) to determine gubernacular length, width, cell density in distal growth centre, middle and proximal cremaster muscle. Analysis of measurements enabled gubernacular growth modelling under variable growth centre sizes/densities, assuming no apoptosis. RESULTS: Modelling found that gubernacular growth occurred mostly within cremaster muscle, rather than primarily in the undifferentiated mesenchymal tip, despite its higher mitotic rate. The growth centre accounted for ≤ 10% of total gubernacular enlargement/elongation. CONCLUSIONS: These results suggest the gubernaculum elongates by proliferation throughout cremaster muscle like a limb bud. The distal undifferentiated tip may provide signalling molecules for growth, which could be a fruitful source for causes of failed migration/elongation in cryptorchidism.

Simple rules for a "simple" nervous system? Molecular and biomathematical approaches to enteric nervous system formation and malformation.

We review morphogenesis of the enteric nervous system from migratory neural crest cells, and defects of this process such as Hirschsprung disease, centering on cell motility and assembly, and cell adhesion and extracellular matrix molecules, along with cell proliferation and growth factors. We then review continuum and agent-based (cellular automata) models with rules of cell movement and logistical proliferation. Both movement and proliferation at the individual cell level are modeled with stochastic components from which stereotyped outcomes emerge at the population level. These models reproduced the wave-like colonization of the intestine by enteric neural crest cells, and several new properties emerged, such as colonization by frontal expansion, which were later confirmed biologically. These models predict a surprising level of clonal heterogeneity both in terms of number and distribution of daughter cells. Biologically, migrating cells form stable chains made up of unstable cells, but this is not seen in the initial model. We outline additional rules for cell differentiation into neurons, axon extension, cell-axon and cell-cell adhesions, chemotaxis and repulsion which can reproduce chain migration. After the migration stage, the cells re-arrange as a network of ganglia. Changes in cell adhesion molecules parallel this, and we describe additional rules based on Steinberg's Differential Adhesion Hypothesis, reflecting changing levels of adhesion in neural crest cells and neurons. This was able to reproduce enteric ganglionation in a model. Mouse mutants with disturbances of enteric nervous system morphogenesis are discussed, and these suggest future refinement of the models. The modeling suggests a relatively simple set of cell behavioral rules could account for complex patterns of morphogenesis. The model has allowed the proposal that Hirschsprung disease is mostly an enteric neural crest cell proliferation defect, not a defect of cell migration. In addition, the model suggests an explanations for zonal and skip segment variants of Hirschsprung disease, and also gives a novel stochastic explanation for the observed discordancy of Hirschsprung disease in identical twins.