Mathematical modelling has a long history in the context of collective cell migration, with applications throughout development, disease and regenerative medicine. The aim of modelling in this context is to provide a framework in which to mathematically encode experimentally derived mechanistic hypotheses, and then to test and validate them to provide new insights and understanding. Traditionally, mathematical models have consisted of systems of partial differential equations that model the evolution of cell density over time, together with the dynamics of any associated biochemical signals or the underlying substrate. The various terms in the model are usually chosen to provide simplified, phenomenological descriptions of the underlying biology, and follow long-standing conventions in the field. However, with the recent development of a plethora of new experimental technologies that provide quantitative data on collective cell migration processes, we now have the opportunity to leverage statistical and machine learning tools to determine mathematical models directly from the data. This perspectives article aims to provide an overview of recently developed data-driven modelling approaches, outlining the main methodologies and the challenges involved in using them to interrogate real-world data relating to collective cell migration.

Mechanistic mathematical models of biological systems usually contain a number of unknown parameters whose values need to be estimated from available experimental data in order for the models to be validated and used to make quantitative predictions. This requires that the models are practically identifiable, that is, the values of the parameters can be confidently determined, given available data. A well-designed experiment can produce data that are much more informative for the purpose of inferring parameter values than a poorly designed experiment. It is, therefore, of great interest to optimally design experiments such that the resulting data maximise the practical identifiability of a chosen model. Experimental design is also useful for model discrimination, where we seek to distinguish between multiple distinct, competing models of the same biological system in order to determine which model better reveals insight into the underlying biological mechanisms. In many cases, an external stimulus can be used as a control input to probe the behaviour of the system. In this paper, we will explore techniques for optimally designing such a control for a given experiment, in order to maximise parameter identifiability and model discrimination, and demonstrate these techniques in the context of commonly applied ordinary differential equation models. We use a profile likelihood-based approach to assess parameter identifiability. We then show how the problem of optimal experimental design for model discrimination can be formulated as an optimal control problem, which can be solved efficiently by applying Pontryagin's Maximum Principle.

In eukaryotic cell chemotaxis, many cells extend and retract transient actin-driven protrusions at their membrane that facilitate both the detection of external chemical gradients and directional movement via cell-matrix coupling and traction. Although extensive experimental work has detailed how cellular protrusions and morphology vary under different environmental conditions, the mechanistic principles linking protrusive activity to these factors remain poorly understood. Here, we model the extension of actin-based protrusions in chemotaxis as an optimisation problem, wherein cells balance the detection of chemical gradients with the energetic cost of protrusion formation. Our model, built on the assumption of maximal energetic efficiency during migration, provides a framework that qualitatively reproduces experimentally observed patterns of protrusive activity across a range of biological systems and environmental conditions, suggesting that energetic efficiency may underpin the morphology and chemotactic behaviour of motile eukaryotic cells. Additionally, we leverage the model to generate novel predictions regarding cellular responses to other, experimentally untested environmental perturbations, providing testable hypotheses for future experimental work that may be used to validate and refine the model presented here.

Advances in experimental techniques allow the collection of high-resolution spatio-temporal data that track individual motile entities. These tracking data can be used to calibrate mathematical models describing the motility of individual entities. The challenges in calibrating models for single-agent motion derive from the intrinsic characteristics of experimental data, collected at discrete time steps and with measurement noise. We consider the motion of individual agents that can be described by velocity-jump models in one spatial dimension. These agents transition between a network of n states, in which each state is associated with a fixed velocity and fixed rates of switching to every other state. Exploiting approximate solutions to the resultant stochastic process, we develop a Bayesian inference framework to calibrate these models to discrete-time noisy data. We first demonstrate that the framework can be used to effectively recover the model parameters of data simulated from two-state and three-state models. Finally, we explore the question of model selection first using simulated data and then using experimental data tracking mRNA transport inside Drosophila neurons. Overall, our results demonstrate that the framework is effective and efficient in calibrating and selecting between velocity-jump models, and it can be applied to a range of motion processes.

