Bridging the gap between individual-based and continuum models of growing cell populations

Agent-Based and Continuum Models for Spatial Dynamics of Infection by Oncolytic Viruses

The use of oncolytic viruses as cancer treatment has received considerable attention in recent years, however the spatial dynamics of this viral infection is still poorly understood. We present here a stochastic agent-based model describing infected and uninfected cells for solid tumours, which interact with viruses in the absence of an immune response. Two kinds of movement, namely undirected random and pressure-driven movements, are considered: the continuum limit of the models is derived and a systematic comparison between the systems of partial differential equations and the individual-based model, in one and two dimensions, is carried out. In the case of undirected movement, a good agreement between agent-based simulations and the numerical and well-known analytical results for the continuum model is possible. For pressure-driven motion, instead, we observe a wide parameter range in which the infection of the agents remains confined to the center of the tumour, even though the continuum model shows traveling waves of infection; outcomes appear to be more sensitive to stochasticity and uninfected regions appear harder to invade, giving rise to irregular, unpredictable growth patterns. Our results show that the presence of spatial constraints in tumours' microenvironments limiting free expansion has a very significant impact on virotherapy. Outcomes for these tumours suggest a notable increase in variability. All these aspects can have important effects when designing individually tailored therapies where virotherapy is included.

The Impact of Phenotypic Heterogeneity on Chemotactic Self-Organisation

The capacity to aggregate through chemosensitive movement forms a paradigm of self-organisation, with examples spanning cellular and animal systems. A basic mechanism assumes a phenotypically homogeneous population that secretes its own attractant, with the well known system introduced more than five decades ago by Keller and Segel proving resolutely popular in modelling studies. The typical assumption of population phenotypic homogeneity, however, often lies at odds with the heterogeneity of natural systems, where populations may comprise distinct phenotypes that vary according to their chemotactic ability, attractant secretion, etc. To initiate an understanding into how this diversity can impact on autoaggregation, we propose a simple extension to the classical Keller and Segel model, in which the population is divided into two distinct phenotypes: those performing chemotaxis and those producing attractant. Using a combination of linear stability analysis and numerical simulations, we demonstrate that switching between these phenotypic states alters the capacity of a population to self-aggregate. Further, we show that switching based on the local environment (population density or chemoattractant level) leads to diverse patterning and provides a route through which a population can effectively curb the size and density of an aggregate. We discuss the results in the context of real world examples of chemotactic aggregation, as well as theoretical aspects of the model such as global existence and blow-up of solutions.

A global method for fast simulations of molecular dynamics in multiscale agent-based models of biological tissues

Agent-based models (ABMs) are a natural platform for capturing the multiple time and spatial scales in biological processes. However, these models are computationally expensive, especially when including molecular-level effects. The traditional approach to simulating this type of multiscale ABM is to solve a system of ordinary differential equations for the molecular events per cell. This significantly adds to the computational cost of simulations as the number of agents grows, which contributes to many ABMs being limited to around10 5 cells. We propose an approach that requires the same computational time independent of the number of agents. This speeds up the entire simulation by orders of magnitude, allowing for more thorough explorations of ABMs with even larger numbers of agents. We use two systems to show that the new method strongly agrees with the traditionally used approach. This computational strategy can be applied to a wide range of biological investigations.

Extinction of Bistable Populations is Affected by the Shape of their Initial Spatial Distribution

The question of whether biological populations survive or are eventually driven to extinction has long been examined using mathematical models. In this work, we study population survival or extinction using a stochastic, discrete lattice-based random walk model where individuals undergo movement, birth and death events. The discrete model is defined on a two-dimensional hexagonal lattice with periodic boundary conditions. A key feature of the discrete model is that crowding effects are introduced by specifying two different crowding functions that govern how local agent density influences movement events and birth/death events. The continuum limit description of the discrete model is a nonlinear reaction-diffusion equation, and we focus on crowding functions that lead to linear diffusion and a bistable source term that is often associated with the strong Allee effect. Using both the discrete and continuum modelling tools, we explore the complicated relationship between the long-term survival or extinction of the population and the initial spatial arrangement of the population. In particular, we study different spatial arrangements of initial distributions: (i) a well-mixed initial distribution where the initial density is independent of position in the domain; (ii) a vertical strip initial distribution where the initial density is independent of vertical position in the domain; and, (iii) several forms of two-dimensional initial distributions where the initial population is distributed in regions with different shapes. Our results indicate that the shape of the initial spatial distribution of the population affects extinction of bistable populations. All software required to solve the discrete and continuum models used in this work are available on GitHub .

Derivation of continuum models from discrete models of mechanical forces in cell populations

