Simulating BRAFV600E-MEK-ERK signalling dynamics in response to vertical inhibition treatment strategies.

In vertical inhibition treatment strategies, multiple components of an intracellular pathway are simultaneously inhibited. Vertical inhibition of the BRAFV600E-MEK-ERK signalling pathway is a standard of care for treating BRAFV600E-mutated melanoma where two targeted cancer drugs, a BRAFV600E-inhibitor, and a MEK inhibitor, are administered in combination. Targeted therapies have been linked to early onsets of drug resistance, and thus treatment strategies of higher complexities and lower doses have been proposed as alternatives to current clinical strategies. However, finding optimal complex, low-dose treatment strategies is a challenge, as it is possible to design more treatment strategies than are feasibly testable in experimental settings. To quantitatively address this challenge, we develop a mathematical model of BRAFV600E-MEK-ERK signalling dynamics in response to combinations of the BRAFV600E-inhibitor dabrafenib (DBF), the MEK inhibitor trametinib (TMT), and the ERK-inhibitor SCH772984 (SCH). From a model of the BRAFV600E-MEK-ERK pathway, and a set of molecular-level drug-protein interactions, we extract a system of chemical reactions that is parameterised by in vitro data and converted to a system of ordinary differential equations (ODEs) using the law of mass action. The ODEs are solved numerically to produce simulations of how pathway-component concentrations change over time in response to different treatment strategies, i.e., inhibitor combinations and doses. The model can thus be used to limit the search space for effective treatment strategies that target the BRAFV600E-MEK-ERK pathway and warrant further experimental investigation. The results demonstrate that DBF and DBF-TMT-SCH therapies show marked sensitivity to BRAFV600E concentrations in silico, whilst TMT and SCH monotherapies do not.

Stochastic differential equation modelling of cancer cell migration and tissue invasion.

Invasion of the surrounding tissue is a key aspect of cancer growth and spread involving a coordinated effort between cell migration and matrix degradation, and has been the subject of mathematical modelling for almost 30 years. In this current paper we address a long-standing question in the field of cancer cell migration modelling. Namely, identify the migratory pattern and spread of individual cancer cells, or small clusters of cancer cells, when the macroscopic evolution of the cancer cell colony is dictated by a specific partial differential equation (PDE). We show that the usual heuristic understanding of the diffusion and advection terms of the PDE being one-to-one responsible for the random and biased motion of the solitary cancer cells, respectively, is not precise. On the contrary, we show that the drift term of the correct stochastic differential equation scheme that dictates the individual cancer cell migration, should account also for the divergence of the diffusion of the PDE. We support our claims with a number of numerical experiments and computational simulations.

A novel nonlocal partial differential equation model of endothelial progenitor cell cluster formation during the early stages of vasculogenesis.

The formation of new vascular networks 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, known as 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 and solve the equations numerically. We then conduct a parametric analysis of the numerical solutions of the one-dimensional problem to investigate the role of underlying dynamics on the speed of cluster formation and the size of clusters, measured via appropriate metrics for the cluster width and compactness. We verify the key results of the parametric analysis with simulations of the two-dimensional 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, suggesting chemotaxis is responsible for the cluster topology while matrix degradation is responsible for the speed of cluster formation. Our results also indicate that the nonlocal cell-to-cell adhesion term in our model, even though it initially causes cells to aggregate, is not sufficient to ensure clusters are stable over long time periods. Consequently, new modelling strategies for cell-to-cell adhesion are required to stabilise in silico clusters. We end the paper with a thorough discussion of promising, fruitful future modelling and experimental research perspectives.

Quantifying ERK activity in response to inhibition of the BRAFV600E-MEK-ERK cascade using mathematical modelling.

