Integrating Quantitative Assays with Biologically Based Mathematical Modeling for Predictive Oncology.

We provide an overview on the use of biological assays to calibrate and initialize mechanism-based models of cancer phenomena. Although artificial intelligence methods currently dominate the landscape in computational oncology, mathematical models that seek to explicitly incorporate biological mechanisms into their formalism are of increasing interest. These models can guide experimental design and provide insights into the underlying mechanisms of cancer progression. Historically, these models have included a myriad of parameters that have been difficult to quantify in biologically relevant systems, limiting their practical insights. Recently, however, there has been much interest calibrating biologically based models with the quantitative measurements available from (for example) RNA sequencing, time-resolved microscopy, and in vivo imaging. In this contribution, we summarize how a variety of experimental methods quantify tumor characteristics from the molecular to tissue scales and describe how such data can be directly integrated with mechanism-based models to improve predictions of tumor growth and treatment response.

Introduction

Although there is an ever-increasing body of clinical research that has sought to characterize cancer development and evaluate viable therapies, there is still a fundamental lack of understanding of underlying factors involved in cancer development and progression. This is no doubt owing to the complex and dynamic interplay of heterogeneous cell populations, signaling factors, and microenvironments that exist within tumors (Hanahan and Coussens, 2012; Hanahan and Weinberg, 2011). In an attempt to capture and characterize these phenomena in a well-controlled setting, scientists have developed and utilized a myriad of in vitro and in vivo biological assays to uncover the molecular and cellular drivers of tumor growth and treatment response. However, the vast majority of analyses applied to these in vitro and in vivo assays reduce the observations to a single measurement or snapshot in time, ignoring the complex spatially and temporally evolving nature of cancer. For example, angiogenesis is linked to mechanical properties of the tumor microenvironment, cell-cell interactions between tumor and vascular cells, as well as the presence of soluble growth factors (De Palma et al., 2017). However, common angiogenic assays, such as tube formation assays, reduce this phenomena to a measurement of expression of proangiogenic cytokines (e.g., vascular endothelial growth factor) or endothelial cell proliferation (Staton et al., 2004). Although 2D and 3D models enable the study of more complex interactions (e.g., cell-cell or cell-ECM interactions [Boussommier-Calleja, 2020; Gadde et al., 2018; Hoarau-Véchot et al., 2018; Michna et al., 2018]) and in vivo models can investigate tumor progression in the context of the microenvironment, the ability to quantitatively characterize the spatiotemporal evolution of cancer is quite limited within current in vitro and in vivo experimental paradigms.

We (and others) have posited that without a mathematical theory enabling the systematic analysis of the numerous spatially and temporally evolving parameters and variables that characterize cancer, we are left with trial and error (Coveney et al., 2016; Yankeelov et al, 2013, 2015). Mathematical models based on biological processes can be used to abstract the fundamental phenomena of a developing tumor and allow for “in silico” experimentation to systematically investigate cancer processes. Biological assays involve large investments in both time and money and often entail numerous conditions to evaluate quantities of interest against appropriate controls. Mathematical models can guide experimental design by simulating effects of the phenomena under investigation and narrowing the scope of biological measurements needed, thereby saving resources. Additionally, they can provide mechanistic insight into cancer processes that are difficult to isolate and measure experimentally. Development and validation of these mathematical models with appropriate biological measurements may yield a comprehensive theory of cancer across spatial and temporal scales that can be used to fundamentally shift how cancers are treated. Furthermore, application of biologically based, mathematical models that can be parameterized for individual patients—and not “trained” on large population datasets—may yield patient-specific predictions and have the potential to optimize therapeutic interventions on an individual patient basis (Rockne et al., 2019; Yankeelov et al, 2013, 2015).

Historically there has been a lack of integration of mathematical theory and experimental data in cancer biology and oncology (Yankeelov et al., 2013). We posit that this critical obstacle to personalized predictions can now be overcome with quantitative biological assays combined with appropriate mathematical models designed to incorporate them. One fundamental issue in attempting to bridge this divide is that biological assays have not always yielded quantitative measures that could be used to parameterize mechanism-based models of cancer progression and response. Furthermore, the existing mathematical models were not structured in a way that allowed for incorporation of data from such assays. However, there has been much recent growth in developing assays that can provide quantitative biological data, as well as new efforts to develop mechanism-based mathematical models that enable incorporation of available biological data (Rockne et al., 2019).

In this review, we introduce biological assays appropriate for quantitatively probing cancer dynamics in vitro, ex vivo, and in vivo and highlight the opportunities for their integration into mechanism-based, mathematical models of tumor growth and treatment response at the subcellular, cellular, and tissue scales (please see Figure 1). This developing field is illustrated by presenting representative examples of studies that have integrated quantitative data into mathematical models of solid tumor malignancies. We also discuss the experimental and computational limitations to modeling cancer at each scale to highlight the need for robust integration between mathematical and biological experiments.

Figure 1

Overview of Biological Scales Described in This Review