Melanoma is an aggressive form of skin cancer. Survival rates are excellent if it is detected early but fall markedly if it metastasises. A key step in early tumour progression is the formation of cell clusters, which can promote metastasis. However, the mechanisms driving cell clustering, and the role of phenotypic heterogeneity in the dynamics of these clusters, remain poorly understood. In this work, we propose a system of ordinary differential equations that models cluster formation dynamics within a coagulation-fragmentation-proliferation framework. Using Bayesian inference, we fit this model to in vitro time-lapse microscopy data from two melanoma phenotypes-proliferative and invasive-to uncover the predominant mechanisms driving cluster formation and how these differ between phenotypes. Additionally, we provide preliminary insights into how clustering behaviour in co-cultures contrasts with that observed in monocultures. The model quantifies phenotypic differences in clustering dynamics: invasive cells in monoculture exhibit nearly threefold higher coagulation rates than proliferative cells, whereas proliferative cells display slightly higher proliferation rates. These differences align with known gene expression profiles. When applied to co-culture data, the model predicts hybrid coagulation behaviour of the clusters influenced by both proliferative and invasive cells but dominated by the invasive cells, and an elevated proliferation rate, suggesting a mutually beneficial effect of phenotypic heterogeneity on cell proliferation.

Interpreting data with mathematical models is an important aspect of real-world applied mathematical modeling. Very often we are interested to understand the extent to which a particular data set informs and constrains model parameters. This question is closely related to the concept of parameter identifiability, and in this article we present a series of computational exercises to introduce tools that can be used to assess parameter identifiability, estimate parameters and generate model predictions. Taking a likelihood-based approach, we show that very similar ideas and algorithms can be used to deal with a range of different mathematical modelling frameworks. The exercises and results presented in this article are supported by a suite of open access codes that can be accessed on GitHub.

Understanding the interactions between cells and the extracellular matrix (ECM) during collective cell invasion is crucial for advancements in tissue engineering, cancer therapies, and regenerative medicine. This study focuses on the roles of contact guidance and ECM remodelling in directing cell behaviour, with a particular emphasis on exploring how differences in cell phenotype impact collective cell invasion. We present a computationally tractable two-dimensional hybrid model of collective cell migration within the ECM, where cells are modelled as individual entities and collagen fibres as a continuous tensorial field. Our model incorporates random motility, contact guidance, cell-cell adhesion, volume filling, and the dynamic remodelling of collagen fibres through cellular secretion and degradation. Through a comprehensive parameter sweep, we provide valuable insights into how differences in the cell phenotype, in terms of the ability of the cell to migrate, secrete, degrade, and respond to contact guidance cues from the ECM, impacts the characteristics of collective cell invasion.

Collective cell migration plays a crucial role in numerous biological processes, including tumour growth, wound healing, and the immune response. Often, the migrating population consists of cells with various different phenotypes. This study derives a general mathematical framework for modelling cell migration in the local environment, which is coarse-grained from an underlying individual-based model that captures the dynamics of cell migration that are influenced by the phenotype of the cell, such as random movement, proliferation, phenotypic transitions, and interactions with the local environment. The resulting, flexible, and general model provides a continuum, macroscopic description of cell invasion, which represents the phenotype of the cell as a continuous variable and is much more amenable to simulation and analysis than its individual-based counterpart when considering a large number of phenotypes. We showcase the utility of the generalised framework in three biological scenarios: range expansion; cell invasion into the extracellular matrix; and T cell exhaustion. The results highlight how phenotypic structuring impacts the spatial and temporal dynamics of cell populations, demonstrating that different environmental pressures and phenotypic transition mechanisms significantly influence migration patterns, a phenomenon that would be computationally very expensive to explore using an individual-based model alone. This framework provides a versatile and robust tool for understanding the role of phenotypic heterogeneity in collective cell migration, with potential applications in optimising therapeutic strategies for diseases involving cell migration.

