Research overview

I am an applied mathematician working in Mathematical Biology and I have been particularly interested in the formulation and study of mathematical models of the spatio-temporal and evolutionary dynamics of cell populations in cancer and development. I have been focussing on continuous, deterministic models of tissue-level cell population dynamics, which translate mathematically into systems of nonlinear, and often nonlocal, partial differential equations (PDEs).
The nonlinear, nonlocal and at times stiff nature of these PDEs poses a series of interesting analytical and numerical challenges, which can be tackled by means of formal asymptotic methods (e.g. Hamilton-Jacobi formalism), linear stability analyses and appropriate numerical schemes preventing the emergence of spurious oscillations (e.g. mixed finite difference and finite volume schemes, flux limiting schemes). Moreover, model integration with experimental data poses a series of statistical challenges, that may be tackled with likelihood-maximisation and bootstrapping algorithms.
These models are capable of shedding light on the hidden mechanisms responsible for the spatial sorting of cell populations at the tissue scale, whether this may arise due to the adaptive dynamics of cancer cells and their nonlinear interaction with nutrients entering from the blood vessels, or from more complex forms of cell movement mediated by intracellular signalling and their dynamic interaction with the extracellular matrix. Ultimately, these models may complement experimental research and be of help for better treatment design.
Scroll down for a list of my research outputs