BACKGROUND: Simultaneous inhibition of multiple components of the BRAF-MEK-ERK cascade (vertical inhibition) has become a standard of care for treating BRAF-mutant melanoma. However, the molecular mechanism of how vertical inhibition synergistically suppresses intracellular ERK activity, and consequently cell proliferation, are yet to be fully elucidated. METHODS: We develop a mechanistic mathematical model that describes how the mutant BRAF inhibitor, dabrafenib, and the MEK inhibitor, trametinib, affect BRAFV600E-MEK-ERK signalling. The model is based on a system of chemical reactions that describes cascade signalling dynamics. Using mass action kinetics, the chemical reactions are re-expressed as ordinary differential equations that are parameterised by in vitro data and solved numerically to obtain the temporal evolution of cascade component concentrations. RESULTS: The model provides a quantitative method to compute how dabrafenib and trametinib can be used in combination to synergistically inhibit ERK activity in BRAFV600E-mutant melanoma cells. The model elucidates molecular mechanisms of vertical inhibition of the BRAFV600E-MEK-ERK cascade and delineates how elevated BRAF concentrations generate drug resistance to dabrafenib and trametinib. The computational simulations further suggest that elevated ATP levels could be a factor in drug resistance to dabrafenib. CONCLUSIONS: The model can be used to systematically motivate which dabrafenib-trametinib dose combinations, for treating BRAFV600E-mutated melanoma, warrant experimental investigation.

Targeting Cellular DNA Damage Responses in Cancer: An In Vitro-Calibrated Agent-Based Model Simulating Monolayer and Spheroid Treatment Responses to ATR-Inhibiting Drugs.

We combine a systems pharmacology approach with an agent-based modelling approach to simulate LoVo cells subjected to AZD6738, an ATR (ataxia-telangiectasia-mutated and rad3-related kinase) inhibiting anti-cancer drug that can hinder tumour proliferation by targeting cellular DNA damage responses. The agent-based model used in this study is governed by a set of empirically observable rules. By adjusting only the rules when moving between monolayer and multi-cellular tumour spheroid simulations, whilst keeping the fundamental mathematical model and parameters intact, the agent-based model is first parameterised by monolayer in vitro data and is thereafter used to simulate treatment responses in in vitro tumour spheroids subjected to dynamic drug delivery. Spheroid simulations are subsequently compared to in vivo data from xenografts in mice. The spheroid simulations are able to capture the dynamics of in vivo tumour growth and regression for approximately 8 days post-tumour injection. Translating quantitative information between in vitro and in vivo research remains a scientifically and financially challenging step in preclinical drug development processes. However, well-developed in silico tools can be used to facilitate this in vitro to in vivo translation, and in this article, we exemplify how data-driven, agent-based models can be used to bridge the gap between in vitro and in vivo research. We further highlight how agent-based models, that are currently underutilised in pharmaceutical contexts, can be used in preclinical drug development.

Calibrating models of cancer invasion: parameter estimation using approximate Bayesian computation and gradient matching.

We present two different methods to estimate parameters within a partial differential equation model of cancer invasion. The model describes the spatio-temporal evolution of three variables-tumour cell density, extracellular matrix density and matrix degrading enzyme concentration-in a one-dimensional tissue domain. The first method is a likelihood-free approach associated with approximate Bayesian computation; the second is a two-stage gradient matching method based on smoothing the data with a generalized additive model (GAM) and matching gradients from the GAM to those from the model. Both methods performed well on simulated data. To increase realism, additionally we tested the gradient matching scheme with simulated measurement error and found that the ability to estimate some model parameters deteriorated rapidly as measurement error increased.

Mechanical Models of Pattern and Form in Biological Tissues: The Role of Stress-Strain Constitutive Equations.

Mechanical and mechanochemical models of pattern formation in biological tissues have been used to study a variety of biomedical systems, particularly in developmental biology, and describe the physical interactions between cells and their local surroundings. These models in their original form consist of a balance equation for the cell density, a balance equation for the density of the extracellular matrix (ECM), and a force-balance equation describing the mechanical equilibrium of the cell-ECM system. Under the assumption that the cell-ECM 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 pattern formation process and the type of biological tissue considered, other constitutive models of linear viscoelasticity may be better suited. In this paper, 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 and the dispersion relations derived therefrom support the idea that fluid-like constitutive models, such as the Maxwell model and the Jeffrey model, have a pattern formation potential much higher than solid-like models, such as the Kelvin-Voigt model and the standard linear solid model. This is confirmed by the results of numerical simulations, which demonstrate that, all else being equal, spatial patterns emerge in the case where the Maxwell model is used to represent the stress-strain relation of the ECM, while no patterns are observed when the Kelvin-Voigt model is employed. 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 mechanical and mechanochemical models of pattern formation with stress-strain constitutive equations for the ECM that provide a more faithful representation of the underlying tissue rheology.

