Search

Ligand-receptor interactions are fundamental to many biological processes. For example in antibody-based immunotherapies, the dynamics of an antibody binding with its target antigen directly influence the potency and efficacy of monoclonal antibody (mAb) therapies. In this paper, we present an asymptotic analysis of an ordinary differential equation (ODE) model of bivalent antibody-antigen binding in the context of mAb cancer therapies, highlighting the added complexity associated with bivalency of the antibody. To understand what drives the complex temporal dynamics of bivalent antibody-antigen binding, we construct asymptotic approximations to the model's solutions at different timescales and antibody concentrations that are in good agreement with numerical simulations of the full model. We show how the dynamics differ between two scenarios; a region where unbound antigens are abundant, and one where the number of unbound antigens is small such that the dominant balance within the model equations changes. Of particular importance to the potency and efficacy of mAb treatments are the values of quantities such as antigen occupancy and bound antibody number. We use the results of our asymptotic analysis to approximate the long-time values of these quantities that could be combined with experimental data to facilitate parameter estimation.

The logistic equation is ubiquitous in applied mathematics as a minimal model of saturating growth. Here, we examine a generalisation of the logistic model to fitness-structured populations, motivated by examples that range from the ageing of individuals in a species to immune cell exhaustion by cancerous tissue. Through exploration of a range of concrete examples and a general analysis of polynomial fitness functions, we derive necessary and sufficient conditions for the dependence of the kinetics on structure to result in closed, low-dimensional moment equations that are exact. Further, we showcase how coarse-grained moment information can be used to elucidate the details of structured dynamics, with immediate potential for model selection and hypothesis testing.

Epigenetics plays a key role in cellular differentiation and maintaining cell identity, enabling cells to regulate their genetic activity without altering the DNA sequence. Epigenetic regulation occurs within the context of hierarchically folded chromatin, yet the interplay between the dynamics of epigenetic modifications and chromatin architecture remains poorly understood. In addition, it remains unclear what mechanisms drive the formation of rugged epigenetic patterns, characterised by alternating genomic regions enriched in activating and repressive marks. In this study, we focus on post-translational modifications of histone H3 tails, particularly H3K27me3, H3K4me3, and H3K27ac. We introduce a mesoscopic stochastic model that incorporates chromatin architecture and competition of histone-modifying enzymes into the dynamics of epigenetic modifications in small genomic loci comprising several nucleosomes. Our approach enables us to investigate the mechanisms by which epigenetic patterns form on larger scales of chromatin organisation, such as loops and domains. Through bifurcation analysis and stochastic simulations, we demonstrate that the model can reproduce uniform chromatin states (open, closed, and bivalent) and generate previously unexplored rugged profiles. Our results suggest that enzyme competition and chromatin conformations with high-frequency interactions between distant genomic loci can drive the emergence of rugged epigenetic landscapes. Additionally, we hypothesise that bivalent chromatin can act as an intermediate state, facilitating transitions between uniform and rugged landscapes. This work offers a powerful mathematical framework for understanding the dynamic interactions between chromatin architecture and epigenetic regulation, providing new insights into the formation of complex epigenetic patterns., Comment: 39 pages, graphical abstract, 8 figures, 1 table. Submitted for publication to the Journal of Theoretical Biology

Macrophages in atherosclerotic lesions exhibit a spectrum of behaviours or phenotypes. The phenotypic distribution of monocyte-derived macrophages (MDMs), its correlation with MDM lipid content, and relation to blood lipoprotein densities are not well understood. Of particular interest is the balance between low density lipoproteins (LDL) and high density lipoproteins (HDL), which carry bad and good cholesterol respectively. To address these issues, we have developed a mathematical model for early atherosclerosis in which the MDM population is structured by phenotype and lipid content. The model admits a simpler, closed subsystem whose analysis shows how lesion composition becomes more pathological as the blood density of LDL increases relative to the HDL capacity. We use asymptotic analysis to derive a power-law relationship between MDM phenotype and lipid content at steady-state. This relationship enables us to understand why, for example, lipid-laden MDMs have a more inflammatory phenotype than lipid-poor MDMs when blood LDL lipid density greatly exceeds HDL capacity. We show further that the MDM phenotype distribution always attains a local maximum, while the lipid content distribution may be unimodal, adopt a quasi-uniform profile or decrease monotonically. Pathological lesions exhibit a local maximum in both the phenotype and lipid content MDM distributions, with the maximum at an inflammatory phenotype and near the lipid content capacity respectively. These results illustrate how macrophage heterogeneity arises in early atherosclerosis and provide a framework for future model validation through comparison with single-cell RNA sequencing data.

Differential equation models are crucial to scientific processes. The values of model parameters are important for analyzing the behaviour of solutions. A parameter is called globally identifiable if its value can be uniquely determined from the input and output functions. To determine if a parameter estimation problem is well-posed for a given model, one must check if the model parameters are globally identifiable. This problem has been intensively studied for ordinary differential equation models, with theory and several efficient algorithms and software packages developed. A comprehensive theory of algebraic identifiability for PDEs has hitherto not been developed due to the complexity of initial and boundary conditions. Here, we provide theory and algorithms, based on differential algebra, for testing identifiability of polynomial PDE models. We showcase this approach on PDE models arising in the sciences.

