A 3D Hybrid Model for Tissue Growth: The Interplay between Cell Population and Mass Transport Dynamics

To provide theoretical guidance for the design and in vitro cultivation of bioartificial tissues, we have developed a multiscale computational model that can describe the complex interplay between cell population and mass transport dynamics that governs the growth of tissues in three-dimensional scaffolds. The model has three components: a transient partial differential equation for the simultaneous diffusion and consumption of a limiting nutrient; a cellular automaton describing cell migration, proliferation, and collision; and equations that quantify how the varying nutrient concentration modulates cell division and migration. The hybrid discrete-continuous model was parallelized and solved on a distributed-memory multicomputer to study how transport limitations affect tissue regeneration rates under conditions encountered in typical bioreactors. Simulation results show that the severity of transport limitations can be estimated by the magnitude of two dimensionless groups: the Thiele modulus and the Biot number. Key parameters including the initial seeding mode, cell migration speed, and the hydrodynamic conditions in the bioreactor are shown to affect not only the overall rate, but also the pattern of tissue growth. This study lays the groundwork for more comprehensive models that can handle mixed cell cultures, multiple nutrients and growth factors, and other cellular processes, such as cell death.     

Diffusion models for chemotaxis: a statistical analysis of noninteractive unicellular movement.

A program is developed for applying stochastic differential equations to models for chemotaxis. First a few of the experimental and theoretical models for chemotaxis both for swimming bacteria and for cells migrating along a substrate are reviewed. In physical and biological models of deterministic systems, finite difference equations are often replaced by a limiting differential equation in order to take advantage of the ease in the use of calculus. A similar but more intricate methodology is developed here for stochastic models for chemotaxis. This exposition is possible because recent work in probability theory gives ease in the use of the stochastic calculus for diffusions and broad applicability in the convergence of stochastic difference equations to a stochastic differential equation. Stochastic differential equations suggest useful data for the model and provide statistical tests. We begin with phenomenological considerations as we analyze a one-dimensional model proposed by Boyarsky, Noble, and Peterson in their study of human granulocytes. In this context, a theoretical model consists in identifying which diffusion best approximates a model for cell movement based upon theoretical considerations of cell physiology. Such a diffusion approximation theorem is presented along with discussion of the relationship between autocovariance and persistence. Both the stochastic calculus and the diffusion approximation theorem are described in one dimension. Finally, these tools are extended to multidimensional models and applied to a three-dimensional experimental setup of spherical symmetry.     

The cellular basis of cell sorting kinetics

Cell sorting is a dynamical cooperative phenomenon that is fundamental for tissue morphogenesis and tissue homeostasis. According to Steinberg's differential adhesion hypothesis, the structure of sorted cell aggregates is determined by physical characteristics of the respective tissues, the tissue surface tensions. Steinberg postulated that tissue surface tensions result from quantitative differences in intercellular adhesion. Several experiments in cell cultures as well as in developing organisms support this hypothesis.The question of how tissue surface tension might result from differential adhesion was addressed in some theoretical models. These models describe the cellular interdependence structure once the temporal evolution has stabilized. In general, these models are capable of reproducing sorted patterns. However, the model dynamics at the cellular scale are defined implicitly and are not well-justified. The precise mechanism describing how differential adhesion generates the observed sorting kinetics at the tissue level is still unclear.It is necessary to formulate the concepts of cell level kinetics explicitly. Only then it is possible to understand the temporal development at the cellular and tissue scales. Here we argue that individual cell mobility is reduced the more the cells stick to their neighbors. We translate this assumption into a precise mathematical model which belongs to the class of stochastic interacting particle systems. Analyzing this model, we are able to predict the emergent sorting behavior at the population level. We describe qualitatively the geometry of cell segregation depending on the intercellular adhesion parameters. Furthermore, we derive a functional relationship between intercellular adhesion and surface tension and highlight the role of cell mobility in the process of sorting. We show that the interaction between the cells and the boundary of a confining vessel has a major impact on the sorting geometry.     