A novel 3D atomistic-continuum cancer invasion model: In silico simulations of an in vitro organotypic invasion assay.

We develop a three-dimensional genuinely hybrid atomistic-continuum model that describes the invasive growth dynamics of individual cancer cells in tissue. The framework explicitly accounts for phenotypic variation by distinguishing between cancer cells of an epithelial-like and a mesenchymal-like phenotype. It also describes mutations between these cell phenotypes in the form of epithelial-mesenchymal transition (EMT) and its reverse process mesenchymal-epithelial transition (MET). The proposed model consists of a hybrid system of partial and stochastic differential equations that describe the evolution of epithelial-like and mesenchymal-like cancer cells, respectively, under the consideration of matrix-degrading enzyme concentrations and the extracellular matrix density. With the help of inverse parameter estimation and a sensitivity analysis, this three-dimensional model is then calibrated to an in vitro organotypic invasion assay experiment of oral squamous cell carcinoma cells.

Development of a coupled simulation toolkit for computational radiation biology based on Geant4 and CompuCell3D.

Understanding and designing clinical radiation therapy is one of the most important areas of state-of-the-art oncological treatment regimens. Decades of research have gone into developing sophisticated treatment devices and optimization protocols for schedules and dosages. In this paper, we presented a comprehensive computational platform that facilitates building of the sophisticated multi-cell-based model of how radiation affects the biology of living tissue. We designed and implemented a coupled simulation method, including a radiation transport model, and a cell biology model, to simulate the tumor response after irradiation. The radiation transport simulation was implemented through Geant4 which is an open-source Monte Carlo simulation platform that provides many flexibilities for users, as well as low energy DNA damage simulation physics, Geant4-DNA. The cell biology simulation was implemented using CompuCell3D (CC3D) which is a cell biology simulation platform. In order to couple Geant4 solver with CC3D, we developed a 'bridging' module, RADCELL, that extracts tumor cellular geometry of the CC3D simulation (including specification of the individual cells) and ported it to the Geant4 for radiation transport simulation. The cell dose and cell DNA damage distribution in multicellular system were obtained using Geant4. The tumor response was simulated using cell-based tissue models based on CC3D, and the cell dose and cell DNA damage information were fed back through RADCELL to CC3D for updating the cell properties. By merging two powerful and widely used modeling platforms, CC3D and Geant4, we delivered a novel tool that can give us the ability to simulate the dynamics of biological tissue in the presence of ionizing radiation, which provides a framework for quantifying the biological consequences of radiation therapy. In this introductory methods paper, we described our modeling platform in detail and showed how it can be applied to study the application of radiotherapy to a vascularized tumor.

A hybrid discrete-continuum approach to model Turing pattern formation.