Interstitial fluid flow is a feature of many solid tumours. In vitro Experiments have shown that such fluid flow can direct tumour cell movement upstream or downstream depending on the balance between the competing mechanisms of tensotaxis and autologous chemotaxis. In this work we develop a probabilistic-continuum, two-phase model for cell migration in response to interstitial flow. We use a kinetic description for the cell-velocity probability density function, and model the flow-dependent stimuli as forcing terms which bias cell migration upstream and downstream. Using velocity-space averaging, we reformulate the model as a system of continuum equations for the spatio-temporal evolution of the cell volume fraction and flux, in response to forcing terms which depend on the local direction and magnitude of the mechanochemical cues. We specialise our model to describe a one-dimensional cell layer subject to fluid flow. Using a combination of numerical simulations and asymptotic analysis, we delineate the parameter regime where transitions from downstream to upstream cell migration occur. As has been observed experimentally, the model predicts downstream-oriented, chemotactic migration at low cell volume fractions, and upstream-oriented, tensotactic migration at larger volume fractions. We show that the locus of the critical volume fraction, at which the system transitions from downstream to upstream migration, is dominated by the ratio of the rate of chemokine secretion and advection. Our model also predicts that, because the tensotactic stimulus depends strongly on the cell volume fraction, upstream, tensotaxis-dominated migration occurs only transiently when the cells are initially seeded, and transitions to downstream, chemotaxis-dominated migration occur at later times due to the dispersive effect of cell diffusion., Comment: 23 pages, 6 figures. Accepted for publication in Biophysical Journal

Topological data analysis (TDA) is an active field of mathematics for quantifying shape in complex data. Standard methods in TDA such as persistent homology (PH) are typically focused on the analysis of data consisting of a single entity (e.g., cells or molecular species). However, state-of-the-art data collection techniques now generate exquisitely detailed multispecies data, prompting a need for methods that can examine and quantify the relations among them. Such heterogeneous data types arise in many contexts, ranging from biomedical imaging, geospatial analysis, to species ecology. Here, we propose two methods for encoding spatial relations among different data types that are based on Dowker complexes and Witness complexes. We apply the methods to synthetic multispecies data of a tumor microenvironment and analyze topological features that capture relations between different cell types, e.g., blood vessels, macrophages, tumor cells, and necrotic cells. We demonstrate that relational topological features can extract biological insight, including the dominant immune cell phenotype (an important predictor of patient prognosis) and the parameter regimes of a data-generating model. The methods provide a quantitative perspective on the relational analysis of multispecies spatial data, overcome the limits of traditional PH, and are readily computable.

The complex and dynamic crosstalk between tumour and immune cells results in tumours that can exhibit distinct qualitative behaviours - elimination, equilibrium, and escape - and intricate spatial patterns, yet share similar cell configurations in the early stages. We offer a topological approach to analyse time series of spatial data of cell locations (including tumour cells and macrophages) in order to predict malignant behaviour. We propose four topological vectorisations specialised to such cell data: persistence images of Vietoris-Rips and radial filtrations at static time points, and persistence images for zigzag filtrations and persistence vineyards varying in time. To demonstrate the approach, synthetic data are generated from an agent-based model with varying parameters. We compare the performance of topological summaries in predicting - with logistic regression at various time steps - whether tumour niches surrounding blood vessels are present at the end of the simulation, as a proxy for metastasis (i.e., tumour escape). We find that both static and time-dependent methods accurately identify perivascular niche formation, significantly earlier than simpler markers such as the number of tumour cells and the macrophage phenotype ratio. We find additionally that dimension 0 persistence applied to macrophage data, representing multi-scale clusters of the spatial arrangement of macrophages, performs best at this classification task at early time steps, prior to full tumour development, and performs even better when time-dependent data are included; in contrast, topological measures capturing the shape of the tumour, such as tortuosity and punctures in the cell arrangement, perform best at intermediate and later stages. The logistic regression coefficients reveal detailed shape differences between the classes., Comment: 29 pages, 12 figures

Angiogenesis is the process wherein endothelial cells (ECs) form sprouts that elongate from the pre-existing vasculature to create new vascular networks. In addition to its essential role in normal development, angiogenesis plays a vital role in pathologies such as cancer, diabetes and atherosclerosis. Mathematical and computational modelling has contributed to unravelling its complexity. Many existing theoretical models of angiogenic sprouting are based on the 'snail-trail' hypothesis. This framework assumes that leading ECs positioned at sprout tips migrate towards low-oxygen regions while other ECs in the sprout passively follow the leaders' trails and proliferate to maintain sprout integrity. However, experimental results indicate that, contrary to the snail-trail assumption, ECs exchange positions within developing vessels, and the elongation of sprouts is primarily driven by directed migration of ECs. The functional role of cell rearrangements remains unclear. This review of the theoretical modelling of angiogenesis is the first to focus on the phenomenon of cell mixing during early sprouting. We start by describing the biological processes that occur during early angiogenesis, such as phenotype specification, cell rearrangements and cell interactions with the microenvironment. Next, we provide an overview of various theoretical approaches that have been employed to model angiogenesis, with particular emphasis on recent in silico models that account for the phenomenon of cell mixing. Finally, we discuss when cell mixing should be incorporated into theoretical models and what essential modelling components such models should include in order to investigate its functional role., Comment: 27 pages, 9 figures, 1 table