Preprints

  • 12. T. Lorenzi, K.J. Painter, C. Villa, Phenotype structuring in collective cell migration: a tutorial of mathematical models and methods, 2024
    Preprint: arXiv:2410.13629
    Abstract
    Populations are heterogeneous, deviating in numerous ways. Phenotypic diversity refers to the range of traits or characteristics across a population, where for cells this could be the levels of signalling, movement and growth activity, etc. Clearly, the phenotypic distribution - and how this changes over time and space - could be a major determinant of population-level dynamics. For instance, across a cancerous population, variations in movement, growth, and ability to evade death may determine its growth trajectory and response to therapy. In this review, we discuss how classical partial differential equation (PDE) approaches for modelling cellular systems and collective cell migration can be extended to include phenotypic structuring. The resulting non-local models - which we refer to as phenotype-structured partial integro-differential equations (PS-PIDEs) - form a sophisticated class of models with rich dynamics. We set the scene through a brief history of structured population modelling, and then review the extension of several classic movement models - including the Fisher-KPP and Keller-Segel equations - into a PS-PIDE form. We proceed with a tutorial-style section on derivation, analysis, and simulation techniques. First, we show a method to formally derive these models from underlying agent-based models. Second, we recount travelling waves in PDE models of spatial spread dynamics and concentration phenomena in non-local PDE models of evolutionary dynamics, and combine the two to deduce phenotypic structuring across travelling waves in PS-PIDE models. Third, we discuss numerical methods to simulate PS-PIDEs, illustrating with a simple scheme based on the method of lines and noting the finer points of consideration. We conclude with a discussion of future modelling and mathematical challenges.
  • 11. A.P. Browning, R. Crossley, C. Villa, P. K. Maini, A.L. Jenner, T. Cassidy and S. Hamis, Identifiability of heterogeneous phenotype adaptation from low-cell-count experiments and a stochastic model, 2024
    Abstract
    Adaptive resistance contributes significantly to treatment failure in many cancers. Despite the increased prevalence of experimental studies that interrogate this phenomenon, there remains a lack of applicable quantitative tools to characterise data, and importantly to distinguish between resistance as a discrete phenotype and a (potentially heterogeneous) continuous distribution of phenotypes. To address this, we develop a stochastic individual-based model of adaptive resistance in low-cell-count proliferation assays. That our model corresponds probabilistically to common partial differential equation models of resistance allows us to formulate a likelihood that captures the intrinsic noise ubiquitous to such experiments. We apply our framework to assess the identifiability of key model parameters in several population-level data collection regimes; in particular, parameters relating to the adaptation velocity and within-population heterogeneity. Significantly, we find that heterogeneity is practically non-identifiable from both cell count and proliferation marker data, implying that population-level behaviours may be well characterised by homogeneous ordinary differential equation models. Additionally, we demonstrate that population-level data are insufficient to distinguish resistance as a discrete phenotype from a continuous distribution of phenotypes. Our results inform the design of both future experiments and future quantitative analyses that probe adaptive resistance in cancer.
  • 10. S. Hamis, A.P. Browning, A.L. Jenner, C. Villa, P. K. Maini and T. Cassidy, Growth rate-driven modelling reveals how phenotypic adaptation drives drug resistance in BRAFV600E-mutant melanoma, 2024
    Abstract
    Phenotypic adaptation, the ability of cells to change phenotype in response to external pressures, has been identified as a driver of drug resistance in cancer. To quantify phenotypic adaptation in BRAFV600E-mutant melanoma, we develop a theoretical model that emerges from data analysis of WM239A-BRAFV600E growth rates in response to drug challenge with the BRAF-inhibitor encorafenib. Our model constitutes a cell population model in which each cell is individually described by one of multiple discrete and plastic phenotype states that are directly linked to drug-dependent net growth rates and, by extension, drug resistance. Data-matched simulations reveal that phenotypic adaptation in the cells is directed towards states of high net growth rates, which enables evasion of drug-effects. The model subsequently provides an explanation for when and why intermittent treatments outperform continuous treatments in vitro, and demonstrates the benefits of not only targeting, but also leveraging, phenotypic adaptation in treatment protocols.
  • 9. C. Villa, P. K. Maini, A.P. Browning, A.L. Jenner, S. Hamis and T. Cassidy, Reducing phenotype-structured PDE models of cancer evolution to systems of ODEs: a generalised moment dynamics approach, 2024
    Abstract
    Intratumour phenotypic heterogeneity is nowadays understood to play a critical role in disease progression and treatment failure. Accordingly, there has been increasing interest in the development of mathematical models capable of capturing its role in cancer cell adaptation. This can be systematically achieved by means of models comprising phenotype-structured nonlocal partial differential equations, tracking the evolution of the phenotypic density distribution of the cell population, which may be compared to gene and protein expression distributions obtained experimentally. Nevertheless, given the high analytical and computational cost of solving these models, much is to be gained from reducing them to systems of ordinary differential equations for the moments of the distribution. We propose a generalised method of model-reduction, relying on the use of a moment generating function, Taylor series expansion and truncation closure, to reduce a nonlocal reaction-advection-diffusion equation, with general phenotypic drift and proliferation rate functions, to a system of moment equations up to arbitrary order. Our method extends previous results in the literature, which we address via two examples, by removing any a priori assumption on the shape of the distribution, and provides a flexible framework for mathematical modellers to account for the role of phenotypic heterogeneity in cancer adaptive dynamics, in a simpler mathematical framework.
  • 8. L. Almeida, A. Poulain, A. Pourtier, C. Villa, Mathematical modelling of the contribution of senescent fibroblasts to basement membrane digestion during carcinoma invasion, 2024
    Preprint: hal-04574340
    Abstract
    Senescent cells have been recognized to play major roles in tumor progression and are nowadays included in the hallmarks of cancer. Our work aims to develop a mathematical model capable of capturing a pro-invasion effect of senescent fibroblasts located in the conjunctive tissue. We focus in the present article on the first moments of the invasion cascade. Considering a localized epithelial tumor, we model the digestion of the collagen fibers of the basement membrane by the proteolytic enzyme MMP-2. The activation of MMP-2 is modelled in detail, as MT1-MMPs bound to the surface of tumor cells interact with proMMPs and TIMPs, proteins enriched in the secretome of senescent Cancer-Associated Fibroblasts, along with its inhibition by TIMPs. Using numerical simulations of the model, calibrated via an extensive literature search, reproducing biologically relevant scenarios, we test the model’s suitability to investigate the effect on basement membrane digestion of fibroblasts presenting a senescence-associated secretory phenotype. Via model reduction, steady state and global sensitivity analyses, we identify the most influential parameters in view of their calibration with empirical data. We conclude the paper discussing mathematical and interdisciplinary perspectives.