Stochastic simulations on a model of circadian rhythm generation

Biological phenomena are often modeled by differential equations, where states of a model system are described by continuous real values. When we consider concentrations of molecules as dynamical variables for a set of biochemical reactions, we implicitly assume that numbers of the molecules are large enough so that their changes can be regarded as continuous and they are described deterministically. However, for a system with small numbers of molecules, changes in their numbers are apparently discrete and molecular noises become significant. In such cases, models with deterministic differential equations may be inappropriate, and the reactions must be described by stochastic equations. In this study, we focus a clock gene expression for a circadian rhythm generation, which is known as a system involving small numbers of molecules. Thus it is appropriate for the system to be modeled by stochastic equations and analyzed by methodologies of stochastic simulations. The interlocked feedback model proposed by Ueda et al. as a set of deterministic ordinary differential equations provides a basis of our analyses. We apply two stochastic simulation methods, namely Gillespie's direct method and the stochastic differential equation method also by Gillespie, to the interlocked feedback model. To this end, we first reformulated the original differential equations back to elementary chemical reactions. With those reactions, we simulate and analyze the dynamics of the model using two methods in order to compare them with the dynamics obtained from the original deterministic model and to characterize dynamics how they depend on the simulation methodologies.     

Macroscopic coherent structures in a stochastic neural network: from interface dynamics to coarse-grained bifurcation analysis

Journal of Mathematical Biology - We study coarse pattern formation in a cellular automaton modelling a spatially-extended stochastic neural network. The model, originally proposed by Gong and...     

Integrating intracellular dynamics using CompuCell3D and Bionetsolver: applications to multiscale modelling of cancer cell growth and invasion

In this paper we present a multiscale, individual-based simulation environment that integrates CompuCell3D for lattice-based modelling on the cellular level and Bionetsolver for intracellular modelling. CompuCell3D or CC3D provides an implementation of the lattice-based Cellular Potts Model or CPM (also known as the Glazier-Graner-Hogeweg or GGH model) and a Monte Carlo method based on the metropolis algorithm for system evolution. The integration of CC3D for cellular systems with Bionetsolver for subcellular systems enables us to develop a multiscale mathematical model and to study the evolution of cell behaviour due to the dynamics inside of the cells, capturing aspects of cell behaviour and interaction that is not possible using continuum approaches. We then apply this multiscale modelling technique to a model of cancer growth and invasion, based on a previously published model of Ramis-Conde et al. (2008) where individual cell behaviour is driven by a molecular network describing the dynamics of E-cadherin and β-catenin. In this model, which we refer to as the centre-based model, an alternative individual-based modelling technique was used, namely, a lattice-free approach. In many respects, the GGH or CPM methodology and the approach of the centre-based model have the same overall goal, that is to mimic behaviours and interactions of biological cells. Although the mathematical foundations and computational implementations of the two approaches are very different, the results of the presented simulations are compatible with each other, suggesting that by using individual-based approaches we can formulate a natural way of describing complex multi-cell, multiscale models. The ability to easily reproduce results of one modelling approach using an alternative approach is also essential from a model cross-validation standpoint and also helps to identify any modelling artefacts specific to a given computational approach.     

Efficient parametric analysis of the chemical master equation through model order reduction

ABSTRACT: BACKGROUND: Stochastic biochemical reaction networks are commonly modelled by the chemical master equation, and can be simulated as first order linear differential equations through a finite state projection. Due to the very high state space dimension of these equations, numerical simulations are computationally expensive. This is a particular problem for analysis tasks requiring repeated simulations for different parameter values. Such tasks are computationally expensive to the point of infeasibility with the chemical master equation. RESULTS: In this article, we apply parametric model order reduction techniques in order to construct accurate low-dimensional parametric models of the chemical master equation. These surrogate models can be used in various parametric analysis task such as identifiability analysis, parameter estimation, or sensitivity analysis. As biological examples, we consider two models for gene regulation networks, a bistable switch and a network displaying stochastic oscillations. CONCLUSIONS: The results show that the parametric model reduction yields efficient models of stochastic biochemical reaction networks, and that these models can be useful for systems biology applications involving parametric analysis problems such as parameter exploration, optimization, estimation or sensitivity analysis.     