Cancer is a heterogeneous disease and tumours of the same type can differ greatly at the genetic and phenotypic levels. Understanding how these differences impact sensitivity to treatment is an essential step towards patient-specific treatment design. In this paper, we investigate how two different mechanisms for growth control may affect tumour cell responses to fractionated radiotherapy (RT) by extending an existing ordinary differential equation model of tumour growth. In the absence of treatment, this model distinguishes between growth arrest due to nutrient insufficiency and competition for space and exhibits three growth regimes: nutrient-limited (NL), space limited (SL) and bistable (BS), where both mechanisms for growth arrest coexist. We study the effect of RT for tumours in each regime, finding that tumours in the SL regime typically respond best to RT, while tumours in the BS regime typically respond worst to RT. For tumours in each regime, we also identify the biological processes that may explain positive and negative treatment outcomes and the dosing regimen which maximises the reduction in tumour burden., Comment: 33 pages, 22 figures

Organoids are multi-cellular structures which are cultured in vitro from stem cells to resemble specific organs (e.g., brain, liver) in their three-dimensional composition. Dynamic changes in the shape and composition of these model systems can be used to understand the effect of mutations and treatments in health and disease. In this paper, we propose a new technique in the field of topological data analysis for DEtecting Temporal shape changes with the Euler Characteristic Transform (DETECT). DETECT is a rotationally invariant signature of dynamically changing shapes. We demonstrate our method on a data set of segmented videos of mouse small intestine organoid experiments and show that it outperforms classical shape descriptors. We verify our method on a synthetic organoid data set and illustrate how it generalises to 3D. We conclude that DETECT offers rigorous quantification of organoids and opens up computationally scalable methods for distinguishing different growth regimes and assessing treatment effects.

Atherosclerotic plaques form in artery walls due to a chronic inflammatory response driven by lipid accumulation. A key component of the inflammatory response is the interaction between monocyte-derived macrophages and extracellular lipid. Although concentrations of low-density lipoprotein (LDL) and high-density lipoprotein (HDL) particles in the blood are known to affect plaque progression, their impact on the lipid load of plaque macrophages remains unexplored. In this paper, we develop a lipid-structured mathematical model to investigate the impact of blood LDL/HDL levels on plaque composition, and lipid distribution in plaque macrophages. A reduced subsystem, derived by summing the equations of the full model, describes the dynamics of biophysical quantities relating to plaque composition (e.g. total number of macrophages, total amount of intracellular lipid). We also derive a continuum approximation of the model to facilitate analysis of the macrophage lipid distribution. The results, which include time-dependent numerical solutions and asymptotic analysis of the unique steady state solution, indicate that plaque lipid content is sensitive to the influx of LDL relative to HDL capacity. The macrophage lipid distribution evolves in a wave-like manner towards an equilibrium profile which may be monotone decreasing, quasi-uniform or unimodal, attaining its maximum value at a non-zero lipid level. Our model also reveals that macrophage uptake may be severely impaired by lipid accumulation. We conclude that lipid accumulation in plaque macrophages may serve as a partial explanation for the defective uptake of apoptotic cells (efferocytosis) often reported in atherosclerotic plaques., Comment: 44 pages, 13 figures

Atherosclerosis is an inflammatory disease characterised by the formation of plaques, which are deposits of lipids and cholesterol-laden macrophages that form in the artery wall. The inflammation is often non-resolving, due in large part to changes in normal macrophage anti-inflammatory behaviour that are induced by the toxic plaque microenvironment. These changes include higher death rates, defective efferocytic uptake of dead cells, and reduced rates of emigration. We develop a free boundary multiphase model for early atherosclerotic plaques, and we use it to investigate the effects of impaired macrophage anti-inflammatory behaviour on plaque structure and growth. We find that high rates of cell death relative to efferocytic uptake results in a plaque populated mostly by dead cells. We also find that emigration can potentially slow or halt plaque growth by allowing material to exit the plaque, but this is contingent on the availability of live macrophage foam cells in the deep plaque. Finally, we introduce an additional bead species to model macrophage tagging via microspheres, and we use the extended model to explore how high rates of cell death and low rates of efferocytosis and emigration prevent the clearance of macrophages from the plaque.

The Lotka-Volterra model is widely used to model interactions between two species. Here, we generate synthetic data mimicking competitive, mutualistic and antagonistic interactions between two tumor cell lines, and then use the Lotka-Volterra model to infer the interaction type. Structural identifiability of the Lotka-Volterra model is confirmed, and practical identifiability is assessed for three experimental designs: (a) use of a single data set, with a mixture of both cell lines observed over time, (b) a sequential design where growth rates and carrying capacities are estimated using data from experiments in which each cell line is grown in isolation, and then interaction parameters are estimated from an experiment involving a mixture of both cell lines, and (c) a parallel experimental design where all model parameters are fitted to data from two mixtures simultaneously. In addition to assessing each design for practical identifiability, we investigate how the predictive power of the model-i.e., its ability to fit data for initial ratios other than those to which it was calibrated-is affected by the choice of experimental design. The parallel calibration procedure is found to be optimal and is further tested on in silico data generated from a spatially-resolved cellular automaton model, which accounts for oxygen consumption and allows for variation in the intensity level of the interaction between the two cell lines. We use this study to highlight the care that must be taken when interpreting parameter estimates for the spatially-averaged Lotka-Volterra model when it is calibrated against data produced by the spatially-resolved cellular automaton model, since baseline competition for space and resources in the CA model may contribute to a discrepancy between the type of interaction used to generate the CA data and the type of interaction inferred by the LV model., Comment: 25 pages, 18 figures