Journal articles

  • 7. F. Padovano, C. Villa, The development of drug resistance in metastatic tumours under chemotherapy: an evolutionary perspective, Journal of Theoretical Biology, 595(1):111957, 2024
    Abstract
    We present a mathematical model of the evolutionary dynamics of a metastatic tumour under chemotherapy, comprising non-local partial differential equations for the phenotype-structured cell populations in the primary tumour and its metastasis. These equations are coupled with a physiologically-based pharmacokinetic model of drug delivery, implementing a realistic delivery schedule. The model is carefully calibrated from the literature, focusing on BRAF-mutated melanoma treated with Dabrafenib as a case study. By means of long-time asymptotic analysis, global sensitivity analysis and numerical simulations, we explore the impact of cell migration from the primary to the metastatic site, physiological aspects of the tumour sites and drug dose on the development of drug resistance and treatment efficacy. Our findings provide a possible explanation for empirical evidence indicating that chemotherapy may foster metastatic spread and that metastatic sites may be less impacted by chemotherapy.
  • 6. L. Almeida, J.A. Denis, N. Ferrand, T. Lorenzi, A. Prunet, M. Sabbah, C. Villa, Evolutionary dynamics of glucose-deprived cancer cells: insights from experimentally-informed mathematical modelling, Journal of the Royal Society Interface, 21:20230587, 2024
    Abstract
    Glucose is a primary energy source for cancer cells. Several lines of evidence support the idea that monocarboxylate transporters, such as MCT1, elicit metabolic reprogramming of cancer cells in glucose-poor environments, allowing them to re-use lactate, a by-product of glucose metabolism, as an alternative energy source with serious consequences for disease progression. We employ a synergistic experimental and mathematical modelling approach to explore the evolutionary processes at the root of cancer cell adaptation to glucose deprivation, with particular focus on the mechanisms underlying the increase in MCT1 expression observed in glucose-deprived aggressive cancer cells. Data from in vitro experiments on breast cancer cells are used to inform and calibrate a mathematical model that comprises a partial integro-differential equation for the dynamics of a population of cancer cells structured by the level of MCT1 expression. Analytical and numerical results of this model suggest that environment-induced changes in MCT1 expression mediated by lactate-associated signalling pathways enable a prompt adaptive response of glucose-deprived cancer cells, while fluctuations in MCT1 expression due to epigenetic changes create the substrate for environmental selection to act upon, speeding up the selective sweep underlying cancer cell adaptation to glucose deprivation, and may constitute a long-term bet-hedging mechanism.
  • 5. C. Villa, A. Gerisch, M.A.J. Chaplain, A novel nonlocal partial differential equation model of endothelial progenitor cell cluster formation during the early stages of vasculogenesis, Journal of Theoretical Biology, 534(1):110963, 2022
    Abstract
    Neovascularisation is essential for tissue development and regeneration, in addition to playing a key role in pathological settings such as ischemia and tumour development. Experimental findings in the past two decades have led to the identification of a new mechanism of neovascularisation, cluster-based vasculogenesis, during which endothelial progenitor cells (EPCs) mobilised from the bone marrow are capable of bridging distant vascular beds in a variety of hypoxic settings in vivo. This process is characterised by the formation of EPC clusters during its early stages and, while much progress has been made in identifying various mechanisms underlying cluster formation, we are still far from a comprehensive description of such spatio-temporal dynamics. In order to achieve this, we propose a novel mathematical model of the early stages of cluster-based vasculogenesis, comprising of a system of nonlocal partial differential equations including key mechanisms such as endogenous chemotaxis, matrix degradation, cell proliferation and cell-to-cell adhesion. We conduct a linear stability analysis on the system, solve the equations numerically, conduct a parametric analysis of the numerical solutions of the 1D problem to investigate the role of underlying dynamics on the speed of cluster formation and the size of clusters, and verify the key results of the parametric analysis with simulations of the 2D problem. Our results, which qualitatively compare with data from in vitro experiments, elucidate the complementary role played by endogenous chemotaxis and matrix degradation in the formation of clusters, and they indicate that previous approaches to the nonlocal modelling of cell-to-cell adhesion, while they capture the aggregating effect of cell-to-cell adhesion, are not sufficient to capture its stabilising effect on clusters, and new continuum cell-adhesion modelling strategies are required.
  • 4. F. Mottes, C. Villa, M. Osella, M. Caselle, The impact of whole genome duplications on the human gene regulatory networks, PLOS Computational Biology, 17(12):e1009638, 2021
    Abstract
    This work studies the effects of the two rounds of Whole Genome Duplication (WGD) at the origin of the vertebrate lineage on the architecture of the human gene regulatory networks. We integrate information on transcriptional regulation, miRNA regulation, and protein-protein interactions to comparatively analyse the role of WGD and Small Scale Duplications (SSD) in the structural properties of the resulting multilayer network. We show that complex network motifs, such as combinations of feed-forward loops and bifan arrays, deriving from WGD events are specifically enriched in the network. Pairs of WGD-derived proteins display a strong tendency to interact both with each other and with common partners and WGD-derived transcription factors play a prominent role in the retention of a strong regulatory redundancy. Combinatorial regulation and synergy between different regulatory layers are in general enhanced by duplication events, but the two types of duplications contribute in different ways. Overall, our findings suggest that the two WGD events played a substantial role in increasing the multi-layer complexity of the vertebrate regulatory network by enhancing its combinatorial organization, with potential consequences on its overall robustness and ability to perform high-level functions like signal integration and noise control.
  • 3. C. Villa, M.A.J. Chaplain, A. Gerisch, T. Lorenzi, Mechanical models of pattern and form in biological tissues: the role of stress-strain constitutive equations, Bulletin of Mathematical Biology, 83(4):1-38, 2021
    Abstract
    Mechanochemical models of pattern formation in biological tissues have been used to study a variety of biomedical systems and describe the physical interactions between cells and their local surroundings. These models generally consist of a balance equation for the cell density, one for the density of the extracellular matrix (ECM), and a force-balance equation describing the mechanical equilibrium of the cell-ECM system. Assuming this system can be regarded as an isotropic linear viscoelastic material, the force-balance equation is often defined using the Kelvin-Voigt model of linear viscoelasticity to represent the stress-strain relation of the ECM. However, due to the multifaceted bio-physical nature of the ECM constituents, there are rheological aspects that cannot be effectively captured by this model and, therefore, depending on the type of biological tissue considered, other constitutive models of linear viscoelasticity may be better suited. In this work, we systematically assess the pattern formation potential of different stress-strain constitutive equations for the ECM within a mechanical model of pattern formation in biological tissues. The results obtained through linear stability analysis support the idea that constitutive equations capturing viscous flow and permanent set (Maxwell model, Jeffrey model) have a pattern formation potential much higher than the others (Kelvin-Voigt model, standard linear solid model), further confirmed by the results of our numerical simulations. Our findings suggest that further empirical work is required to acquire detailed quantitative information on the mechanical properties of components of the ECM in different biological tissues in order to furnish mechanochemical models of pattern formation with stress-strain constitutive equations for the ECM that provide a more faithful representation of the underlying tissue rheology.
  • 2. C. Villa, M.A.J. Chaplain, T. Lorenzi, Evolutionary dynamics in vascularised tumours under chemotherapy: Mathematical modelling, asymptotic analysis and numerical simulations, Vietnam Journal of Mathematics, 49(1):143–167, 2021
    Abstract
    We consider a mathematical model for the evolutionary dynamics of tumour cells in vascularised tumours under chemotherapy. The model comprises a system of coupled partial integro-differential equations for the phenotypic distribution of tumour cells, the concentration of oxygen and the concentration of a chemotherapeutic agent. In order to disentangle the impact of different evolutionary parameters on the emergence of intra-tumour phenotypic heterogeneity and the development of resistance to chemotherapy, we construct explicit solutions to the equation for the phenotypic distribution of tumour cells and provide a detailed quantitative characterisation of the long-time asymptotic behaviour of such solutions. Analytical results are integrated with numerical simulations of a calibrated version of the model based on biologically consistent parameter values. The results obtained provide a theoretical explanation for the observation that the phenotypic properties of tumour cells in vascularised tumours vary with the distance from the blood vessels. Moreover, we demonstrate that lower oxygen levels may correlate with higher levels of phenotypic variability, which suggests that the presence of hypoxic regions supports intra-tumour phenotypic heterogeneity. Finally, the results of our analysis put on a rigorous mathematical basis the idea, previously suggested by formal asymptotic results and numerical simulations, that hypoxia favours the selection for chemoresistant phenotypic variants prior to treatment. Consequently, this facilitates the development of resistance following chemotherapy.
  • 1. C. Villa, M.A.J. Chaplain, T. Lorenzi, Modelling phenotypic heterogeneity in vascularised tumours, SIAM Journal on Applied Mathematics, 81(2):434–453, 2021
    Abstract
    We present a mathematical study of the emergence of phenotypic heterogeneity in vascularised tumours. Our study is based on formal asymptotic analysis and numerical simulations of a system of non-local parabolic equations that describes the phenotypic evolution of tumour cells and their nonlinear dynamic interactions with the oxygen, which is released from the intratumoural vascular network. Numerical simulations are carried out both in the case of arbitrary distributions of intratumour blood vessels and in the case where the intratumoural vascular network is reconstructed from clinical images obtained using dynamic optical coherence tomography. The results obtained support a more in-depth theoretical understanding of the eco-evolutionary process which underpins the emergence of phenotypic heterogeneity in vascularised tumours. In particular, our results offer a theoretical basis for empirical evidence indicating that the phenotypic properties of cancer cells in vascularised tumours vary with the distance from the blood vessels, and establish a relation between the degree of tumour tissue vascularisation and the level of intratumour phenotypic heterogeneity.