For each scale we will introduce key elements of the relevant tumor biology, as well as how experimental data are integrated with mathematical models to study cancer. At the cellular scale (A), we focus on the development of mathematical models for identification of heterogeneous phenotypes and differentiation of cancer cell populations. At the microenvironmental and tissue scales (B), we focus on in vitro and ex vivo data driven mathematical models for analyzing protein and nutrient gradients, multicellular interactions, and vasculature within the tumor microenvironment. Finally, we explore the tissue and organ scales (C) in the in vivo setting by discussing the coupling of imaging data from animal models and patient tumors to mathematical models for predicting tumor development.

Biology-Based Mathematical Models

Mathematical models of cancer are broadly classified as discrete or continuous. Discrete approaches usually consist of agent-based models (Metzcar et al., 2019; Wang et al., 2015), whereby tumor cells are represented as independent entities that follow a set of rules dictating cellular behavior. Thus, agent-based models are more amenable to in vitro cellular or small-range tissue-scale applications, which can utilize microscopically resolved, genomic, and transcriptomic data. Continuous formulations feature ordinary (ODE) or partial (PDE) differential equations that describe the population-averaged or spatially resolved dynamics of tumor cells, respectively. ODE models are limited to temporal dynamics, but this enables their use in a wide range of cellular and tissue-scale settings that can leverage multiple types of time-resolved data (Benzekry et al., 2014; Jarrett et al., 2019; Johnson et al., 2020; Lorenzo et al., 2019b; Mendoza-Juez et al., 2012; Morken et al., 2014). PDE models further include spatial dynamics and constitute the standard approach for clinical, tissue-scale applications, which usually rely on longitudinal in vivo imaging data (Hormuth et al., 2020; Jarrett et al., 2018; Lorenzo et al., 2019a; Mang et al., 2020; Rockne et al., 2015; Wong et al., 2017). Additionally, hybrid models combine a discrete and continuous formulation, which may require multiscale multimodal data (Metzcar et al., 2019; Phillips et al., 2020; Vilanova et al., 2017; Wang et al., 2015).

In principle, mathematical models of cancer often describe biological mechanisms common across tumor types and therefore have the potential for general application. However, in practice, these models often must be customized to describe disease-specific characteristics related to tumor development and response to therapy. Despite the wealth and diversity of biologically based, mathematical models of cancer, there is a paucity of studies that explicitly link model predictions with experimental data (Altrock et al., 2015; Karolak et al., 2018; Yankeelov et al., 2013). To facilitate future efforts at this interface, we now describe how state-of-the-art biological assays can provide quantitative data to calibrate and validate mathematical models from the cellular (following section) to tissue scales (next two sections).

Mathematical Modeling of Cancer at the Cellular Scale

Cellular Heterogeneity in Cancer

The origin and interactions of genetically and phenotypically distinct cells within a patient's tumor influence its proliferation rate, metastatic potential, and treatment response (Brock et al., 2009; Maley et al., 2017; Pisco and Huang, 2015). Mathematical modeling can be used to describe and predict cellular behavior, to better understand the contribution of individual subpopulations of cells to a cancer's dynamic ecosystem. Experimental assays can be used to quantify subcellular attributes (from genes to protein), as well as identify cellular phenotypes. In Section “Cellular Scale Data Appropriate for Incorporating into Mathematical Models of Cancer,” we will cover a diverse set of quantitative assays that can be used to unveil the cellular phenotypic states. These widely available biological techniques can be used to calibrate, initialize, and validate mathematical models that seek to unveil fundamental mechanisms of tumorigenesis and understand how cancer cell populations interact and evolve (Figure 2).

Figure 2

Quantification of Cellular Properties to Parameterize Mathematical Models of Cancer

(A) Illustration of how cancer cell populations exhibit heterogeneity at multiple levels, described by genomics (represented by color), morphology (represented by shape), and gene expression (represented by receptor status). Experimental assays quantify genomic, transcriptomic, or proteomic differences between cells to reveal the characteristic properties of tumor cell subpopulations.

(B) Example outputs from high-throughput assays such as (1) flow cytometry in which cell populations are identified by surface marker expression, (2) barcode labeling to quantify abundance of clonal populations over time, (3) scRNA-seq to measure differential gene expression in individual cells, and (4) reveal distinct phenotypic clusters of cells.

(C) Display of how biological processes of tumor cell populations can be represented by mathematical models describing the behavior and interactions of distinct cell types. Such models can be informed by the data in (B) to reveal novel insights about the dynamic behavior of the cell population.

Cellular Scale Data Appropriate for Incorporating into Mathematical Models of Cancer