Collective cell migration plays an essential role in vertebrate development, yet the extent to which dynamically changing microenvironments influence this phenomenon remains unclear. Observations of the distribution of the extracellular matrix (ECM) component fibronectin during the migration of loosely connected neural crest cells (NCCs) lead us to hypothesize that NCC remodeling of an initially punctate ECM creates a scaffold for trailing cells, enabling them to form robust and coherent stream patterns. We evaluate this idea in a theoretical setting by developing an individual-based computational model that incorporates reciprocal interactions between NCCs and their ECM. ECM remodeling, haptotaxis, contact guidance, and cell-cell repulsion are sufficient for cells to establish streams in silico, however additional mechanisms, such as chemotaxis, are required to consistently guide cells along the correct target corridor. Further model investigations imply that contact guidance and differential cell-cell repulsion between leader and follower cells are key contributors to robust collective cell migration by preventing stream breakage. Global sensitivity analysis and simulated gain- and loss-of-function experiments suggest that long-distance migration without jamming is most likely to occur when leading cells specialize in creating ECM fibers, and trailing cells specialize in responding to environmental cues by upregulating mechanisms such as contact guidance., Comment: 66 pages, 19 figures (of which 14 are supplementary)

Tumour spheroids have been the focus of a variety of mathematical models, ranging from Greenspan's classical study of the 1970s through to contemporary agent-based models. Of the many factors that regulate spheroid growth, mechanical effects are perhaps some of the least studied, both theoretically and experimentally, though experimental enquiry has established their significance to tumour growth dynamics. In this tutorial, we formulate a hierarchy of mathematical models of increasing complexity to explore the role of mechanics in spheroid growth, all the while seeking to retain desirable simplicity and analytical tractability. Beginning with the theory of morphoelasticity, which combines solid mechanics and growth, we successively refine our assumptions to develop a somewhat minimal model of mechanically regulated spheroid growth that is free from many unphysical and undesirable behaviours. In doing so, we will see how iterating upon simple models can provide rigorous guarantees of emergent behaviour, which are often precluded by existing, more complex modelling approaches. Perhaps surprisingly, we also demonstrate that the final model considered in this tutorial agrees favourably with classical experimental results, highlighting the potential for simple models to provide mechanistic insight whilst also serving as mathematical examples.

Analysis of single-cell transcriptomics often relies on clustering cells and then performing differential gene expression (DGE) to identify genes that vary between these clusters. These discrete analyses successfully determine cell types and markers; however, continuous variation within and between cell types may not be detected. We propose three topologically motivated mathematical methods for unsupervised feature selection that consider discrete and continuous transcriptional patterns on an equal footing across multiple scales simultaneously. Eigenscores ($\text{eig}_i$) rank signals or genes based on their correspondence to low-frequency intrinsic patterning in the data using the spectral decomposition of the Laplacian graph. The multiscale Laplacian score (MLS) is an unsupervised method for locating relevant scales in data and selecting the genes that are coherently expressed at these respective scales. The persistent Rayleigh quotient (PRQ) takes data equipped with a filtration, allowing the separation of genes with different roles in a bifurcation process (e.g., pseudo-time). We demonstrate the utility of these techniques by applying them to published single-cell transcriptomics data sets. The methods validate previously identified genes and detect additional biologically meaningful genes with coherent expression patterns. By studying the interaction between gene signals and the geometry of the underlying space, the three methods give multidimensional rankings of the genes and visualisation of relationships between them., Comment: 32 pages, 15 figures, 1 table. Revised and published in Entropy, special issue Applications of Topological Data Analysis in the Life Sciences

\emph{In vitro} experiments in which tumour cells are seeded in a gelatinous medium, or hydrogel, show how mechanical interactions between tumour cells and the tissue in which they are embedded, together with local levels of an externally-supplied, diffusible nutrient (e.g., oxygen), affect the tumour's growth dynamics. In this article, we present a mathematical model that describes these \emph{in vitro} experiments. We use the model to understand how tumour growth generates mechanical deformations in the hydrogel and how these deformations in turn influence the tumour's growth. The hydrogel is viewed as a nonlinear hyperelastic material and the tumour is modelled as a two-phase mixture, comprising a viscous tumour cell phase and an isotropic, inviscid interstitial fluid phase. Using a combination of numerical and analytical techniques, we show how the tumour's growth dynamics change as the mechanical properties of the hydrogel vary. When the hydrogel is soft, nutrient availability dominates the dynamics: the tumour evolves to a large equilibrium configuration where the proliferation rate of nutrient-rich cells on the tumour boundary balances the death rate of nutrient-starved cells in the central, necrotic core. As the hydrogel stiffness increases, mechanical resistance to growth increases and the tumour's equilibrium size decreases. Indeed, for small tumours embedded in stiff hydrogels, the inhibitory force experienced by the tumour cells may be so large that the tumour is eliminated. Analysis of the model identifies parameter regimes in which the presence of the hydrogel drives tumour elimination., Comment: 26 pages, 7 figures