A model of sequential branching in hierarchical cell fate determination

Multipotent stem or progenitor cells undergo a sequential series of binary fate decisions, which ultimately generate the diversity of differentiated cells. Efforts to understand cell fate control have focused on simple gene regulatory circuits that predict the presence of multiple stable states, bifurcations and switch-like transitions. However, existing gene network models do not explain more complex properties of cell fate dynamics such as the hierarchical branching of developmental paths. Here, we construct a generic minimal model of the genetic regulatory network controlling cell fate determination, which exhibits five elementary characteristics of cell differentiation: stability, directionality, branching, exclusivity, and promiscuous expression. We argue that a modular architecture comprising repeated network elements reproduces these features of differentiation by sequentially repressing selected modules and hence restricting the dynamics to lower dimensional subspaces of the high-dimensional state space. We implement our model both with ordinary differential equations (ODEs), to explore the role of bifurcations in producing the one-way character of differentiation, and with stochastic differential equations (SDEs), to demonstrate the effect of noise on the system. We further argue that binary cell fate decisions are prevalent in cell differentiation due to general features of the underlying dynamical system. This minimal model makes testable predictions about the structural basis for directional, discrete and diversifying cell phenotype development and thus can guide the evaluation of real gene regulatory networks that govern differentiation.     

Modeling of the role of a Bax-activation switch in the mitochondrial apoptosis decision

We performed in silico modeling of the regulatory network of mitochondrial apoptosis through which we examined the role of a Bax-activation switch in governing the mitochondrial apoptosis decision. Two distinct modeling methods were used in this article. One is continuous and deterministic, comprised of a set of ordinary differential equations. The other, carried out in a discrete manner, is based on a cellular automaton, which takes stochastic fluctuations into consideration. We focused on dynamic properties of the mitochondrial apoptosis regulatory network. The roles of Bcl-2 family proteins in cellular responses to apoptotic stimuli were examined. In our simulations, a self-amplification process of Bax-activation is indicated. Further analysis suggests that the core module of Bax-activation is bistable in both deterministic and stochastic models, and this feature is robust to noise and wide ranges of parameter variation. When coupling with Bax-polymerization, it forms a one-way-switch, which governs irreversible behaviors of Bax-activation even with attenuation of apoptotic stimulus. Together with the growing biochemical evidence, we propose a novel molecular switch mechanism embedded in the mitochondrial apoptosis regulatory network and give a plausible explanation for the all-or-none, irreversible character of mitochondrial apoptosis.     

The Eps8/IRSp53/VASP network differentially controls actin capping and bundling in filopodia formation

There is a body of literature that describes the geometry and the physics of filopodia using either stochastic models or partial differential equations and elasticity and coarse-grained theory. Comparatively, there is a paucity of models focusing on the regulation of the network of proteins that control the formation of different actin structures. Using a combination of in-vivo and in-vitro experiments together with a system of ordinary differential equations, we focused on a small number of well-characterized, interacting molecules involved in actin-dependent filopodia formation: the actin remodeler Eps8, whose capping and bundling activities are a function of its ligands, Abi-1 and IRSp53, respectively; VASP and Capping Protein (CP), which exert antagonistic functions in controlling filament elongation. The model emphasizes the essential role of complexes that contain the membrane deforming protein IRSp53, in the process of filopodia initiation. This model accurately accounted for all observations, including a seemingly paradoxical result whereby genetic removal of Eps8 reduced filopodia in HeLa, but increased them in hippocampal neurons, and generated quantitative predictions, which were experimentally verified. The model further permitted us to explain how filopodia are generated in different cellular contexts, depending on the dynamic interaction established by Eps8, IRSp53 and VASP with actin filaments, thus revealing an unexpected plasticity of the signaling network that governs the multifunctional activities of its components in the formation of filopodia.     