Tumors are composed of phenotypically diverse cancer cell populations, which determine tumor behavior and therapeutic sensitivity and response (Hardeman et al., 2017; Pisco et al., 2013). Experimental assays quantifying gene expression and protein content can reveal key phenotypic states and signaling pathways, information that can then be used to parameterize mathematical models of cancer cell dynamics. Quantitative polymerase chain reaction (qPCR) can be used to quantify DNA levels within a sample and determine if a specific mutant allele is present in a cell population (Alvarez-Garcia et al., 2018). Using qPCR measurements from patients with chronic myeloid leukemia (CML), Michor et al. developed a four-compartment model of leukemic cell populations (Michor et al., 2005). Their model suggests that imatinib, a CML treatment, was effective treatment of differentiated leukemic cells but not of leukemic stem cells, yielding insight into CML treatment failure (Michor et al., 2005). DNA sequencing technologies, such as whole genome sequencing, can be used to quantify genetic mutations associated with cancer cell populations. Such data have been used to link levels of specific mutations with tumor growth patterns, to reveal which mutations are driving tumor behavior. For example, in patients with chronic lymphocytic leukemia, Gruber et al. tracked the accumulation of genetic mutations over time and applied Bayesian methods to estimate clonal proportions within the bulk population, then quantifying the fitness benefits conferred to the bulk tumor (Gruber et al., 2019). A number of stochastic and agent-based models have been developed to describe the effects of genetic drift and subclonal evolution on bulk tumor behavior (Beerenwinkel et al., 2007; Chowell et al., 2018; McFarland et al., 2013; Moolgavkar and Knudson, 1981; Waclaw et al., 2015), but there remains a need for integration of experimental data into these frameworks to construct models capable of making predictions that can be directly tested.

RNA sequencing (RNA-seq) data, in combination with bioinformatics techniques can reveal gene expression patterns that drive protein signaling, cellular behavior, and therapeutic response. For example, Bouhaddou et al. integrated high-dimensional multi-omics data, including patient-derived RNA-seq data, to construct an ODE model that integrated stochastic fluctuations in expression of commonly mutated cancer pathways and showed how these fluctuations influence tumor cell proliferation and death (Bouhaddou et al., 2018). With bulk RNA sequencing, gene expression patterns are averaged across cell populations. Single-cell RNA-seq (scRNA-seq) can be used to assess gene expression patterns at the individual cell level to identify unique cellular phenotypes, enabling characterization of the diversity within cancer cell populations (Becht et al., 2019). For example, Johnson et al. calibrated a model with estimates of subpopulation composition from scRNA-seq with longitudinal population data of cell number over time to predict the population response to new treatment regimens and assess the degree of drug-induced resistance (Johnson et al., 2020). Inclusion of scRNA-seq data yielded improved predictive accuracy of treatment response dynamics with a concordance correlation coefficient (CCC) of 0.89, compared with 0.79 without transcriptomic information (Johnson et al., 2020).

Protein expression on the surface of cancer cells allows quick identification of distinct phenotypes that can be monitored and isolated using flow cytometry and fluorescence-activated cell sorting (Fillmore and Kuperwasser, 2008; Gupta et al., 2011; Shipitsin et al., 2007). For example, Gupta et al. parameterized a Markov chain Monte Carlo model of phenotypic state transitions by monitoring the expression of characteristic surface proteins on isolated basal, luminal, and stem-cell like cells from breast cancer cell lines (Gupta et al., 2011). More recently, Morinishi et al. constructed an ODE model to describe cancer cell transitions between a more differentiated and a more stem-like state and experimentally estimated transition rates between the two cell states under different conditions with flow cytometric time course data (Morinishi et al., 2020). Modeling the temporal dynamics of key proteins that identify unique cancer cell states is critical to developing an understanding of the biological mechanisms of cancer. Integration of flow cytometric measurements into Markov and ODE models by Pisco et al. allowed them to reject the null hypothesis that drug resistance in leukemia arises purely from Darwinian selection by differential growth rates, lending support to the hypothesis that a drug-resistant cellular state is driven by chemotherapeutically activated cell state switching (Pisco et al., 2013).

Limitations of Using In Vitro Data for Modeling Cancer

Although there is ever-increasing availability and quality of quantitative biological data at the cellular scale, even the best designed experiments have intrinsic limitations owing to sample origin, choice of methodology, and complexity of data output. Targeted experimental assays quantify the expression of specific genes or proteins of interest, whereas global experimental assays have increased readouts, providing expression information for thousands of genes. Targeted experiments can define cell behavior and states for quantitative modeling but may limit the scope of the cellular states identified and oversimplify a complex state-space. Global experimental approaches can overcome this limitation by revealing the gene networks driving specific phenomena but add significant complexity in data analysis as it can be difficult to identify which signals are most relevant for describing the phenomena of interest. Similarly, experimental assays that pool groups of cells together (prior to, for example, quantifying gene or protein expression) are beneficial for investigating whole population differences but ignore the cellular diversity within groups of cells. This limitation is overcome in single cell analyses, including scRNA-seq and flow cytometry, which provide insight into the full distribution of individual cell phenotypes in a population, but do so at the cost of increased complexity and use of resources.