We analyse mathematical models in order to understand how microstructural features of vascular networks may affect blood-flow dynamics, and to identify particular characteristics that promote the onset of self-sustained oscillations. By focusing on a simple three-node motif, we predict that network "redundancy", in the form of a redundant vessel connecting two main flow-branches, together with differences in haemodynamic resistance in the branches, can promote the emergence of oscillatory dynamics. We use existing mathematical descriptions for blood rheology and haematocrit splitting at vessel branch-points to construct our flow model; we combine numerical simulations and stability analysis to study the dynamics of the three-node network and its relation to the system's multiple steady-state solutions. While, for the case of equal inlet-pressure conditions, a "trivial" equilibrium solution with no flow in the redundant vessel always exists, we find that it is not stable when other, stable, steady-state attractors exist. In turn, these "nontrivial" steady-state solutions may undergo a Hopf bifurcation into an oscillatory state. We use the branch diameter ratio, together with the inlet haematocrit rate, to construct a two-parameter stability diagram that delineates regimes in which such oscillatory dynamics exist. We show that flow oscillations in this network geometry are only possible when the branch diameters are sufficiently different to allow for a sufficiently large flow in the redundant vessel, which acts as the driving force of the oscillations. These microstructural properties, which were found to promote oscillatory dynamics, could be used to explore sources of flow instability in biological microvascular networks., Comment: 30 pages, 8 figures. To be published in Bulletin of Mathematical Biology. Minor text revision are included in the replacement

Disease complications can alter vascular network morphology and disrupt tissue functioning. Diabetic retinopathy, for example, is a complication of types 1 and 2 diabetes mellitus that can cause blindness. Microvascular diseases are assessed by visual inspection of retinal images, but this can be challenging when diseases exhibit silent symptoms or patients cannot attend in-person meetings. We examine the performance of machine learning algorithms in detecting microvascular disease when trained on statistical and topological summaries of segmented retinal vascular images. We apply our methods to three publicly-available datasets and find that, among the 13 total descriptor vectors we consider, either a statistical Box-counting descriptor vector or a topological Flooding descriptor vector achieves the highest accuracy levels on these datasets. We then created a fourth dataset by merging several datasets: the Box-counting vector outperforms all descriptors on this dataset, including the topological Flooding vector which is sensitive to differences in the annotation styles within the combined dataset. Our work represents a first step to establishing which computational methods are most suitable for identifying microvascular disease as well as some of their current limitations. In the longer term, these methods could be incorporated into automated disease assessment tools.

Longitudinal tumour volume data from head-and-neck cancer patients show that tumours of comparable pre-treatment size and stage may respond very differently to the same radiotherapy fractionation protocol. Mathematical models are often proposed to predict treatment outcome in this context, and have the potential to guide clinical decision-making and inform personalised fractionation protocols. Hindering effective use of models in this context is the sparsity of clinical measurements juxtaposed with the model complexity required to produce the full range of possible patient responses. In this work, we present a compartment model of tumour volume and tumour composition, which, despite relative simplicity, is capable of producing a wide range of patient responses. We then develop novel statistical methodology and leverage a cohort of existing clinical data to produce a predictive model of both tumour volume progression and the associated level of uncertainty that evolves throughout a patient's course of treatment. To capture inter-patient variability, all model parameters are patient specific, with a bootstrap particle filter-like Bayesian approach developed to model a set of training data as prior knowledge. We validate our approach against a subset of unseen data, and demonstrate both the predictive ability of our trained model and its limitations.

The MEK/ERK signalling pathway is involved in cell division, cell specialisation, survival and cell death. Here we study a polynomial dynamical system describing the dynamics of MEK/ERK proposed by Yeung et al. with their experimental setup, data and known biological information. The experimental dataset is a time-course of ERK measurements in different phosphorylation states following activation of either wild-type MEK or MEK mutations associated with cancer or developmental defects. We demonstrate how methods from computational algebraic geometry, differential algebra, Bayesian statistics and computational algebraic topology can inform the model reduction, identification and parameter inference of MEK variants, respectively. Throughout, we show how this algebraic viewpoint offers a rigorous and systematic analysis of such models.

We present a mathematical model that describes how tumour heterogeneity evolves in a tissue slice that is oxygenated by a single blood vessel. Phenotype is identified with the stemness level of a cell, $s$, that determines its proliferative capacity, apoptosis propensity and response to treatment. Our study is based on numerical bifurcation analysis and dynamical simulations of a system of coupled non-local (in phenotypic space) partial differential equations that links the phenotypic evolution of the tumour cells to local oxygen levels in the tissue. In our formulation, we consider a 1D geometry where oxygen is supplied by a blood vessel located on the domain boundary and consumed by the tumour cells as it diffuses through the tissue. For biologically relevant parameter values, the system exhibits multiple steady states; in particular, depending on the initial conditions, the tumour is either eliminated ("tumour-extinction") or it persists ("tumour-invasion"). We conclude by using the model to investigate tumour responses to radiotherapy (RT), and focus on establishing which RT strategies can eliminate the tumour. Numerical simulations reveal how phenotypic heterogeneity evolves during treatment and highlight the critical role of tissue oxygen levels on the efficacy of radiation protocols that are commonly used clinically.

In this paper, we carry out a travelling-wave analysis of a model of tumour invasion with degenerate, cross-dependent diffusion. We consider two types of invasive fronts of tumour tissue into extracellular matrix (ECM), which represents healthy tissue. These types differ according to whether the density of ECM far ahead of the wave front is maximal or not. In the former case, we use a shooting argument to prove that there exists a unique travelling wave solution for any positive propagation speed. In the latter case, we further develop this argument to prove that there exists a unique travelling wave solution for any propagation speed greater than or equal to a strictly positive minimal wave speed. Using a combination of analytical and numerical results, we conjecture that the minimal wave speed depends monotonically on the degradation rate of ECM by tumour cells and the ECM density far ahead of the front., Comment: 41 pages (21 for the main paper, 17 for the supplementary material and 3 for the bibliography) with 6 figures (4 in the main paper and 2 in the supplementary material)