Random walks and related spatial stochastic models have been used in a range of application areas including animal and plant ecology, infectious disease epidemiology, developmental biology, wound healing, and oncology. Classical random walk models assume that all individuals in a population behave independently, ignoring local physical and biological interactions. This assumption simplifies the mathematical description of the population considerably, enabling continuum-limit descriptions to be derived and used in model analysis and fitting. However, interactions between individuals can have a crucial impact on population-level behaviour. In recent decades, research has increasingly been directed towards models that include interactions, including physical crowding effects and local biological processes such as adhesion, competition, dispersal, predation and adaptive directional bias. In this article, we review the progress that has been made with models of interacting individuals. We aim to provide an overview that is accessible to researchers in application areas, as well as to specialist modellers. We focus particularly on derivation of asymptotically exact or approximate continuum-limit descriptions and simplified deterministic models of mean-field behaviour and resulting spatial patterns. We provide worked examples and illustrative results of selected models. We conclude with a discussion of current areas of focus and future challenges.

Stochasticity plays a key role in many biological systems, necessitating the calibration of stochastic mathematical models to interpret associated data. For model parameters to be estimated reliably, it is typically the case that they must be structurally identifiable. Yet, while theory underlying structural identifiability analysis for deterministic differential equation models is highly developed, there are currently no tools for the general assessment of stochastic models. In this work, we present a differential algebra-based framework for the structural identifiability analysis of linear and a class of near-linear partially observed stochastic differential equation (SDE) models. Our framework is based on a deterministic recurrence relation that describes the dynamics of the statistical moments of the system of SDEs. From this relation, we iteratively form a series of necessarily satisfied equations involving only the observed moments, from which we are able to establish structurally identifiable parameter combinations. We demonstrate our framework for a suite of linear (two- and n-dimensional) and non-linear (two-dimensional) models. Most importantly, we define the notion of structural identifiability for SDE models and establish the effect of the initial condition on identifiability. We conclude with a discussion on the applicability and limitations of our approach, and potential future research directions in this understudied area.

Real-world cellular invasion processes often take place in curved geometries. Such problems are frequently simplified in models to neglect the curved geometry in favour of computational simplicity, yet doing so risks inaccuracies in any model-based predictions. To quantify the conditions under which neglecting a curved geometry is justifiable, we explore the dynamics of a system of reaction-diffusion equations (RDEs) on a two-dimensional annular geometry analytically. Defining ϵ as the ratio of the annulus thickness δ and radius r 0 we derive, through an asymptotic expansion, the conditions under which it is appropriate to ignore the domain curvature for a general system of reaction-diffusion equations. To highlight the consequences of these results, we simulate solutions to the Fisher-Kolmogorov-Petrovsky-Piskunov (Fisher-KPP) model, a paradigm nonlinear RDE typically used to model spatial invasion, on an annular geometry. Thus, we quantify the size of the deviation from an analogous simulation on the rectangle, and how this deviation changes across the width of the annulus. We further characterise the nature of the solutions through numerical simulations for different values of r 0 and δ . Our results provide insight into when it is appropriate to neglect the domain curvature in studying travelling wave behaviour in RDEs.

Cell heterogeneity plays an important role in patient responses to drug treatments. In many cancers, it is associated with poor treatment outcomes. Many modern drug combination therapies aim to exploit cell heterogeneity, but determining how to optimise responses from heterogeneous cell populations while accounting for multi-drug synergies remains a challenge. In this work, we introduce and analyse a general optimal control framework that can be used to model the treatment response of multiple cell populations that are treated with multiple drugs that mutually interact. In this framework, we model the effect of multiple drugs on the cell populations using a system of coupled semi-linear ordinary differential equations and derive general results for the optimal solutions. We then apply this framework to three canonical examples and discuss the wider question of how to relate mathematical optimality to clinically observable outcomes, introducing a systematic approach to propose qualitatively different classes of drug dosing inspired by optimal control.