Cancer cell populations vary in both space and time, and mathematical models of these dynamics require data that are spatially and/or temporally resolved. Many in vitro assays (including those mentioned in section “Cellular Scale Data Appropriate for Incorporating into Mathematical Models of Cancer”) require analysis of experimental samples at specific time points, often in a destructive or disruptive manner. Although high-throughput assays can be performed at multiple time points, these temporal samplings may miss dynamics occurring more quickly than can be experimentally recovered. Time-resolved microscopy provides a means to overcome some of this limitation by enabling tracking of individual cells over time. This technology has been applied to quantify whole-cell population growth under various conditions and, with the addition of fluorescent (Weissman and Pan, 2015) or nucleotide (Al’Khafaji et al., 2018; Nguyen et al., 2014) labels, can enable tracking of distinct cell phenotypes and specific protein expression within cells. This information can then be used to inform models of population growth and treatment response, as well as cell state changes over time (Hsu et al., 2019; Johnson et al., 2019).

Mathematical Modeling of Cancer at the Microenvironmental and Tissue Scales

The Tumor Microenvironment and Cancer Development

Central to tumor growth is the interaction of cancer cells with non-malignant stromal cell populations and the local environment, collectively referred to as the tumor microenvironment (Hanahan and Coussens, 2012). Stromal cell populations, such as fibroblasts, immune infiltrates, or vascular endothelial cells, are recruited to the tumor site via the cytokine expression of cancer cells (Hanahan and Coussens, 2012; Junttila and de Sauvage, 2013). These cell populations impact the composition of the tumor microenvironment, from local extracellular matrix deposition to gradients in nutrient concentration, which subsequently influence cancer cell protein expression and behavior (Hanahan and Coussens, 2012; Junttila and de Sauvage, 2013). Thus, the tumor microenvironment is an important consideration when modeling tumor development. In this section, we highlight quantitative assays that characterize the spatial distribution of cancer and stromal cell populations as well as environmental quantities such as oxygen and growth factor concentration.

Microenvironmental and Tissue Scale Data Appropriate for Incorporating into Mathematical Models

Microscopy techniques, ranging from confocal and brightfield microscopy to whole slide imaging, offer quantitative in vitro and ex vivo measurements of many properties of the tumor microenvironment (Gurcan et al., 2009; Heindl et al., 2015; Saucedo et al., 2018). Used in combination with immunohistochemical or immunofluorescent staining, microscopy can measure the spatial distribution of different cancer and stromal cell types, as well as microenvironment components such as collagen fibers or microvasculature (Gurcan et al., 2009; Heindl et al., 2015; Saucedo et al., 2018). For example, immunohistochemical staining techniques can be used to identify the spatial distribution of proliferating or apoptotic cancer cell subpopulations, which can then be used to calibrate models of tumor growth. Confocal imaging of normal and mutant breast epithelial spheroids was employed by Rejniak et al. to stain for caspase-3 (apoptosis) and Ki-67 (proliferation) to calibrate an in silico model of tumor spheroid formation (Rejniak et al., 2010). The model was then used to determine potential mechanisms of morphological and molecular changes that lead to development of mutant breast epithelial cells (Rejniak et al., 2010). Ahmadzadeh et al. employed two-photon microscopy to measure spatial variations in collagen distribution within the extracellular matrix surrounding tumor spheroids (Ahmadzadeh et al., 2017). These data were then used to calibrate a mechanochemical free-energy based model that sought to describe the feedback loop between cell contractility and realignment of the extracellular matrix fibers (Ahmadzadeh et al., 2017).

Techniques coupling immunofluorescent staining with time-resolved microscopy can be used to longitudinally track fluorescently labeled cells. Studies have applied such in vitro data to develop in silico models of tumor growth and invasion (Cristini et al., 2005; Enmon et al., 2001; Frieboes et al., 2006; Kam et al., 2009; Liedekerke et al., 2019; Lorenzo et al., 2016; Macklin and Lowengrub, 2007; Sanga et al., 2007), cellular migration driven by the tumor microenvironment (Kim et al., 2018; Mark et al., 2018; Shamloo et al., 2016; Stokes et al., 1991), vasculogenesis (Merks et al., 2006; Serini et al., 2003; Shamloo et al., 2016), and drug response (Kozusko et al., 2001; Maffei et al., 2014; Raghavan et al., 2016). In particular, our group has developed a 3D in vitro vascularized microfluidic tumor platform to study angiogenic sprouting of an endothelial vessel (Gadde et al., 2020). We utilized confocal microscopy images of vessel sprouting in combination with an agent-based model (Phillips et al., 2020) to determine parameter values that govern the driving dynamics of the model (e.g., vessel growth rate, production rate of pro-angiogenic factors) (Figure 3). This calibrated model can then be used to make predictions of vascular structure in response to growth factors secreted by different tumor cell distributions in the in vitro platform.

Figure 3

Model Calibration and Validation Pipeline