Conference proceedings

  • P1. T. Lorenzi, F.R. Macfarlane, C. Villa, Discrete and continuum models for the evolutionary and spatial dynamics of cancer: a very short introduction through two case studies, (pp. 359-380) in Trends in Biomathematics: Modeling Cells, Flows, Epidemics, and the Environment, Ed. R. Mondaini, Springer, Cham, 2019
    Abstract
    We give a very short introduction to discrete and continuum models for the evolutionary and spatial dynamics of cancer through two case studies: a model for the evolutionary dynamics of cancer cells under cytotoxic therapy and a model for the mechanical interaction between healthy and cancer cells during tumour growth. First we develop the discrete models, whereby the dynamics of single cells are described through a set of rules that result in branching random walks. Then we present the corresponding continuum models, which are formulated in terms of non-local and nonlinear partial differential equations, and we summarise the key properties of their solutions. Finally, we carry out numerical simulations of the discrete models and we construct numerical solutions of the corresponding continuum models. The biological implications of the results obtained are briefly discussed.

Theses

  • T1. C. Villa, Partial differential equation modelling in cancer and development, PhD thesis, University of St Andrews, St Andrews, 2022
    Abstract
    This thesis explores various partial differential equation (PDE) models of the spatiotemporal and evolutionary dynamics of cell populations in different problems in cancer and development. In particular, these models are used to investigate: (i) the emergence of intratumour phenotypic heterogeneity and the development of chemotherapeutic resistance in vascularised tumours; (ii) the formation of endothelial progenitor cell clusters during the early stages of vasculogenesis; (iii) mechanical pattern formation under different linear viscoelasticity assumptions for the extracellular matrix. The mathematical models proposed for these problems are formulated as systems of nonlinear and nonlocal PDEs, which provide a mean-field representation of the underlying cellular dynamics and pose a series of interesting analytical and numerical challenges. These are tackled by means of formal asymptotic methods, linear stability analyses and appropriate numerical schemes preventing the emergence of spurious oscillations. Numerical simulations, relying on parameter values drawn from the extant literature, complement the analytical results and are employed for in silico investigations qualitatively testing the model assumptions against empirical observations. The obtained results help us shed light on the hidden mechanisms responsible for the emergence of the studied phenomena in biology and medicine, suggesting promising research perspectives.

Publicly available code

Note: open access versions of preprints and published articles are available on HAL and ArXiv or bioRxiv