Since its introduction in 1952, with a further refinement in 1972 by Gierer and Meinhardt, Turing's (pre-)pattern theory (the chemical basis of morphogenesis) has been widely applied to a number of areas in developmental biology, where evolving cell and tissue structures are naturally observed. The related pattern formation models normally comprise a system of reaction-diffusion equations for interacting chemical species (morphogens), whose heterogeneous distribution in some spatial domain acts as a template for cells to form some kind of pattern or structure through, for example, differentiation or proliferation induced by the chemical pre-pattern. Here we develop a hybrid discrete-continuum modelling framework for the formation of cellular patterns via the Turing mechanism. In this framework, a stochastic individual-based model of cell movement and proliferation is combined with a reaction-diffusion system for the concentrations of some morphogens. As an illustrative example, we focus on a model in which the dynamics of the morphogens are governed by an activator-inhibitor system that gives rise to Turing pre-patterns. The cells then interact with the morphogens in their local area through either of two forms of chemically-dependent cell action: Chemotaxis and chemically-controlled proliferation. We begin by considering such a hybrid model posed on static spatial domains, and then turn to the case of growing domains. In both cases, we formally derive the corresponding deterministic continuum limit and show that that there is an excellent quantitative match between the spatial patterns produced by the stochastic individual-based model and its deterministic continuum counterpart, when sufficiently large numbers of cells are considered. This paper is intended to present a proof of concept for the ideas underlying the modelling framework, with the aim to then apply the related methods to the study of specific patterning and morphogenetic processes in the future.

Learning-induced switching costs in a parasitoid can maintain diversity of host aphid phenotypes although biocontrol is destabilized under abiotic stress.

Aphid populations frequently include phenotypes that are resistant to parasitism by hymenopterous parasitoid wasps, which is often attributed to the presence of 'protective' facultative endosymbionts residing in aphid tissues, particularly Hamiltonella defensa. In field conditions, under parasitoid pressure, the observed coexistence of aphids with and without protective symbionts cannot be explained by their difference in fitness alone. Using the cereal aphid Rhopalosiphum padi as a model, we propose an alternative mechanism whereby parasitoids are more efficient at finding common phenotypes of aphid and experience a fitness cost when switching to the less common phenotype. We construct a model based on delay differential equations and parameterize and validate the model with values within the ranges obtained from experimental studies. We then use it to explore the possible effects on system dynamics under conditions of environmental stress, using our existing data on the effects of drought stress in crops as an example. We show the 'switching penalty' incurred by parasitoids leads to stable coexistence of aphids with and without H. defensa and provides a potential mechanism for maintaining phenotypic diversity among host organisms. We show that drought-induced reduction in aphid development time has little impact. However, greater reduction in fecundity on droughted plants of symbiont-protected aphids can cause insect population cycles when the system would be stable in the absence of drought stress. The stabilizing effect of the increased efficiency in dealing with more commonly encountered host phenotypes is applicable to a broad range of consumer-resource systems and could explain stable coexistence in competitive environments. The loss of stable coexistence when drought has different effects on the competing aphid phenotypes highlights the importance of scenario testing when considering biocontrol for pest management.

A multiscale model of virus pandemic: Heterogeneous interactive entities in a globally connected world.

This paper is devoted to the multidisciplinary modelling of a pandemic initiated by an aggressive virus, specifically the so-called SARS-CoV-2 Severe Acute Respiratory Syndrome, corona virus n.2. The study is developed within a multiscale framework accounting for the interaction of different spatial scales, from the small scale of the virus itself and cells, to the large scale of individuals and further up to the collective behaviour of populations. An interdisciplinary vision is developed thanks to the contributions of epidemiologists, immunologists and economists as well as those of mathematical modellers. The first part of the contents is devoted to understanding the complex features of the system and to the design of a modelling rationale. The modelling approach is treated in the second part of the paper by showing both how the virus propagates into infected individuals, successfully and not successfully recovered, and also the spatial patterns, which are subsequently studied by kinetic and lattice models. The third part reports the contribution of research in the fields of virology, epidemiology, immune competition, and economy focussed also on social behaviours. Finally, a critical analysis is proposed looking ahead to research perspectives.

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

Continuum models for the spatial dynamics of growing cell populations have been widely used to investigate the mechanisms underpinning tissue development and tumour invasion. These models consist of nonlinear partial differential equations that describe the evolution of cellular densities in response to pressure gradients generated by population growth. Little prior work has explored the relation between such continuum models and related single-cell-based models. We present here a simple stochastic individual-based model for the spatial dynamics of multicellular systems whereby cells undergo pressure-driven movement and pressure-dependent proliferation. We show that nonlinear partial differential equations commonly used to model the spatial dynamics of growing cell populations can be formally derived from the branching random walk that underlies our discrete model. Moreover, we carry out a systematic comparison between the individual-based model and its continuum counterparts, both in the case of one single cell population and in the case of multiple cell populations with different biophysical properties. The outcomes of our comparative study demonstrate that the results of computational simulations of the individual-based model faithfully mirror the qualitative and quantitative properties of the solutions to the corresponding nonlinear partial differential equations. Ultimately, these results illustrate how the simple rules governing the dynamics of single cells in our individual-based model can lead to the emergence of complex spatial patterns of population growth observed in continuum models.