The left-most column shows the evolution of our agent-based model of tumor angiogenesis and the measured vessel segmented from confocal microscopy images (right column). We utilize the microscopy data from days 4 and 8 to calibrate key model parameters (e.g., production rate of VEGF) and then predict forward to day 12 where we can directly compare the predictions to the experimentally observed data. If the model is found to be invalid (i.e., if the predicted vasculature at day 12 does not appropriately recapitulate the observed data), the model must be modified by amending terms, parameters, or rules based on established biological principles. Once the model is calibrated and passes a validity test, the vasculature at day 16 is predicted with the calibrated parameters and compared with the data.

Histopathology analysis yields high-resolution brightfield or immunofluorescent images of tissue sections ex vivo (Gurcan et al., 2009). These sections are stained for markers of biological interest, including (for example) vessel maturation, immune infiltration, and proliferation (Gurcan et al., 2009; Heindl et al., 2015). Such measures can then be used to parameterize or validate models of tumor growth. Macklin et al. developed a patient-specific, agent-based model parameterized by clinically available histopathology data and accurately predicted the growth of ductal carcinoma in situ based on underlying cellular phenomena (Macklin et al., 2012). Similarly, Jarrett et al. developed a model of HER2+ cancer that accounted for the immune response to targeted therapy. In their model, they employed histological data to assign values for tumor necrosis within the model and then compared the predicted immune response with immunofluorescent images of immune infiltration (Jarrett et al., 2019). This mathematical-experimental approach was able to accurately model tumor response, with a CCC = 0.90 for treated tumor volumes, as well as identify interactions between tumor response, immune cells, and vascular development, interactions that would be challenging to elucidate only through experimentation.

Although immunohistochemistry staining of in vitro tumor models and ex vivo tumor samples can be used to quantify nutrient, cytokine, or drug distributions across a tissue, ELISAs, western blots, and protein arrays offer more accessible acquisition of similar information from bulk tissue samples. By employing time resolved microscopy and cytokine multiplex array in a microfluidic tumor model, Lee et al. developed an in silico model that couples chemokine-mediated signaling with mechanosensing to determine the influence of interstitial flow on migration of macrophages in the tumor microenvironment (Lee et al., 2020). Göttlich et al. used in vitro data from immunofluorescence, ELISA, western blot, and phosphokinase arrays from a 3D lung model to generate cell specific in silico models for use in drug effect prediction and patient stratification in treatment of lung cancer (Göttlich et al., 2018). A similar approach was utilized by Baur et al. for colorectal cancer application (Baur et al., 2019).

Limitations of Using In Vitro and Ex Vivo Data for Modeling Cancer

Calibrating biologically based, mathematical models often requires data obtained at multiple time points, whereas many biological assays measure quantities of interest at few or single time points. This limitation can be overcome, in part, by employing dynamic assays, such as time-resolved microscopy, but there are still constraints. For example, fluorescent assays often stain for only one or two markers of interest at a time and the continuous exposure of cells to the lasers used in confocal microscopy can result in cell death via phototoxicity. Additionally, biological data points defined using small sample sizes will have greater uncertainty in their measurements. For example, tumor heterogeneity provides a significant challenge in the calibration of tissue scale tumor models, as many mathematical parameters (e.g., cell growth rate) can vary depending on cell type, and heterogeneity can exist within cell populations of a cell line (Almendro et al., 2013). Furthermore, bridging the gap between the in vitro and in vivo settings is non-trivial owing to the fundamental differences in in vitro and in vivo experimental design. Thus, it is an open problem in the field as to how mathematical models that are calibrated with in vitro data can provide insight into in vivo systems. Adaption of newer technologies, such as microfluidic tumor- or organ-on-a-chip systems and 3D bioprinting, present means by which researchers can perform in vitro and ex vivo work in settings that are more comparable with the in vivo system. In particular, microfluidic platforms have revealed insight into metastatic cascade, cell migration, transport of anticancer agents, and even in predictive application for determination of tumor response to cancer treatments (Belgodere et al., 2018; Roberts et al., 2019; Sun et al., 2019; Trujillo-de Santiago et al., 2019). 3D bioprinting allows for fine control over the recreation of the tumor microenvironment (including the structural and compositional heterogeneity [Belgodere et al., 2018; Heinrich et al., 2019; Langer et al., 2019; Ma et al., 2018]) resulting in in vitro models that are realistic representations of an in vivo tumor.

Mathematical Modeling of Cancer at the Tissue and Organ Scales

Cancer Physiology at the Tissue and Organ Scale

Although the spatial scales of in vivo assessments of tumor physiology are typically limited to 0.25–5 mm (Frangioni, 2008; Ramamonjisoa and Ackerstaff, 2017), these techniques do enable measurement of temporal dynamics within the tumor microenvironment (Yankeelov et al., 2016b), such as vascular perfusion (Yankeelov et al., 2014; Yankeelov and Gore, 2009), metabolic activity(Castell and Cook, 2008), or interstitial flow (Kingsmore et al., 2018; LoCastro et al., 2020). These data types can also be acquired in 3D to evaluate an entire tumor volume, permitting investigation of regional differences within a tumor (e.g., regions of necrosis or proliferation, gradients of oxygen or nutrients) (Kim and Gatenby, 2017; Syed et al., 2020). Additionally, in vivo measurements can be acquired longitudinally within the same subject, thereby providing temporally and spatially resolved data critical for the calibration of spatiotemporal models of tumor progression (Hormuth et al., 2019b). The next section introduces several established techniques for imaging of tissue and incorporating such data into biologically based models of cancer.