In certain discrete models of populations of biological cells, the mechanical forces between the cells are center based or vertex based on the microscopic level where each cell is individually represented. The cells are circular or spherical in a center based model and polygonal or polyhedral in a vertex based model. On a higher, macroscopic level, the time evolution of the density of the cells is described by partial differential equations (PDEs). We derive relations between the modelling on the micro and macro levels in one, two, and three dimensions by regarding the micro model as a discretization of a PDE for conservation of mass on the macro level. The forces in the micro model correspond on the macro level to a gradient of the pressure scaled by quantities depending on the cell geometry. The two levels of modelling are compared in numerical experiments in one and two dimensions.

Modeling codelivery of CD73 inhibitor and dendritic cell-based vaccines in cancer immunotherapy

Dendritic cells (DCs) are the dominant class of antigen-presenting cells in humans; therefore, a range of DC-based approaches have been established to promote an immune response against cancer cells. The efficacy of DC-based immunotherapeutic approaches is markedly affected by the immunosuppressive factors related to the tumor microenvironment, such as adenosine. In this paper, based on immunological theories and experimental data, a hybrid model is designed that offers some insights into the effects of DC-based immunotherapy combined with adenosine inhibition. The model combines an individual-based model for describing tumor-immune system interactions with a set of ordinary differential equations for adenosine modeling. Computational simulations of the proposed model clarify the conditions for the onset of a successful immune response against cancer cells. Global and local sensitivity analysis of the model highlights the importance of adenosine blockage for strengthening effector cells. The model is used to determine the most effective suppressive mechanism caused by adenosine, proper vaccination time, and the appropriate time interval between injections.

Learning differential equation models from stochastic agent-based model simulations

Agent-based models provide a flexible framework that is frequently used for modelling many biological systems, including cell migration, molecular dynamics, ecology and epidemiology. Analysis of the model dynamics can be challenging due to their inherent stochasticity and heavy computational requirements. Common approaches to the analysis of agent-based models include extensive Monte Carlo simulation of the model or the derivation of coarse-grained differential equation models to predict the expected or averaged output from the agent-based model. Both of these approaches have limitations, however, as extensive computation of complex agent-based models may be infeasible, and coarse-grained differential equation models can fail to accurately describe model dynamics in certain parameter regimes. We propose that methods from the equation learning field provide a promising, novel and unifying approach for agent-based model analysis. Equation learning is a recent field of research from data science that aims to infer differential equation models directly from data. We use this tutorial to review how methods from equation learning can be used to learn differential equation models from agent-based model simulations. We demonstrate that this framework is easy to use, requires few model simulations, and accurately predicts model dynamics in parameter regions where coarse-grained differential equation models fail to do so. We highlight these advantages through several case studies involving two agent-based models that are broadly applicable to biological phenomena: a birth-death-migration model commonly used to explore cell biology experiments and a susceptible-infected-recovered model of infectious disease spread.

Comparative analysis of continuum angiogenesis models