A Mathematical Framework for Modelling the Metastatic Spread of Cancer.

Cancer is a complex disease that starts with mutations of key genes in one cell or a small group of cells at a primary site in the body. If these cancer cells continue to grow successfully and, at some later stage, invade the surrounding tissue and acquire a vascular network, they can spread to distant secondary sites in the body. This process, known as metastatic spread, is responsible for around 90% of deaths from cancer and is one of the so-called hallmarks of cancer. To shed light on the metastatic process, we present a mathematical modelling framework that captures for the first time the interconnected processes of invasion and metastatic spread of individual cancer cells in a spatially explicit manner-a multigrid, hybrid, individual-based approach. This framework accounts for the spatiotemporal evolution of mesenchymal- and epithelial-like cancer cells, membrane-type-1 matrix metalloproteinase (MT1-MMP) and the diffusible matrix metalloproteinase-2 (MMP-2), and for their interactions with the extracellular matrix. Using computational simulations, we demonstrate that our model captures all the key steps of the invasion-metastasis cascade, i.e. invasion by both heterogeneous cancer cell clusters and by single mesenchymal-like cancer cells; intravasation of these clusters and single cells both via active mechanisms mediated by matrix-degrading enzymes (MDEs) and via passive shedding; circulation of cancer cell clusters and single cancer cells in the vasculature with the associated risk of cell death and disaggregation of clusters; extravasation of clusters and single cells; and metastatic growth at distant secondary sites in the body. By faithfully reproducing experimental results, our simulations support the evidence-based hypothesis that the membrane-bound MT1-MMP is the main driver of invasive spread rather than diffusible MDEs such as MMP-2.

Blackboard to Bedside: A Mathematical Modeling Bottom-Up Approach Toward Personalized Cancer Treatments.

Cancers present with high variability across patients and tumors; thus, cancer care, in terms of disease prevention, detection, and control, can highly benefit from a personalized approach. For a comprehensive personalized oncology practice, this personalization should ideally consider data gathered from various information levels, which range from the macroscale population level down to the microscale tumor level, without omission of the central patient level. Appropriate data mined from each of these levels can significantly contribute in devising personalized treatment plans tailored to the individual patient and tumor. Mathematical models of solid tumors, combined with patient-specific tumor profiles, present a unique opportunity to personalize cancer treatments after detection using a bottom-up approach. Here, we discuss how information harvested from mathematical models and from corresponding in silico experiments can be implemented in preclinical and clinical applications. To conceptually illustrate the power of these models, one such model is presented, and various pertinent tumor and treatment scenarios are demonstrated in silico. The presented model, specifically a multiscale, hybrid cellular automaton, has been fully validated in vitro using multiple cell-line-specific data. We discuss various insights provided by this model and other models like it and their role in designing predictive tools that are both patient, and tumor specific. After refinement and parametrization with appropriate data, such in silico tools have the potential to be used in a clinical setting to aid in treatment protocols and decision making.

Spatial-Stochastic modelling of synthetic gene regulatory networks.