Tissue and Organ Scale Data Appropriate for Incorporating into Mathematical Models

Magnetic resonance imaging (MRI), positron emission tomography (PET), and X-ray computed tomography (CT) are in vivo imaging techniques that characterize a myriad of tumor features. These include measures of size and shape (used clinically to assess treatment response [Nishino et al., 2010]), as well as more advanced measurements of cellularity (Padhani et al., 2009), vascularity (Yankeelov and Gore, 2009), and metabolic activity (Castell and Cook, 2008; Rajendran and Krohn, 2015) that can report more specifically on tumor status. As these data can provide longitudinal and 3D evaluations of tumors, they are excellent candidates for initializing, parameterizing, and validating mechanism-based models of tumor growth (Hormuth et al., 2019a).

Contrast-enhanced CT and MRI techniques (i.e., studies in which a contrast agent is injected into the patient to highlight differences between diseased and healthy tissue) as well as anatomical MRI (Figure 4A) are considered anatomical imaging techniques. These data types are commonly used to assess the extent of clinical (or bulk) tumor burden (e.g., longest dimension and volume), subclinical (or invasive) tumor burden, and vasogenic edema (Eisenhauer et al., 2009; Mabray et al., 2015). Additionally, healthy tissues such as white matter, gray matter, muscle, adipose, or fibro glandular tissue also may be identified from these images. Such data have been used to parameterize models of tumor growth and metastasis (Hormuth et al., 2019b; Neal et al., 2013; Yamamoto et al., 2019).

Figure 4

Representative Images from a Murine Model of Glioma

(A) Contrast-enhanced magnetic resonance image with the brain indicated by the dashed box. (B and C) Transformation of DW-MRI estimates of ADC (B) to cellularity (C). Likewise, (E) and (F) illustrate the transformation of DCE-MRI data (E) to blood volume fraction maps (F). Specifically, (E) shows a representative voxel signal intensity time course from within the tumor, which can then be analyzed to estimate the blood volume fraction. These measurements are then used to parameterize a mathematical model of tumor growth and angiogenesis. Using imaging time points not included in model calibration, the error is then assessed between predicted and measured tumor growth (D and G).

Diffusion-weighted (DW) MRI can provide both anatomical information (e.g., white matter fiber tracks) and physiological measures (e.g., cellularity, necrosis, edema) (Padhani et al., 2009). Briefly, DW-MRI techniques return estimates of the diffusion of water molecules in tissues called the apparent diffusion coefficient (ADC). As the ADC is heavily restricted by cell membranes, it has been used as a surrogate for cellularity within the tumor tissue (Hormuth et al., 2020) (Figure 4B). Measures of cellularity derived from DW-MRI have been employed to initialize and validate models of preclinical glioma (Atuegwu et al., 2012; Hormuth et al., 2017) and clinical breast (Atuegwu et al., 2013; Weis et al., 2015) tumor growth. For example, Weis et al. used DW-MRI estimates of cellularity before and following one cycle of therapy to calibrate a model to predict the pathological complete response of breast tumors to neoadjuvant therapy at the time of surgery (Weis et al., 2015).

18F-flourodeoxyglucose (18F-FDG) is a glucose analog that is metabolically trapped when internalized and phosphorylated by cells, resulting in an accumulation in cells actively consuming glucose. Thus,18F-FDG-PET imaging is used to assess tumor glucose metabolism (Castell and Cook, 2008) and has been employed to assign and calibrate cell proliferation rates in models of tumor growth (Liu et al., 2014). Liu et al. used 18F-FDG to estimate tumor cell proliferation in a mathematical model of pancreatic tumor growth (Liu et al., 2014). Images collected at baseline and at the first follow-up were used to calibrate model parameters, and then the individualized parameters were used to predict tumor growth at the time of the second follow-up, showing a high degree of accuracy between the predicted and observed tumor at the second follow-up for each patient.

Dynamic contrast-enhanced (DCE) MRI can provide measures of healthy and diseased vasculature (Yankeelov and Gore, 2009). In DCE-MRI, several images are collected before, during, and after the injection of a contrast agent, resulting in changes in signal intensity as a function of the concentration of the contrast agent. These signal intensity dynamics can be analyzed with pharmacokinetic models to estimate perfusion, vascular permeability, blood volume, blood flow, and tissue volume fractions. These measures of tumor vasculature (Figure 4E) can be incorporated into models of tumor growth and angiogenesis (Hormuth et al, 2019, 2020). For example, Hormuth et al. used DCE-MRI estimates of blood volume to initialize and calibrate a coupled biophysical model of tumor growth and angiogenesis in a murine model of glioma (Hormuth et al., 2019). They were able to spatially and temporally evolve proliferation and carrying capacity values in response to changes in the local blood volume fraction, achieving a high degree of agreement with less than 10% error between voxel-level predictions and experimental measurements.