Although discrete approaches are increasingly employed to model biological phenomena, it remains unclear how complex, population-level behaviours in such frameworks arise from the rules used to represent interactions between individuals. Discrete-to-continuum approaches, which are used to derive systems of coarse-grained equations describing the mean-field dynamics of a microscopic model, can provide insight into such emergent behaviour. Coarse-grained models often contain nonlinear terms that depend on the microscopic rules of the discrete framework, however, and such nonlinearities can make a model difficult to mathematically analyse. By contrast, models developed using phenomenological approaches are typically easier to investigate but have a more obscure connection to the underlying microscopic system. To our knowledge, there has been little work done to compare solutions of phenomenological and coarse-grained models. Here we address this problem in the context of angiogenesis (the creation of new blood vessels from existing vasculature). We compare asymptotic solutions of a classical, phenomenological “snail-trail” model for angiogenesis to solutions of a nonlinear system of partial differential equations (PDEs) derived via a systematic coarse-graining procedure (Pillay et al. in Phys Rev E 95(1):012410, 2017. https://doi.org/10.1103/PhysRevE.95.012410). For distinguished parameter regimes corresponding to chemotaxis-dominated cell movement and low branching rates, both continuum models reduce at leading order to identical PDEs within the domain interior. Numerical and analytical results confirm that pointwise differences between solutions to the two continuum models are small if these conditions hold, and demonstrate how perturbation methods can be used to determine when a phenomenological model provides a good approximation to a more detailed coarse-grained system for the same biological process.

Invading and Receding Sharp-Fronted Travelling Waves

Biological invasion, whereby populations of motile and proliferative individuals lead to moving fronts that invade vacant regions, is routinely studied using partial differential equation models based upon the classical Fisher-KPP equation. While the Fisher-KPP model and extensions have been successfully used to model a range of invasive phenomena, including ecological and cellular invasion, an often-overlooked limitation of the Fisher-KPP model is that it cannot be used to model biological recession where the spatial extent of the population decreases with time. In this work, we study the Fisher-Stefan model, which is a generalisation of the Fisher-KPP model obtained by reformulating the Fisher-KPP model as a moving boundary problem. The nondimensional Fisher-Stefan model involves just one parameter, [Formula: see text], which relates the shape of the density front at the moving boundary to the speed of the associated travelling wave, c. Using numerical simulation, phase plane and perturbation analysis, we construct approximate solutions of the Fisher-Stefan model for both slowly invading and receding travelling waves, as well as for rapidly receding travelling waves. These approximations allow us to determine the relationship between c and [Formula: see text] so that commonly reported experimental estimates of c can be used to provide estimates of the unknown parameter [Formula: see text]. Interestingly, when we reinterpret the Fisher-KPP model as a moving boundary problem, many overlooked features of the classical Fisher-KPP phase plane take on a new interpretation since travelling waves solutions with [Formula: see text] are normally disregarded. This means that our analysis of the Fisher-Stefan model has both practical value and an inherent mathematical value.

Evaluating snail-trail frameworks for leader-follower behavior with agent-based modeling

Branched networks constitute a ubiquitous structure in biology, arising in plants, lungs, and the circulatory system; however, the mechanisms behind their creation are not well understood. A commonly used model for network morphogenesis proposes that sprouts develop through interactions between leader (tip) cells and follower (stalk) cells. In this description, tip cells emerge from existing structures, travel up chemoattractant gradients, and form new networks by guiding the movement of stalk cells. Such dynamics have been mathematically represented by continuum "snail-trail" models in which the tip cell flux contributes to the stalk cell proliferation rate. Although snail-trail models constitute a classical depiction of leader-follower behavior, their accuracy has yet to be evaluated in a rigorous quantitative setting. Here, we extend the snail-trail modeling framework to two spatial dimensions by introducing a novel multiplicative factor to the stalk cell rate equation, which corrects for neglected network creation in directions other than that of the migrating front. Our derivation of this factor demonstrates that snail-trail models are valid descriptions of cell dynamics when chemotaxis dominates cell movement. We confirm that our snail-trail model accurately predicts the dynamics of tip and stalk cells in an existing agent-based model (ABM) for network formation [Pillay et al., Phys. Rev. E 95, 012410 (2017)10.1103/PhysRevE.95.012410]. We also derive conditions for which it is appropriate to use a reduced, one-dimensional snail-trail model to analyze ABM results. Our analysis identifies key metrics for cell migration that may be used to anticipate when simple snail-trail models will accurately describe experimentally observed cell dynamics in network formation.

Population dynamics with spatial structure and an Allee effect

Population dynamics including a strong Allee effect describe the situation where long-term population survival or extinction depends on the initial population density. A simple mathematical model of an Allee effect is one where initial densities below the threshold lead to extinction, whereas initial densities above the threshold lead to survival. Mean-field models of population dynamics neglect spatial structure that can arise through short-range interactions, such as competition and dispersal. The influence of non-mean-field effects has not been studied in the presence of an Allee effect. To address this, we develop an individual-based model that incorporates both short-range interactions and an Allee effect. To explore the role of spatial structure we derive a mathematically tractable continuum approximation of the IBM in terms of the dynamics of spatial moments. In the limit of long-range interactions where the mean-field approximation holds, our modelling framework recovers the mean-field Allee threshold. We show that the Allee threshold is sensitive to spatial structure neglected by mean-field models. For example, there are cases where the mean-field model predicts extinction but the population actually survives. Through simulations we show that our new spatial moment dynamics model accurately captures the modified Allee threshold in the presence of spatial structure.

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

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

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.

A stochastic individual-based model to explore the role of spatial interactions and antigen recognition in the immune response against solid tumours

Spatial interactions between cancer and immune cells, as well as the recognition of tumour antigens by cells of the immune system, play a key role in the immune response against solid tumours. The existing mathematical models generally focus only on one of these key aspects. We present here a spatial stochastic individual-based model that explicitly captures antigen expression and recognition. In our model, each cancer cell is characterised by an antigen profile which can change over time due to either epimutations or mutations. The immune response against the cancer cells is initiated by the dendritic cells that recognise the tumour antigens and present them to the cytotoxic T cells. Consequently, T cells become activated against the tumour cells expressing such antigens. Moreover, the differences in movement between inactive and active immune cells are explicitly taken into account by the model. Computational simulations of our model clarify the conditions for the emergence of tumour clearance, dormancy or escape, and allow us to assess the impact of antigenic heterogeneity of cancer cells on the efficacy of immune action. Ultimately, our results highlight the complex interplay between spatial interactions and adaptive mechanisms that underpins the immune response against solid tumours, and suggest how this may be exploited to further develop cancer immunotherapies.