Vertebrates have evolved great diversity in the number of segments dividing the trunk body, however the developmental origin of the evolvability of this trait is poorly understood. The number of segments is thought to be determined in embryogenesis as a product of morphogenesis of the pre-somitic mesoderm (PSM) and the periodicity of a molecular oscillator active within the PSM known as the segmentation clock. Here we explore whether the clock and PSM morphogenesis exhibit developmental modularity, as independent evolution of these two processes may explain the high evolvability of segment number. Using a computational model of the clock and PSM parameterised for zebrafish, we find that the clock is broadly robust to variation in morphogenetic processes such as cell ingression, motility, compaction, and cell division. We show that this robustness is in part determined by the length of the PSM and the strength of phase coupling in the clock. As previous studies report no changes to morphogenesis upon perturbing the clock, we suggest that the clock and morphogenesis of the PSM exhibit developmental modularity.

BackgroundCell migration and invasion are well-coordinated in development and disease but remain poorly understood. We previously showed that the neural crest (NC) cell migratory wavefront shares a 45-gene panel with other cell invasion phenomena. To rapidly and systematically identify critical genes, we performed a high-throughput siRNA screen and statistical and deep learning analyses to determine changes in NC- versus non-NC-derived human cell line behaviors.ResultsWe find 14 out of 45 genes significantly reduced c8161 melanoma cell migration; four of the 14 genes altered leader cell motility (BMP4, ITGB1, KCNE3, and RASGRP1). Deep learning identified marked disruptions in cell-neighbor interactions after BMP4 or RASGRP1 knockdown in c8161 cells. Recombinant proteins added to the culture media revealed five out of the 11 known secreted molecules stimulated c8161 cell migration. BMP4 knockdown severely reduced c8161 in vivo invasion in a chick embryo transplant model. Addition of BMP4 protein to the culture media of BMP4-siRNA-treated c8161 cells rescued cell migratory ability.ConclusionHigh-throughput screening and deep learning distilled a 45-gene panel to a small subset of genes critical to melanoma and warrant deeper in vivo functional analysis for their role and potential synergies in driving NC cell migration and invasion.

Neuroblastoma is a paediatric extracranial solid cancer that arises from the developing sympathetic nervous system and is characterised by an abnormal distribution of cell types in tumours compared to healthy infant tissues. In this paper, we propose a new mathematical model of cell differentiation during sympathoadrenal development. By performing Bayesian inference of the model parameters using clinical data from patient samples, we show that the model successfully accounts for the observed differences in cell type heterogeneity among healthy adrenal tissues and four common types of neuroblastomas. Using a phenotypically structured model, we show that alterations in healthy differentiation dynamics are related to cell malignancy, and tumour volume growth. We use this model to analyse the evolution of malignant traits in a tumour. Our findings suggest that normal development dynamics make the embryonic sympathetic nervous system more robust to perturbations and accumulation of malignancies, and that the diversity of differentiation dynamics found in the neuroblastoma subtypes lead to unique risk profiles for neuroblastoma relapse after treatment.

Advances in experimental techniques allow the collection of high-resolution spatio-temporal data that track individual motile entities over time. These tracking data motivate the use of mathematical models to characterise the motion observed. In this paper, we aim to describe the solutions of velocity-jump models for single-agent motion in one spatial dimension, characterised by successive Markovian transitions within a finite network of n states, each with a specified velocity and a fixed rate of switching to every other state. In particular, we focus on obtaining the solutions of the model subject to noisy, discrete-time, observations, with no direct access to the agent state. The lack of direct observation of the hidden state makes the problem of finding the exact distributions generally intractable. Therefore, we derive a series of approximations for the data distributions. We verify the accuracy of these approximations by comparing them to the empirical distributions generated through simulations of four example model structures. These comparisons confirm that the approximations are accurate given sufficiently infrequent state switching relative to the imaging frequency. The approximate distributions computed can be used to obtain fast forwards predictions, to give guidelines on experimental design, and as likelihoods for inference and model selection.