Finally, there are several additional MRI and PET techniques not typically collected clinically that measure physiological characteristics that may be valuable for modeling approaches. For example, there exist approaches that measure proliferation (Woolf et al., 2014), pH (Rata et al., 2014), metabolic hallmarks (Sinha and Sinha, 2009), receptor status (Mortimer et al., 2018), and mechanical properties (Glaser et al., 2012). 18F-fluoromisonidazole (18F-FMISO) is a PET tracer (Rajendran and Krohn, 2015) used to visualize tissue hypoxia, a phenomena highly associated with resistance to radiation therapy (Gilkes, 2019). Rockne et al. incorporated measures of hypoxia from 18F-FMISO-PET imaging into patient-specific models of glioblastoma growth to inform the spatial variation in radiation sensitivity for informing more accurate predictions of tumor response to radiation therapy (Rockne et al., 2015). Integration of these MRI and PET-derived measurements with mathematical modeling will provide a better understanding of tumor dynamics in vivo.

Limitations of Using In Vivo Data for Modeling Cancer

One primary limitation of using these data to calibrate mechanism-based mathematical models can be lack of clarity associated with their interpretations. For example, in anatomical imaging, changes in steroid dose or treatment-related reactions can result in transient changes in tumor appearance that are related to neither progression nor regression (Hygino da Cruz et al., 2011). Similarly, although several studies have shown that the ADC correlates inversely with tissue cellularity (Anderson et al., 2000; Barnes et al., 2015; Sugahara et al., 1999), other factors such as cell size, cell permeability, and tissue tortuosity (all of which can change during treatment) also influence ADC values. Similar limitations exist for FDG-PET (false positives can result because of macrophage infiltration [Strauss, 1996]) and for FMISO-PET (false negatives seen within necrotic, non-viable tissue [Hirata et al., 2019]).

A second fundamental limitation of employing in vivo imaging to calibrate mechanism-based models of tumor growth and treatment response is due to the achievable spatial resolution of these techniques. In the clinical setting, anatomical MRI typically produces images with a voxel resolution on the order of 1 mm3, whereas the more quantitative approaches (e.g., DW- and DCE-MRI) may only achieve a resolution of approximately 10 mm3. Similarly, the voxels in PET images are approximately 500 mm3 (Frangioni, 2008). These spatial resolutions limit our ability to quantitatively interrogate tumor features in detail. For example, MRI is unable to capture the microscopic invasion of a glioma (Figure 4) at the tumor margins. Additionally, with a voxel size of 10 mm3, small vessels simply cannot be detected, which will directly affect generated maps of tumor-associated vasculature as well as the characterization of perfusion. Thus, if these parameters are used to approximate (for example) systemic drug delivery, it may result in predictions that over- or underestimate actual tumor response owing to the coarseness of the estimate at which the drug distribution profile can be estimated.

Emerging Applications for Practical Mathematical Modeling

Given a mathematical model that is capable of recapitulating the spatial and temporal development of a tumor undergoing treatment, the model can then be used to optimize the delivery of treatment with the goal of optimizing patient outcomes. This can be achieved by systematically varying (in silico) treatment options and schedules, dosing regimens, or the effect of mechanisms of interest. These in silico tests are a practical and direct method for replacing the trial-and-error approach that is the current paradigm for clinical oncology with a more rigorous approach. Indeed, without a mathematical theory, we are left with trial and error—and given the complex landscape of variables and parameters required to devise an effective treatment plan, there are simply not enough patients or resources to run all the clinical trials that would be required to adequately search that space. The integration of quantitative biological assays with biologically based mathematical modeling is one practical way to move the field from population-based trial and error to patient-specific optimized therapy. Furthermore, optimal control theory has been gaining increasing attention to formulate personalized therapeutic regimens by integrating a biologically based mathematical model of cancer with knowledge of drug pharmacodynamics (Jarrett et al., 2020a). This approach was leveraged by Angaroni et al. to determine optimal personalized regimens of imatinib for CML, aiming to minimize drug exposure and cancer stem cell burden (Angaroni et al., 2020). To this end, the authors developed an ODE model of CML, calibrated with longitudinal measurements of cancer cell burden, as well as personalized pharmacokinetic and pharmacodynamic models of imatinib. Their results suggest that the optimal imatinib therapies can render an improved control of the cancer cell burden with similar or lower drug exposures with respect to standard scheduling for most patients in the study cohort.