Transcription factors are important molecules which control the levels of mRNA and proteins within cells by modulating the process of transcription (the mechanism by which mRNA is produced within cells) and hence translation (the mechanism by which proteins are produced within cells). Transcription factors are part of a wider family of molecular interaction networks known as gene regulatory networks (GRNs) which play an important role in key cellular processes such as cell division and apoptosis (e.g. the p53-Mdm2, NFκB pathways). Transcription factors exert control over molecular levels through feedback mechanisms, with proteins binding to gene sites in the nucleus and either up-regulating or down-regulating production of mRNA. In many GRNs, there is a negative feedback in the network and the transcription rate is reduced. Typically, this leads to the mRNA and protein levels oscillating over time and also spatially between the nucleus and cytoplasm. When experimental data for such systems is analysed, it is observed to be noisy and in many cases the actual numbers of molecules involved are quite low. In order to model such systems accurately and connect with the data in a quantitative way, it is therefore necessary to adopt a stochastic approach as well as take into account the spatial aspect of the problem. In this paper, we extend previous work in the area by formulating and analysing stochastic spatio-temporal models of synthetic GRNs e.g. repressilators and activator-repressor systems.

Quantitative Predictive Modelling Approaches to Understanding Rheumatoid Arthritis: A Brief Review.

Rheumatoid arthritis is a chronic autoimmune disease that is a major public health challenge. The disease is characterised by inflammation of synovial joints and cartilage erosion, which lead to chronic pain, poor life quality and, in some cases, mortality. Understanding the biological mechanisms behind the progression of the disease, as well as developing new methods for quantitative predictions of disease progression in the presence/absence of various therapies is important for the success of therapeutic approaches. The aim of this study is to review various quantitative predictive modelling approaches for understanding rheumatoid arthritis. To this end, we start by briefly discussing the biology of this disease and some current treatment approaches, as well as emphasising some of the open problems in the field. Then, we review various mathematical mechanistic models derived to address some of these open problems. We discuss models that investigate the biological mechanisms behind the progression of the disease, as well as pharmacokinetic and pharmacodynamic models for various drug therapies. Furthermore, we highlight models aimed at optimising the costs of the treatments while taking into consideration the evolution of the disease and potential complications.

The role of spatial variations of abiotic factors in mediating intratumour phenotypic heterogeneity.

We present here a space- and phenotype-structured model of selection dynamics between cancer cells within a solid tumour. In the framework of this model, we combine formal analyses with numerical simulations to investigate in silico the role played by the spatial distribution of abiotic components of the tumour microenvironment in mediating phenotypic selection of cancer cells. Numerical simulations are performed both on the 3D geometry of an in silico multicellular tumour spheroid and on the 3D geometry of an in vivo human hepatic tumour, which was imaged using computerised tomography. The results obtained show that inhomogeneities in the spatial distribution of oxygen, currently observed in solid tumours, can promote the creation of distinct local niches and lead to the selection of different phenotypic variants within the same tumour. This process fosters the emergence of stable phenotypic heterogeneity and supports the presence of hypoxic cells resistant to cytotoxic therapy prior to treatment. Our theoretical results demonstrate the importance of integrating spatial data with ecological principles when evaluating the therapeutic response of solid tumours to cytotoxic therapy.

Modelling the Immune Response to Cancer: An Individual-Based Approach Accounting for the Difference in Movement Between Inactive and Activated T Cells.

A growing body of experimental evidence indicates that immune cells move in an unrestricted search pattern if they are in the pre-activated state, whilst they tend to stay within a more restricted area upon activation induced by the presence of tumour antigens. This change in movement is not often considered in the existing mathematical models of the interactions between immune cells and cancer cells. With the aim to fill such a gap in the existing literature, in this work we present a spatially structured individual-based model of tumour-immune competition that takes explicitly into account the difference in movement between inactive and activated immune cells. In our model, a Lévy walk is used to capture the movement of inactive immune cells, whereas Brownian motion is used to describe the movement of antigen-activated immune cells. The effects of activation of immune cells, the proliferation of cancer cells and the immune destruction of cancer cells are also modelled. We illustrate the ability of our model to reproduce qualitatively the spatial trajectories of immune cells observed in experimental data of single-cell tracking. Computational simulations of our model further clarify the conditions for the onset of a successful immune action against cancer cells and may suggest possible targets to improve the efficacy of cancer immunotherapy. Overall, our theoretical work highlights the importance of taking into account spatial interactions when modelling the immune response to cancer cells.