Habitat fragmentation and extinction thresholds on fractal landscapes

Habitat fragmentation is a potentially critical factor in determining population persistence. In this paper, we explore the effect of fragmentation when the fragmentation follows a fractal pattern. The habitat is divided into patches, each of which is suitable or unsuitable. Suitable patches are either occupied or unoccupied, and change state depending on rates of colonization and local extinction. We compare the behaviour of two models: a spatially implicit patch-occupancy (PO) model and a spatially explicit cellular automaton (CA) model. The PO model has two fixed points: extinction, and a stable equilibrium with a fixed proportion of occupied patches. Global extinction results when habitat destruction reduces the proportion of suitable patches below a critical threshold. The PO model successfully recreates the extinction patterns found in other models. We translated the PO model into a stochastic cellular automaton. Fractal arrangements of suitable and unsuitable patches were used to simulate habitat fragmentation. We found that: (i) a population on a fractal landscape can tolerate more habitat destruction than predicted by the patch-occupancy model, and (ii) the extinction threshold decreases as the fractal dimension of the landscape decreases. These effects cannot be seen in spatially implicit models. Landscape struc-ture plays a vital role in mediating the effects of habitat fragmentation on persistence.     

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.     

Mathematical model of the morphogenesis checkpoint in budding yeast

The morphogenesis checkpoint in budding yeast delays progression through the cell cycle in response to stimuli that prevent bud formation. Central to the checkpoint mechanism is Swe1 kinase: normally inactive, its activation halts cell cycle progression in G2. We propose a molecular network for Swe1 control, based on published observations of budding yeast and analogous control signals in fission yeast. The proposed Swe1 network is merged with a model of cyclin-dependent kinase regulation, converted into a set of differential equations and studied by numerical simulation. The simulations accurately reproduce the phenotypes of a dozen checkpoint mutants. Among other predictions, the model attributes a new role to Hsl1, a kinase known to play a role in Swe1 degradation: Hsl1 must also be indirectly responsible for potent inhibition of Swe1 activity. The model supports the idea that the morphogenesis checkpoint, like other checkpoints, raises the cell size threshold for progression from one phase of the cell cycle to the next.     

Stochastic models for virus and immune system dynamics

New stochastic models are developed for the dynamics of a viral infection and an immune response during the early stages of infection. The stochastic models are derived based on the dynamics of deterministic models. The simplest deterministic model is a well-known system of ordinary differential equations which consists of three populations: uninfected cells, actively infected cells, and virus particles. This basic model is extended to include some factors of the immune response related to Human Immunodeficiency Virus-1 (HIV-1) infection. For the deterministic models, the basic reproduction number, R₀, is calculated and it is shown that if R₀<1, the disease-free equilibrium is locally asymptotically stable and is globally asymptotically stable in some special cases. The new stochastic models are systems of stochastic differential equations (SDEs) and continuous-time Markov chain (CTMC) models that account for the variability in cellular reproduction and death, the infection process, the immune system activation, and viral reproduction. Two viral release strategies are considered: budding and bursting. The CTMC model is used to estimate the probability of virus extinction during the early stages of infection. Numerical simulations are carried out using parameter values applicable to HIV-1 dynamics. The stochastic models provide new insights, distinct from the basic deterministic models. For the case R₀>1, the deterministic models predict the viral infection persists in the host. But for the stochastic models, there is a positive probability of viral extinction. It is shown that the probability of a successful invasion depends on the initial viral dose, whether the immune system is activated, and whether the release strategy is bursting or budding.     

A Computational Model Predicting Disruption of Blood Vessel Development