Hybrid multiscale modelling has emerged as a useful framework for modelling complex biological phenomena. However, when accounting for stochasticity in the internal dynamics of agents, these models frequently become computationally expensive. Traditional techniques to reduce the computational intensity of such models can lead to a reduction in the richness of the dynamics observed, compared to the original system. Here we use large deviation theory to decrease the computational cost of a spatially-extended multi-agent stochastic system with a region of multi-stability by coarse-graining it to a continuous time Markov chain on the state space of stable steady states of the original system. Our technique preserves the original description of the stable steady states of the system and accounts for noise-induced transitions between them. We apply the method to a bistable system modelling phenotype specification of cells driven by a lateral inhibition mechanism. For this system, we demonstrate how the method may be used to explore different pattern configurations and unveil robust patterns emerging on longer timescales. We then compare the full stochastic, coarse-grained and mean-field descriptions via pattern quantification metrics and in terms of the numerical cost of each method. Our results show that the coarse-grained system exhibits the lowest computational cost while preserving the rich dynamics of the stochastic system. The method has the potential to reduce the computational complexity of hybrid multiscale models, making them more tractable for analysis, simulation and hypothesis testing.

In this work, we use a combination of formal upscaling and data-driven machine learning for explicitly closing a nonlinear transport and reaction process in a multiscale tissue. The classical effectiveness factor model is used to formulate the macroscale reaction kinetics. We train a multilayer perceptron network using training data generated by direct numerical simulations over microscale examples. Once trained, the network is used for numerically solving the upscaled (coarse-grained) differential equation describing mass transport and reaction in two example tissues. The network is described as being explicit in the sense that the network is trained using macroscale concentrations and gradients of concentration as components of the feature space. Network training and solutions to the macroscale transport equations were computed for two different tissues. The two tissue types (brain and liver) exhibit markedly different geometrical complexity and spatial scale (cell size and sample size). The upscaled solutions for the average concentration are compared with numerical solutions derived from the microscale concentration fields by a posteriori averaging. There are two outcomes of this work of particular note: 1) we find that the trained network exhibits good generalizability, and it is able to predict the effectiveness factor with high fidelity for realistically-structured tissues despite the significantly different scale and geometry of the two example tissue types; and 2) the approach results in an upscaled PDE with an effectiveness factor that is predicted (implicitly) via the trained neural network. This latter result emphasizes our purposeful connection between conventional averaging methods with the use of machine learning for closure; this contrasts with some machine learning methods for upscaling where the exact form of the macroscale equation remains unknown.

In cancer, treatment failure and disease recurrence have been associated with small subpopulations of cancer cells with a stem-like phenotype. In this paper, we develop and investigate a phenotype-structured model of solid tumour growth in which cells are structured by a stemness level, which varies continuously between stem-like and terminally differentiated behaviours. Cell evolution is driven by proliferation and apoptosis, as well as advection and diffusion with respect to the stemness structure variable. We use the model to investigate how the environment, in particular oxygen levels, affects the tumour's population dynamics and composition, and its response to radiotherapy. We use a combination of numerical and analytical techniques to quantify how under physiological oxygen levels the cells evolve to a differentiated phenotype and under low oxygen level (i.e., hypoxia) they de-differentiate. Under normoxia, the proportion of cancer stem cells is typically negligible and the tumour may ultimately become extinct whereas under hypoxia cancer stem cells comprise a dominant proportion of the tumour volume, enhancing radio-resistance and favouring the tumour's long-term survival. We then investigate how such phenotypic heterogeneity impacts the tumour's response to treatment with radiotherapy under normoxia and hypoxia. Of particular interest is establishing how the presence of radio-resistant cancer stem cells can facilitate a tumour's regrowth following radiotherapy. We also use the model to show how radiation-induced changes in tumour oxygen levels can give rise to complex re-growth dynamics. For example, transient periods of hypoxia induced by damage to tumour blood vessels may rescue the cancer cell population from extinction and drive secondary regrowth. Further model extensions to account for spatial variation are also discussed briefly.

Angiogenesis is the process by which blood vessels form from pre-existing vessels. It plays a key role in many biological processes, including embryonic development and wound healing, and contributes to many diseases including cancer and rheumatoid arthritis. The structure of the resulting vessel networks determines their ability to deliver nutrients and remove waste products from biological tissues. Here we simulate the Anderson-Chaplain model of angiogenesis at different parameter values and quantify the vessel architectures of the resulting synthetic data. Specifically, we propose a topological data analysis (TDA) pipeline for systematic analysis of the model. TDA is a vibrant and relatively new field of computational mathematics for studying the shape of data. We compute topological and standard descriptors of model simulations generated by different parameter values. We show that TDA of model simulation data stratifies parameter space into regions with similar vessel morphology. The methodologies proposed here are widely applicable to other synthetic and experimental data including wound healing, development, and plant biology.