Biology-based models enable the transformation of measured data samples into mechanistically interpretable parameters and variables, which help to quantitatively analyze the investigated cancer mechanisms and further guide new research directions (Ayuso et al., 2017; Jarrett et al., 2020b; Johnson et al., 2019, Johnson et al., 2020; Lorenzo et al, 2019b, 2020; Merkher et al., 2020; Morken et al., 2014; Pérez-García et al., 2019; Wang et al., 2006). Moreover, mathematical models parameterized by biological data need not be comprehensive to make to gain insight into underlying tumor dynamics and make valid predictions. For example, studies investigating scaling laws (Pérez-García et al., 2020) and metastasis (Yamamoto et al., 2019) made use of only volumetric measurements from imaging data to model longitudinal tumor growth and response. Additionally, serum prostate-specific antigen (PSA) is a standard clinical measure for patients with prostate cancer undergoing treatment and has been used to parameterize patient-specific models of prostate cancer response. For instance, Brady-Nicholls et al. modeled PSA dynamics for adaptive response prediction of patients with prostate cancer who may be resistant to intermittent androgen deprivation therapy (Brady-Nicholls et al., 2020). Lorenzo et al. proposed an ODE model for the dynamics of PSA after external beam radiotherapy, leveraging clinical PSA data series to calibrate their model and extract a reduced set of patient-specific parameters describing cancer cell proliferation and radiation effects (Lorenzo et al., 2019b). These parameters were then used to define biomarkers to successfully identify tumor recurrence. This study suggests that model-inspired biomarkers may enable the identification of recurring tumors earlier than using standard increasing PSA criteria, thereby extending the window of curability with a secondary treatment.

Bridging the data from cellular to tissue scale models remains a fundamental challenge in biologically based mathematical modeling (Alber et al., 2019; Rockne et al., 2019). Multiscale models seek to simultaneously capture cancer dynamics at both the micro- and macroscale by coupling submodels representing different features of cancer mechanisms at either scale (Deisboeck et al., 2011; Yankeelov et al., 2016a). This strategy has been exploited to model angiogenesis and the motility of invasive cancer cells resulting from epithelial-to-mesenchymal transition (Deisboeck et al., 2011; Vilanova et al., 2017). In both cases, multiscale models track the movement and action of specific cells in response to the local values of macroscopic variables, which usually describe the dynamics of the tumor bulk and key substances (e.g., VEGF, nutrients). In turn, the behavior of these individually tracked cells alter the dynamics of these macroscopic variables. For example, in multiscale cancer models featuring angiogenesis, tumors grow and hypoxic tumor regions progressively vanish as nutrient concentrations increase owing to the sprouting angiogenic vasculature driven by tip endothelial cells, which are activated and consume VEGF released by hypoxic tumor cells. Still, this approach faces two central limitations: the requirement of extensive computational resources and the large number of free parameters whose precise calibration may prove challenging. A primary barrier in the development of these models is finding the balance between mathematical complexity and modeling of relevant processes, particularly processes that can be parameterized with acquirable biological data. Thus, it is of central importance to work with the most parsimonious mathematical model as identified through established model selection schemes (Benzekry et al., 2014; Lima et al., 2016).

Conclusion

We have provided an overview of a number of common experimental assays that quantitatively distinguish tumor characteristics from the molecular to tissue scales and how such data can be directly integrated into biologically based, mathematical models to forecast tumor growth and treatment response in time and space. A summary of the relative strengths and limitations of integrating experimental data and mechanism-based, mathematical modeling of each biological scale is summarized in Table 1. Ultimately, when data and models are successfully combined—in ways that account for the uncertainty in the biological and mathematical strategies—it yields an experimental-mathematical framework with the ability to lower experiment costs, shorten experiment timelines, and maximize both scientific understanding and predictive capabilities. Although the field is still quite young, the early successes strongly suggest that this is an extremely promising direction of investigation that can yield practical insights into both experimental design for cancer biology and patient care for clinical oncology.

Table 1

Strengths and Limitations of Integrating Experimental Data and Mathematical Modeling at Each Biological Scale

Scale	Cellular	Microenvironmental to Tissue Scales	Tissue to Organ Scales
Measures of	• Cancer cell phenotypes and populations	• Interactions between cancer cells and tumor microenvironment	• Physiological characteristics of tumor tissue
Experimental assays	• Whole-genome sequencing • RNA sequencing • Flow cytometry	• Immunofluorescent or immunohistochemical staining and microscopy	• DW-MRI • DCE-MRI • 18F-FDG PET • 18F-FMISO PET
Strengths	• Can gain insight into cell population dynamics under various conditions	• Accounts for role of tumor microenvironment in tumor progression • Offers insight into spatial interaction of tumor microenvironment and cancer cell populations • Cellular-level spatial resolution	• In vivo and noninvasive • 3D, whole tumor imaging • Enables longitudinal measurements of tumor
Limitations	• Limited temporal resolution • Single-cell analysis techniques necessary to account for cellular diversity • Lacking context of tumor microenvironment and other influences in vivo	• Limited temporal resolution • Limited sample sizes	• Limited spatial resolution • Mixed measures of multiple biological phenomena