Vascular development is a complex process regulated by dynamic biological networks that vary in topology and state across different tissues and developmental stages. Signals regulating de novo blood vessel formation (vasculogenesis) and remodeling (angiogenesis) come from a variety of biological pathways linked to endothelial cell (EC) behavior, extracellular matrix (ECM) remodeling and the local generation of chemokines and growth factors. Simulating these interactions at a systems level requires sufficient biological detail about the relevant molecular pathways and associated cellular behaviors, and tractable computational models that offset mathematical and biological complexity. Here, we describe a novel multicellular agent-based model of vasculogenesis using the CompuCell3D (http://www.compucell3d.org/) modeling environment supplemented with semi-automatic knowledgebase creation. The model incorporates vascular endothelial growth factor signals, pro- and anti-angiogenic inflammatory chemokine signals, and the plasminogen activating system of enzymes and proteases linked to ECM interactions, to simulate nascent EC organization, growth and remodeling. The model was shown to recapitulate stereotypical capillary plexus formation and structural emergence of non-coded cellular behaviors, such as a heterologous bridging phenomenon linking endothelial tip cells together during formation of polygonal endothelial cords. Molecular targets in the computational model were mapped to signatures of vascular disruption derived from in vitro chemical profiling using the EPA''s ToxCast high-throughput screening (HTS) dataset. Simulating the HTS data with the cell-agent based model of vascular development predicted adverse effects of a reference anti-angiogenic thalidomide analog, 5HPP-33, on in vitro angiogenesis with respect to both concentration-response and morphological consequences. These findings support the utility of cell agent-based models for simulating a morphogenetic series of events and for the first time demonstrate the applicability of these models for predictive toxicology.     

Spatial Stochastic Dynamics Enable Robust Cell Polarization

Although cell polarity is an essential feature of living cells, it is far from being well-understood. Using a combination of computational modeling and biological experiments we closely examine an important prototype of cell polarity: the pheromone-induced formation of the yeast polarisome. Focusing on the role of noise and spatial heterogeneity, we develop and investigate two mechanistic spatial models of polarisome formation, one deterministic and the other stochastic, and compare the contrasting predictions of these two models against experimental phenotypes of wild-type and mutant cells. We find that the stochastic model can more robustly reproduce two fundamental characteristics observed in wild-type cells: a highly polarized phenotype via a mechanism that we refer to as spatial stochastic amplification, and the ability of the polarisome to track a moving pheromone input. Moreover, we find that only the stochastic model can simultaneously reproduce these characteristics of the wild-type phenotype and the multi-polarisome phenotype of a deletion mutant of the scaffolding protein Spa2. Significantly, our analysis also demonstrates that higher levels of stochastic noise results in increased robustness of polarization to parameter variation. Furthermore, our work suggests a novel role for a polarisome protein in the stabilization of actin cables. These findings elucidate the intricate role of spatial stochastic effects in cell polarity, giving support to a cellular model where noise and spatial heterogeneity combine to achieve robust biological function.     

The impact of exclusion processes on angiogenesis models

Angiogenesis is the process by which new blood vessels form from existing vessels. During angiogenesis, tip cells migrate via diffusion and chemotaxis, new tip cells are introduced through branching, loops form via tip-to-tip and tip-to-sprout anastomosis, and a vessel network forms as endothelial cells, known as stalk cells, follow the paths of tip cells (a process known as the snail-trail). Using a mean-field approximation, we systematically derive one-dimensional non-linear continuum models from a lattice-based cellular automaton model of angiogenesis in the corneal assay, explicitly accounting for cell volume. We compare our continuum models and a well-known phenomenological snail-trail model that is linear in the diffusive, chemotactic and branching terms, with averaged cellular automaton simulation results to distinguish macroscale volume exclusion effects and determine whether linear models can capture them. We conclude that, in general, both linear and non-linear models can be used at low cell densities when single or multi-species exclusion effects are negligible at the macroscale. When cell densities increase, our non-linear model should be used to capture non-linear tip cell behavior that occurs when single-species exclusion effects are pronounced, and alternative models should be derived for non-negligible multi-species exclusion effects.