Autoimmune myocarditis, or cardiac muscle inflammation, is a rare but frequently fatal side-effect of immune checkpoint inhibitors (ICIs), a class of cancer therapies. Despite the dangers that side-effects such as these pose to patients, they are rarely, if ever, included explicitly when mechanistic mathematical modelling of cancer therapy is used for optimization of treatment. In this paper, we develop a two-compartment mathematical model which incorporates the impact of ICIs on both the heart and the tumour. Such a model can be used to inform the conditions under which autoimmune myocarditis may develop as a consequence of treatment. We use this model in an optimal control framework to design optimized dosing schedules for three types of ICI therapy that balance the positive and negative effects of treatment. We show that including the negative side-effects of ICI treatment explicitly within the mathematical framework significantly impacts the predictions for the optimized dosing schedule, thus stressing the importance of a holistic approach to optimizing cancer therapy regimens.

Mathematical modelling has traditionally relied on detailed system knowledge to construct mechanistic models. However, the advent of large-scale data collection and advances in machine learning have led to an increasing use of data-driven approaches. Recently, hybrid models have emerged that combine both paradigms: well-understood system components are modelled mechanistically, while unknown parts are inferred from data. Here, we focus on one such class: universal differential equations (UDEs), where neural networks are embedded within differential equations to approximate unknown dynamics. When fitted to data, these networks act as universal function approximators, learning missing functional components. In this work, we note that UDE identifiability, i.e. our ability to identify true system properties, can be split into parametric and functional identifiability (assessing identifiability for the mechanistic and data-driven model parts, respectively). Next, we investigate how UDE properties, such as neural network numbers and constraints, affect parametric and functional identifiability. Notably, we show that across a wide range of models, the generalisation of a fully mechanistic model to a UDE has little impact on the mechanistic components' parametric identifiability. Finally, we note that hybrid modelling through the fitting of unknown functions (as achieved by UDEs) is particularly well-suited to chemical reaction network (CRN) modelling. Here, CRNs are used in fields ranging from systems biology, chemistry, and pharmacology to epidemiology and population dynamics, making them highly relevant for study. By showcasing how CRN-based UDE models can be highly interpretable, we demonstrate that this hybrid approach is a promising avenue for future applications.

BACKGROUND:In vertebrate embryogenesis, cranial neural crest cells (CNCCs) migrate along discrete pathways. Analyses in the chick have identified key molecular candidates for the confinement of CNCC migration to stereotypical pathways as Colec12, Trail, and Dan. The effects of these factors on CNCCs in vitro are known, but how they confine migration to discrete streams in vivo remains poorly understood. Here, we propose and test several hypothetical mechanisms by which these factors confine cell streams and maintain coherent migration, simulating an expanded agent-based model for collective CNCC migration. RESULTS:Model simulations suggest that Trail enhances adhesion between CNCCs, facilitating movement towards stereotypical migratory pathways, whereas Colec12 confines CNCCs by inducing longer, branched filopodia that facilitate movement down Colec12 gradients and re-connections with streams. Moreover, we find that Trail and Colec12 facilitate the exchange of CNCCs and the formation of CNCC bridges between adjacent streams that are observed in vivo but poorly understood mechanistically. Finally, we predict that Dan increases the coherence of streams by modulating the speed of CNCCs at the leading edge of collectives to prevent escape. CONCLUSIONS:Our work highlights the importance of Trail, Colec12, and Dan in CNCC migration and predicts novel mechanisms for the confinement of CNCCs to stereotypical pathways in vivo.

Over the past decades, nonlocal models have been widely used to describe aggregation phenomena in biology, physics, engineering, and the social sciences. These are often derived as mean-field limits of attraction-repulsion agent-based models, and consist of systems of nonlocal partial differential equations. Using differential adhesion between cells as a biological case study, we introduce a novel local model of aggregation-diffusion phenomena. This system of local aggregation-diffusion equations is fourth-order, resembling thin-film or Cahn-Hilliard type equations. In this framework, cell sorting phenomena are explained through relative surface tensions between distinct cell types. The local model emerges as a limiting case of short-range interactions, providing a significant simplification of earlier nonlocal models, while preserving the same phenomenology. This simplification makes the model easier to implement numerically and more amenable to calibration to quantitative data. Additionally, we discuss recent analytical results based on the gradient-flow structure of the model, along with open problems and future research directions.