Advances in imaging techniques enable high resolution 3D visualisation of vascular networks over time and reveal abnormal structural features such as twists and loops, and their quantification is an active area of research. Here we showcase how topological data analysis (TDA), the mathematical field that studies `shape' of data, can characterise the geometric, spatial and temporal organisation of vascular networks. We propose two topological lenses to study vasculature, which capture inherent multi-scale features and vessel connectivity, and surpass the single scale analysis of existing methods. We analyse images collected using intravital and ultramicroscopy modalities and quantify spatio-temporal variation of twists, loops, and avascular regions (voids) in 3D vascular networks. This topological approach validates and quantifies known qualitative trends such as dynamic changes in tortuosity and loops in response to antibodies that modulate vessel sprouting; furthermore, it quantifies the effect of radiotherapy on vessel architecture.

We present a mechanical model of tissue homeostasis that is specialised to the intestinal crypt. Growth and deformation of the crypt, idealised as a line of cells on a substrate, are modelled using morphoelastic rod theory. Alternating between Lagrangian and Eulerian mechanical descriptions enables us precisely to characterise the dynamic nature of tissue homeostasis, whereby the proliferative structure and morphology are static in the Eulerian frame, but there is active migration of Lagrangian material points out of the crypt. Assuming mechanochemical growth, we identify the necessary conditions for homeostasis, reducing the full, time-dependent system to a static boundary value problem characterising a spatially heterogeneous "treadmilling" state. We extract essential features of crypt homeostasis, such as the morphology, the proliferative structure, the migration velocity, and the sloughing rate. We also derive closed-form solutions for growth and sloughing dynamics in homeostasis, and show that mechanochemical growth is sufficient to generate the observed proliferative structure of the crypt. Key to this is the concept of threshold-dependent mechanical feedback, that regulates an established Wnt signal for biochemical growth. Numerical solutions demonstrate the importance of crypt morphology on homeostatic growth, migration, and sloughing, and highlight the value of this framework as a foundation for studying the role of mechanics in homeostasis., Comment: 21 pages, 7 figures

Deterministic continuum models formulated in terms of non-local partial differential equations for the evolutionary dynamics of populations structured by phenotypic traits have been used recently to address open questions concerning the adaptation of asexual species to periodically fluctuating environmental conditions. These deterministic continuum models are usually defined on the basis of population-scale phenomenological assumptions and cannot capture adaptive phenomena that are driven by stochastic variability in the evolutionary paths of single individuals. In this paper, we develop a stochastic individual-based model for the coevolution between two competing phenotype-structured cell populations that are exposed to time-varying nutrient levels and undergo spontaneous, heritable phenotypic variations with different probabilities. The evolution of every cell is described by a set of rules that result in a discrete-time branching random walk on the space of phenotypic states. We formally show that the deterministic continuum counterpart of this model comprises a system of non-local partial differential equations for the cell population density functions coupled with an ordinary differential equation for the nutrient concentration. We compare the individual-based model and its continuum analogue, focussing on scenarios whereby the predictions of the two models differ. Our results clarify the conditions under which significant differences between the two models can emerge due to stochastic effects associated with small population levels. These differences arise in the presence of low probabilities of phenotypic variation, and become more apparent when the two populations are characterised by less fit initial mean phenotypes and smaller initial levels of phenotypic heterogeneity., Comment: 15 pages, 9 figures

We study a five-compartment mathematical model originally proposed by Kuznetsov et al. (1994) to investigate the effect of nonlinear interactions between tumour and immune cells in the tumour microenvironment, whereby immune cells may induce tumour cell death, and tumour cells may inactivate immune cells. Exploiting a separation of timescales in the model, we use the method of matched asymptotics to derive a new two-dimensional, long-timescale, approximation of the full model, which differs from the quasi-steady-state approximation introduced by Kuznetsov et al. (1994), but is validated against numerical solutions of the full model. Through a phase-plane analysis, we show that our reduced model is excitable, a feature not traditionally associated with tumour-immune dynamics. Through a systematic parameter sensitivity analysis, we demonstrate that excitability generates complex bifurcating dynamics in the model. These are consistent with a variety of clinically observed phenomena, and suggest that excitability may underpin tumour-immune interactions. The model exhibits the three stages of immunoediting - elimination, equilibrium, and escape, via stable steady states with different tumour cell concentrations. Such heterogeneity in tumour cell numbers can stem from variability in initial conditions and/or model parameters that control the properties of the immune system and its response to the tumour. We identify different biophysical parameter targets that could be manipulated with immunotherapy in order to control tumour size, and we find that preferred strategies may differ between patients depending on the strength of their immune systems, as determined by patient-specific values of associated model parameters., Comment: 29 pages, 15 figures

Atherosclerosis is characterised by the growth of fatty plaques in the inner (intimal) layer of the artery wall. In mature plaques, vascular smooth muscle cells (SMCs) are recruited from the adjacent medial layer to deposit a cap of fibrous collagen over the fatty plaque core. The fibrous cap isolates the thrombogenic content of the plaque from the bloodstream and prevents the formation of blood clots that cause myocardial infarction or stroke. Despite the important protective role of the cap, the mechanisms that regulate cap formation and maintenance are not well understood. It remains unclear why certain caps become stable, while others become vulnerable to rupture. We develop a multiphase PDE model with non-standard boundary conditions to investigate collagen cap formation by SMCs in response to growth factor signals from the endothelium. Diffusible platelet-derived growth factor (PDGF) stimulates SMC migration, proliferation and collagen degradation, while diffusible transforming growth factor (TGF)-$\beta$ stimulates SMC collagen synthesis and inhibits collagen degradation. The model SMCs respond haptotactically to gradients in the collagen phase and have reduced rates of migration and proliferation in dense collagenous tissue. The model, which is parameterised using a range of in vivo and in vitro experimental data, reproduces several observations from studies of plaque growth in atherosclerosis-prone mice. Numerical simulations and model analysis demonstrate that a stable cap can be formed by a relatively small SMC population and emphasise the critical role of TGF-$\beta$ in effective cap formation and maintenance. These findings provide unique insight into the cellular and biochemical mechanisms that may lead to plaque destabilisation and rupture. This work represents an important step towards the development of a comprehensive in silico plaque.

Understanding how the spatial structure of blood vessel networks relates to their function in healthy and abnormal biological tissues could improve diagnosis and treatment for diseases such as cancer. New imaging techniques can generate multiple, high-resolution images of the same tissue region, and show how vessel networks evolve during disease onset and treatment. Such experimental advances have created an exciting opportunity for discovering new links between vessel structure and disease through the development of mathematical tools that can analyse these rich datasets. Here we explain how topological data analysis (TDA) can be used to study vessel network structures. TDA is a growing field in the mathematical and computational sciences, that consists of algorithmic methods for identifying global and multi-scale structures in high-dimensional data sets that may be noisy and incomplete. TDA has identified the effect of ageing on vessel networks in the brain and more recently proposed to study blood flow and stenosis. Here we present preliminary work which shows how TDA of spatial network structure can be used to characterise tumour vasculature.

Living species, ranging from bacteria to animals, exist in environmental conditions that exhibit spatial and temporal heterogeneity which requires them to adapt. Risk-spreading through spontaneous phenotypic variations is a known concept in ecology, which is used to explain how species may survive when faced with the evolutionary risks associated with temporally varying environments. In order to support a deeper understanding of the adaptive role of spontaneous phenotypic variations in fluctuating environments, we consider a system of non-local partial differential equations modelling the evolutionary dynamics of two competing phenotype-structured populations in the presence of periodically oscillating nutrient levels. The two populations undergo spontaneous phenotypic variations at different rates. The phenotypic state of each individual is represented by a continuous variable, and the phenotypic landscape of the populations evolves in time due to variations in the nutrient level. Exploiting the analytical tractability of our model, we study the long-time behaviour of the solutions to obtain a detailed mathematical depiction of evolutionary dynamics. The results suggest that when nutrient levels undergo small and slow oscillations, it is evolutionarily more convenient to rarely undergo spontaneous phenotypic variations. Conversely, under relatively large and fast periodic oscillations in the nutrient levels, which bring about alternating cycles of starvation and nutrient abundance, higher rates of spontaneous phenotypic variations confer a competitive advantage. We discuss the implications of our results in the context of cancer metabolism., Comment: 33 pages, 8 figures

Motivation: The scratch assay is a standard experimental protocol used to characterize cell migration. It can be used to identify genes that regulate migration and evaluate the efficacy of potential drugs that inhibit cancer invasion. In these experiments, a scratch is made on a cell monolayer and recolonisation of the scratched region is imaged to quantify cell migration rates. A drawback of this methodology is the lack of its reproducibility resulting in irregular cell-free areas with crooked leading edges. Existing quantification methods deal poorly with such resulting irregularities present in the data. Results: We introduce a new quantification method that can analyse low quality experimental data. By considering in-silico and in-vitro data, we show that the method provides a more accurate statistical classification of the migration rates than two established quantification methods. The application of this method will enable the quantification of migration rates of scratch assay data previously unsuitable for analysis. Availability and Implementation: The source code and the implementation of the algorithm as a GUI along with an example dataset and user instructions, are available in https://bitbucket.org/anavictoria-ponce/local_migration_quantification_scratch_assays/src/master/. The datasets are available in https://ganymed.math.uni-heidelberg.de/~victoria/publications.shtml.

Atherosclerotic plaque growth is characterised by chronic inflammation that promotes accumulation of cellular debris and extracellular fat in the inner artery wall. This material is highly thrombogenic, and plaque rupture can lead to the formation of blood clots that occlude major arteries and cause myocardial infarction or stroke. In advanced plaques, vascular smooth muscle cells (SMCs) migrate from deeper in the artery wall to synthesise a cap of fibrous tissue that stabilises the plaque and sequesters the thrombogenic plaque content from the bloodstream. The fibrous cap provides crucial protection against the clinical consequences of atherosclerosis, but the mechanisms of cap formation are poorly understood. In particular, it is unclear why certain plaques become stable and robust while others become fragile and vulnerable to rupture. We develop a multiphase model with non-standard boundary conditions to investigate early fibrous cap formation in the atherosclerotic plaque. The model is parameterised using a range of in vitro and in vivo data, and includes highly nonlinear mechanisms of SMC proliferation and migration in response to an endothelium-derived chemical signal. We demonstrate that the model SMC population naturally evolves towards a steady-state, and predict a rate of cap formation and a final plaque SMC content consistent with experimental observations in mice. Parameter sensitivity simulations show that SMC proliferation makes a limited contribution to cap formation, and highlight that stable cap formation relies on a critical balance between SMC recruitment to the plaque, SMC migration within the plaque and SMC loss by apoptosis. The model represents the first detailed in silico study of fibrous cap formation in atherosclerosis, and establishes a multiphase modelling framework that can be readily extended to investigate many other aspects of plaque development.