# Integrating Imaging Data into Predictive Biomathematical and Biophysical Models of Cancer

Authors: Thomas E. Yankeelov
Journal: ISRN Biomathematics
 (2012)

---

## Abstract

While there is a mature literature on biomathematical and biophysical modeling in cancer, many of the existing approaches are not of clinical utility, as they require input data that are extremely difficult to obtain in an intact organism, and/or require a large number of assumptions on the free parameters included in the models. Thus, there has only been very limited application of such models to solve problems of clinical import. More recently, however, there has been increased activity at the interface of quantitative, noninvasive imaging data, and tumor mathematical modeling. In addition to reporting on bulk tumor morphology and volume, emerging imaging techniques can quantitatively report on for example tumor vascularity, glucose metabolism, cell density and proliferation, and hypoxia. In this paper, we first motivate the problem of predicting therapy response by highlighting some (acknowledged) shortcomings in existing methods. We then provide introductions to a number of representative quantitative imaging methods and describe how they are currently (and potentially can be) used to initialize and constrain patient specific mathematical and biophysical models of tumor growth and treatment response, thereby increasing the clinical utility of such approaches. We conclude by highlighting some of the exciting research directions when one integrates quantitative imaging and tumor modeling.


---

## Body

## 1. Introduction

The ability to identify—early in the course of therapy—patients that are not responding to a given therapeutic regimen is extraordinarily important. In addition to limiting patients’ exposure to the toxicities associated with unsuccessful therapies, it would allow patients the opportunity to switch to a potentially more efficacious approach. Unfortunately, existing methods of determining early response are inadequate.  In particular, the current standard-of-care imaging assessment of tumor response to treatment is based on the Response Evaluation Criteria in Solid Tumors (RECIST, [1]). RECIST offers a simplified, practical method for extracting the basic morphological features of imaging data by assessing the overall tumor burden at baseline and comparing that measurement to subsequent measurements obtained during the course of therapy. We now present the salient features of RECIST in order to identify weaknesses that the integration of imaging and mathematical modeling can potentially address.High-resolution images (typically computed tomography (CT) or magnetic resonance imaging (MRI)) are acquired at baseline before treatment as commenced. In these image sets, “target lesions” are determined, and the sum of their longest dimensions is recorded. Additional scans are then acquired during or after therapy and similarly analyzed. The change in the sum of the longest diameters from baseline to the follow-up studies are then calculated and then used to divide treatment response into one of the following four categories [1].(1)
Complete response (CR)—disappearance of all target lesions.(2)
Partial response (PR)—at least a 30% decrease in the sum of the longest diameters of the target lesions, taking as the reference the baseline sum longest diameter.(3)
Progressive disease (PD)—at least a 20% increase in the sum of the longest diameters of the target lesions, taking as the reference the smallest sum longest diameter recorded since the baseline measurements or the appearance of one or more new lesions.(4)
Stable Disease (SD)—neither sufficient shrinkage to qualify for partial response nor sufficient increase to qualify for progressive disease.While almost all clinical trials of solid tumors employ the RECIST criteria, it is well recognized that this approach needs to be significantly improved for a number of reasons. For example, the metric for positive response is based on one-dimensional changes (sum of longest diameters of target lesions) which can be grossly misleading in trying to describe a complex object that is almost certainly changing in all three spatial dimensions. Furthermore, this metric is based on anatomical and morphological changes which are (temporally) downstream manifestations of underlying physiological, cellular, or molecular changes; that is, RECIST based evaluations generally do not indicate whether a tumor is responding until the patient has received several treatment cycles. This limitation is particularly problematic in the era of targeted therapies; in fact, it is a matter for debate as to whether the RECIST criteria are even relevant in assessing the activities of non-cytotoxic anti-cancer drugs if changes in morphology may not be the most appropriate method to assess response. What is needed are methods to characterize those underlying physiological, cellular, and molecular changes as they are highly likely to offerearlierand more specific responseto treatment indices than changes in longest dimensions. One approach to this shortcoming of RECIST is to begin to incorporate some of the more quantitative and specific noninvasive imaging methods into clinical trials and practice; indeed, there is much current effort at accomplishing this goal (see, e.g., [2–4]). Some progress has been made on this front. For example, in the most recent version of RECIST (RECIST 1.1, [5]), fluorodeoxyglucose positron emission tomography (PET) has now been included in the guidelines, though only in the limited role of assessing progressive disease. This represents the first time a nonmorphology-based imaging measurement has been incorporated into a consensus response criteria. It is anticipated that, in coming years, there will be more physiological, cellular, and molecular imaging measurements adopted into response criteria to compliment the morphology-based assessments.Independent of the recent developments in the quantitative imaging of cancer, a mature literature has developed in the biomathematical and biophysical modeling of tumor growth (for reviews of this field see, e.g., [6, 7]). The complexity and range of applicability of these models run the gamut from exponential growth of an avascular tumor [8] to elaborate systems of partial differential equations describing cell motility, tumor invasion, and angiogenesis.  While these studies have provided many interesting and exciting insights into fundamental cancer biology, they are, however, fundamentally limited in their ability to be translated to clinical application. This is because the majority of these models rely on measurements of processes that occur at the cellular and/or molecular level (growth factor gradients, gene expression levels, etc.) and are therefore extremely difficult to measure in an intact living system. Thus, the fundamental limitation of such approaches is that they are driven by parameters that are accessible only through invasive methods (surgery, biopsy, or animal sacrifice). This leads to the necessity to make multiple assumptions about in vivo parameter values based on, for example, idealized in vitro or ex vivo studies, thus such approaches do not generally lend themselves to patient specific modeling. However, there exist several emerging biomedical imaging methods that can provide quantitative information noninvasively, in 3D, and at multiple time points, thereby providing the opportunity for integrating imaging data into mathematical models to predict the growth and response of tumors to therapy. More specifically, measurements can be made at the time of diagnosis and early in the course of treatment, and these data can be used to initialize and constrain models that predict treatment outcomes. In this way, imaging allows models to be initialized and constrained with patient specific data; furthermore, only by parameterizing models by data that are specific to the individual, can such an approach find clinical utility. If realized, this would mark a significant improvement over, for example, RECIST.In this paper we aim to show how biomedical imaging data can be used to initialize and constrain patient specific mathematical models of tumor growth to predict tumor status at later time points. To accomplish this task, we first provide brief introductions of a selection of emerging, quantitative, andin vivo imaging measurements which are appropriate for incorporation into predictive biomathematical and biophysical models of tumor growth. We then follow that with a section reviewing the literature where imaging techniques have been successfully incorporated into tumor modeling. The final section includes a discussion of future opportunities at the interface of imaging and modeling and describes three key areas (two theoretical and one experimental) that should be primary goals for investigators working in this exciting and emerging field.


## 2. Data Available from Biomedical Imaging

We limit our discussion of quantitative imaging metrics to the (currently) most widely available MRI and PET measures. Additionally, these imaging measures are those that have already shown encouraging data in their ability to predict treatment response and are, therefore, natural candidates to integrate into biomathematical and biophysical models of tumor growth and treatment response.


### 2.1. Magnetic Resonance Imaging

#### 2.1.1. Dynamic Contrast Enhanced MRI (DCE-MRI)

DCE-MRI involves the serial acquisition of MR images of a tissue of interest (e.g., a tumor) before and after an intravenous injection of paramagnetic contrast agent [9]. As the contrast agent enters into the tissue under investigation, the T1 and T2 values of tissue water decrease to an extent that is determined by the concentration of the contrast; studies that exploit changes in T1 are termed DCE-MRI, while those relying on T2 changes are termed dynamic susceptibility contrast MRI (DSC-MRI; [10]). By considering a set of images acquired before and after a contrast infusion, a region of interest (ROI) or individual voxels will display a characteristic signal intensity time course, which can be related to the concentration of contrast agent. Fitting the DCE-MRI data to an appropriate pharmacokinetic model allows for extraction of physiological parameters related to, for example, tissue perfusion and microvascular vessel wall permeability (denoted by the parameter Ktrans), extracellular volume fraction (denoted by ve), and the plasma volume fraction (denoted by vp). It has been shown that these parameters can change significantly during the course of therapy and be used to predict patients who will have a positive response from those that do not; see, for example, [4, 11–14]. For a detailed review of DCE-MRI and its applications in cancer, the interested reader is referred to [9].


#### 2.1.2. Diffusion Weighted MRI (DW-MRI)

The microscopic thermally-induced behavior of molecules moving in a random pattern is referred to as self-diffusion or Brownian motion [15]. The rate of diffusion in tissues is lower than in free solution and is described by an apparent diffusion coefficient (ADC), which largely depends on the number, permeability and separation of barriers that a diffusing water molecule encounters. MRI methods have been developed to map the ADC, and the variations in ADC have been shown to correlate inversely with tissue cellularity [16]. In regions of high cellularity (e.g., a tumor) there is an increase in cellular membranes thereby providing an increase in the number of barriers encountered by a diffusing water molecule; this results in a decreased ADC. Conversely, when effective cytotoxic therapy reduces the overall cellularity and thereby decreases the number of barriers a diffusing water molecule encounters, there is an increase in the ADC. It has been shown that exposure of tumors to chemotherapy and/or radiotherapy consistently leads to measurable increases in water diffusion in cases of favorable treatment response [17–19]. For a detailed review of DW-MRI and its applications in cancer, the interested reader is referred to [20].


#### 2.1.3. Diffusion Tensor Imaging (DTI)

Diffusion tensor imaging (DTI) is a method to noninvasively characterize the structural connectivity between brain cortical regions. DTI provides, for each voxel, a tensor matrix that describes the constraints on local Brownian motion of water molecules. Originally proposed to assess tissue properties such as diffusion anisotropy [21], DTI has recently evolved to become the primary modality for mapping neuronal fiber tracts in vivo, a technique referred to as DTI tractography. DTI tractography draws upon the principle that the dominant direction of water diffusion coincides with the local tangential direction of fibrous tissue, and integration of the principal diffusion directions permits entire fiber tracts to be delineated [22]. Based on these approaches, a number of streamline-like fiber tracking algorithms have been developed to generate the pathways of fiber connections [22–25]. More recently, a number of investigators have examined both how these fiber tracks affect tumor growth and how the growing tumor affects the fiber tracks [26–28]. For a detailed review of DTI and its applications in cancer, the interested reader is referred to [29].


### 2.2. Positron Emission Tomography

#### 2.2.1. Fluorodeoxyglucose PET (FDG-PET)

Fluorodeoxyglucose is a glucose analogue that accumulates in areas of increased glycolytic activity which is a near universal property of cancer. It is well known that the activity of cell surface glucose transporters GLUT-1 and GLUT-3 and the intracellular enzyme hexokinase are upregulated in malignant cancer cells. The transporters function is to transport FDG into tumor cells where it is phosphorylated by hexokinase and trapped (and hence accumulated) in malignant cells. FDG-PET has been used extensively in both the preclinical and clinical settings to study treatment response [30, 31]. (It is the far and away the most commonly used PET tracer in clinical oncology.) In particular, reductions in the delivery and retention of PET have been shown to correlate to and, in some cases, predict therapeutic response [31–34]. For a detailed review of FDG-PET and its applications in cancer, the interested reader is referred to [35].


#### 2.2.2. Fluorodeoxythymidine PET (PET-FLT)

Fluorodeoxythymidine (FLT) was developed as a surrogate marker of cellular proliferation [34]. Once transported to the tumor intracellular space FLT is phosphorylated by thymidine kinase-1 (TK1), which is known to have a large increase in activity during DNA synthase phase of the cell cycle (S phase). Since FLT monophosphate is not incorporated into DNA and is impermeable to the cell membrane, it is metabolically trapped intracellularly. Thus, rapidly proliferating (tumor) cells can present an increased retention of FLT when compared to background (healthy) tissue. As with FDG, there has been much recent interest in FLT in both preclinical and clinical studies of treatment response [34, 36–39]. The general working hypothesis for FLT-PET studies is that treatments designed to reduce cellular proliferation will show decreased FLT uptake, and this will be evident in responders, whereas nonresponders will not display a decrease in FLT accumulation.  Indeed, many recent studies have displayed just such a relationship. For a detailed review of FLT-PET and its applications in cancer, the interested reader is referred to [40].


#### 2.2.3. Hypoxia Imaging via PET

It is well known that hypoxia in malignant tumors can affect the outcome of anticancer treatments, as malignant tumors are more resistant to chemotherapy and irradiative therapy because of their lack of oxygen, which is a potent radiosensitizer [41]. Thus, there has been much effort in developing imaging methods for characterizing tissue oxygen status [42]; in particular, there have been several approaches employing derivatives of misonidazole. Upon entering a cell misonidazole is reduced and retained in hypoxic cells but freely exits the cell under normoxic conditions. Using 14C-labeled misonidazole,  Chapman was the first to demonstrate the potential of radiolabeled nitroimidazoles for the detection of hypoxia [43]. The clinical detection of hypoxia became available with the introduction of 18F-Fluoromisonidazole [44]. Several studies have since demonstrated the clinical prognostic value of FMISO [45, 46]. FMISO is also being used to assess hypoxia and oxygen-enhancing therapeutics in experimental cancer models [47, 48]. For a detailed review of FMISO-PET and its applications in cancer, the interested reader is referred to [49].Another common PET-based method for studying at hypoxia is based on copper(II)-diacetyl-bis(N4-methylthiosemicarbazone), Cu–ATSM, which can be labeled with a positron emitting isotope of copper (xCu, where x=60,61,62, or 64) [50].  Owing to its high membrane permeability and low redox potential, xCu-ATSM tracers can access the intracellular space and are stable there under normoxia. However, when entering hypoxic tissue the xCu(II) in xCu-ATSM is reduced to xCu(I) and is released from the ATSM chelate, thereby becoming trapped intracellularly and, therefore, xCu-ATSM is preferentially accumulated in hypoxic tissues [51, 52]. These agents have been used to study, for example, the ability to predict treatment response [53, 54] and determine prognosis [55, 56]. For a detailed review of xCu-ATSM-PET and its applications in cancer, the interested reader is referred to [57].


### 2.3. Multimodal Imaging Data

As should be apparent from the above discussions, many of the imaging modalities are complimentary in nature, and this realization has spurred the development of multimodal assessments of cancer. For example, combining the measurements available from DCE-MRI and FMISO-PET can provide insight into the relationship between vascular status and tissue oxygenation. Another obvious paring is DW-MRI and FLT-PET to explore the (temporal and spatial) relationships between cell proliferation and cellularity. To date the overwhelming majority of multimodal studies have been correlative in nature and there have been only minimal efforts at integrating such data into an appropriate biomathematical or biophysical model. (We review such investigations below in Section3.2) However, as multimodality imaging of cancer continues to develop [58–60], and clinical PET-MRI scanners come online [61–63], the ability to acquire quantitative, multimodality imaging data becomes increasingly accessible, and consequently, the ability to initialize and constrain mathematical models of tumor growth becomes increasingly realistic. This is an area for which there is great opportunity for investigations of real clinical impact.


## 2.1. Magnetic Resonance Imaging

### 2.1.1. Dynamic Contrast Enhanced MRI (DCE-MRI)

DCE-MRI involves the serial acquisition of MR images of a tissue of interest (e.g., a tumor) before and after an intravenous injection of paramagnetic contrast agent [9]. As the contrast agent enters into the tissue under investigation, the T1 and T2 values of tissue water decrease to an extent that is determined by the concentration of the contrast; studies that exploit changes in T1 are termed DCE-MRI, while those relying on T2 changes are termed dynamic susceptibility contrast MRI (DSC-MRI; [10]). By considering a set of images acquired before and after a contrast infusion, a region of interest (ROI) or individual voxels will display a characteristic signal intensity time course, which can be related to the concentration of contrast agent. Fitting the DCE-MRI data to an appropriate pharmacokinetic model allows for extraction of physiological parameters related to, for example, tissue perfusion and microvascular vessel wall permeability (denoted by the parameter Ktrans), extracellular volume fraction (denoted by ve), and the plasma volume fraction (denoted by vp). It has been shown that these parameters can change significantly during the course of therapy and be used to predict patients who will have a positive response from those that do not; see, for example, [4, 11–14]. For a detailed review of DCE-MRI and its applications in cancer, the interested reader is referred to [9].


### 2.1.2. Diffusion Weighted MRI (DW-MRI)

The microscopic thermally-induced behavior of molecules moving in a random pattern is referred to as self-diffusion or Brownian motion [15]. The rate of diffusion in tissues is lower than in free solution and is described by an apparent diffusion coefficient (ADC), which largely depends on the number, permeability and separation of barriers that a diffusing water molecule encounters. MRI methods have been developed to map the ADC, and the variations in ADC have been shown to correlate inversely with tissue cellularity [16]. In regions of high cellularity (e.g., a tumor) there is an increase in cellular membranes thereby providing an increase in the number of barriers encountered by a diffusing water molecule; this results in a decreased ADC. Conversely, when effective cytotoxic therapy reduces the overall cellularity and thereby decreases the number of barriers a diffusing water molecule encounters, there is an increase in the ADC. It has been shown that exposure of tumors to chemotherapy and/or radiotherapy consistently leads to measurable increases in water diffusion in cases of favorable treatment response [17–19]. For a detailed review of DW-MRI and its applications in cancer, the interested reader is referred to [20].


### 2.1.3. Diffusion Tensor Imaging (DTI)

Diffusion tensor imaging (DTI) is a method to noninvasively characterize the structural connectivity between brain cortical regions. DTI provides, for each voxel, a tensor matrix that describes the constraints on local Brownian motion of water molecules. Originally proposed to assess tissue properties such as diffusion anisotropy [21], DTI has recently evolved to become the primary modality for mapping neuronal fiber tracts in vivo, a technique referred to as DTI tractography. DTI tractography draws upon the principle that the dominant direction of water diffusion coincides with the local tangential direction of fibrous tissue, and integration of the principal diffusion directions permits entire fiber tracts to be delineated [22]. Based on these approaches, a number of streamline-like fiber tracking algorithms have been developed to generate the pathways of fiber connections [22–25]. More recently, a number of investigators have examined both how these fiber tracks affect tumor growth and how the growing tumor affects the fiber tracks [26–28]. For a detailed review of DTI and its applications in cancer, the interested reader is referred to [29].


## 2.1.1. Dynamic Contrast Enhanced MRI (DCE-MRI)

DCE-MRI involves the serial acquisition of MR images of a tissue of interest (e.g., a tumor) before and after an intravenous injection of paramagnetic contrast agent [9]. As the contrast agent enters into the tissue under investigation, the T1 and T2 values of tissue water decrease to an extent that is determined by the concentration of the contrast; studies that exploit changes in T1 are termed DCE-MRI, while those relying on T2 changes are termed dynamic susceptibility contrast MRI (DSC-MRI; [10]). By considering a set of images acquired before and after a contrast infusion, a region of interest (ROI) or individual voxels will display a characteristic signal intensity time course, which can be related to the concentration of contrast agent. Fitting the DCE-MRI data to an appropriate pharmacokinetic model allows for extraction of physiological parameters related to, for example, tissue perfusion and microvascular vessel wall permeability (denoted by the parameter Ktrans), extracellular volume fraction (denoted by ve), and the plasma volume fraction (denoted by vp). It has been shown that these parameters can change significantly during the course of therapy and be used to predict patients who will have a positive response from those that do not; see, for example, [4, 11–14]. For a detailed review of DCE-MRI and its applications in cancer, the interested reader is referred to [9].


## 2.1.2. Diffusion Weighted MRI (DW-MRI)

The microscopic thermally-induced behavior of molecules moving in a random pattern is referred to as self-diffusion or Brownian motion [15]. The rate of diffusion in tissues is lower than in free solution and is described by an apparent diffusion coefficient (ADC), which largely depends on the number, permeability and separation of barriers that a diffusing water molecule encounters. MRI methods have been developed to map the ADC, and the variations in ADC have been shown to correlate inversely with tissue cellularity [16]. In regions of high cellularity (e.g., a tumor) there is an increase in cellular membranes thereby providing an increase in the number of barriers encountered by a diffusing water molecule; this results in a decreased ADC. Conversely, when effective cytotoxic therapy reduces the overall cellularity and thereby decreases the number of barriers a diffusing water molecule encounters, there is an increase in the ADC. It has been shown that exposure of tumors to chemotherapy and/or radiotherapy consistently leads to measurable increases in water diffusion in cases of favorable treatment response [17–19]. For a detailed review of DW-MRI and its applications in cancer, the interested reader is referred to [20].


## 2.1.3. Diffusion Tensor Imaging (DTI)

Diffusion tensor imaging (DTI) is a method to noninvasively characterize the structural connectivity between brain cortical regions. DTI provides, for each voxel, a tensor matrix that describes the constraints on local Brownian motion of water molecules. Originally proposed to assess tissue properties such as diffusion anisotropy [21], DTI has recently evolved to become the primary modality for mapping neuronal fiber tracts in vivo, a technique referred to as DTI tractography. DTI tractography draws upon the principle that the dominant direction of water diffusion coincides with the local tangential direction of fibrous tissue, and integration of the principal diffusion directions permits entire fiber tracts to be delineated [22]. Based on these approaches, a number of streamline-like fiber tracking algorithms have been developed to generate the pathways of fiber connections [22–25]. More recently, a number of investigators have examined both how these fiber tracks affect tumor growth and how the growing tumor affects the fiber tracks [26–28]. For a detailed review of DTI and its applications in cancer, the interested reader is referred to [29].


## 2.2. Positron Emission Tomography

### 2.2.1. Fluorodeoxyglucose PET (FDG-PET)

Fluorodeoxyglucose is a glucose analogue that accumulates in areas of increased glycolytic activity which is a near universal property of cancer. It is well known that the activity of cell surface glucose transporters GLUT-1 and GLUT-3 and the intracellular enzyme hexokinase are upregulated in malignant cancer cells. The transporters function is to transport FDG into tumor cells where it is phosphorylated by hexokinase and trapped (and hence accumulated) in malignant cells. FDG-PET has been used extensively in both the preclinical and clinical settings to study treatment response [30, 31]. (It is the far and away the most commonly used PET tracer in clinical oncology.) In particular, reductions in the delivery and retention of PET have been shown to correlate to and, in some cases, predict therapeutic response [31–34]. For a detailed review of FDG-PET and its applications in cancer, the interested reader is referred to [35].


### 2.2.2. Fluorodeoxythymidine PET (PET-FLT)

Fluorodeoxythymidine (FLT) was developed as a surrogate marker of cellular proliferation [34]. Once transported to the tumor intracellular space FLT is phosphorylated by thymidine kinase-1 (TK1), which is known to have a large increase in activity during DNA synthase phase of the cell cycle (S phase). Since FLT monophosphate is not incorporated into DNA and is impermeable to the cell membrane, it is metabolically trapped intracellularly. Thus, rapidly proliferating (tumor) cells can present an increased retention of FLT when compared to background (healthy) tissue. As with FDG, there has been much recent interest in FLT in both preclinical and clinical studies of treatment response [34, 36–39]. The general working hypothesis for FLT-PET studies is that treatments designed to reduce cellular proliferation will show decreased FLT uptake, and this will be evident in responders, whereas nonresponders will not display a decrease in FLT accumulation.  Indeed, many recent studies have displayed just such a relationship. For a detailed review of FLT-PET and its applications in cancer, the interested reader is referred to [40].


### 2.2.3. Hypoxia Imaging via PET

It is well known that hypoxia in malignant tumors can affect the outcome of anticancer treatments, as malignant tumors are more resistant to chemotherapy and irradiative therapy because of their lack of oxygen, which is a potent radiosensitizer [41]. Thus, there has been much effort in developing imaging methods for characterizing tissue oxygen status [42]; in particular, there have been several approaches employing derivatives of misonidazole. Upon entering a cell misonidazole is reduced and retained in hypoxic cells but freely exits the cell under normoxic conditions. Using 14C-labeled misonidazole,  Chapman was the first to demonstrate the potential of radiolabeled nitroimidazoles for the detection of hypoxia [43]. The clinical detection of hypoxia became available with the introduction of 18F-Fluoromisonidazole [44]. Several studies have since demonstrated the clinical prognostic value of FMISO [45, 46]. FMISO is also being used to assess hypoxia and oxygen-enhancing therapeutics in experimental cancer models [47, 48]. For a detailed review of FMISO-PET and its applications in cancer, the interested reader is referred to [49].Another common PET-based method for studying at hypoxia is based on copper(II)-diacetyl-bis(N4-methylthiosemicarbazone), Cu–ATSM, which can be labeled with a positron emitting isotope of copper (xCu, where x=60,61,62, or 64) [50].  Owing to its high membrane permeability and low redox potential, xCu-ATSM tracers can access the intracellular space and are stable there under normoxia. However, when entering hypoxic tissue the xCu(II) in xCu-ATSM is reduced to xCu(I) and is released from the ATSM chelate, thereby becoming trapped intracellularly and, therefore, xCu-ATSM is preferentially accumulated in hypoxic tissues [51, 52]. These agents have been used to study, for example, the ability to predict treatment response [53, 54] and determine prognosis [55, 56]. For a detailed review of xCu-ATSM-PET and its applications in cancer, the interested reader is referred to [57].


## 2.2.1. Fluorodeoxyglucose PET (FDG-PET)

Fluorodeoxyglucose is a glucose analogue that accumulates in areas of increased glycolytic activity which is a near universal property of cancer. It is well known that the activity of cell surface glucose transporters GLUT-1 and GLUT-3 and the intracellular enzyme hexokinase are upregulated in malignant cancer cells. The transporters function is to transport FDG into tumor cells where it is phosphorylated by hexokinase and trapped (and hence accumulated) in malignant cells. FDG-PET has been used extensively in both the preclinical and clinical settings to study treatment response [30, 31]. (It is the far and away the most commonly used PET tracer in clinical oncology.) In particular, reductions in the delivery and retention of PET have been shown to correlate to and, in some cases, predict therapeutic response [31–34]. For a detailed review of FDG-PET and its applications in cancer, the interested reader is referred to [35].


## 2.2.2. Fluorodeoxythymidine PET (PET-FLT)

Fluorodeoxythymidine (FLT) was developed as a surrogate marker of cellular proliferation [34]. Once transported to the tumor intracellular space FLT is phosphorylated by thymidine kinase-1 (TK1), which is known to have a large increase in activity during DNA synthase phase of the cell cycle (S phase). Since FLT monophosphate is not incorporated into DNA and is impermeable to the cell membrane, it is metabolically trapped intracellularly. Thus, rapidly proliferating (tumor) cells can present an increased retention of FLT when compared to background (healthy) tissue. As with FDG, there has been much recent interest in FLT in both preclinical and clinical studies of treatment response [34, 36–39]. The general working hypothesis for FLT-PET studies is that treatments designed to reduce cellular proliferation will show decreased FLT uptake, and this will be evident in responders, whereas nonresponders will not display a decrease in FLT accumulation.  Indeed, many recent studies have displayed just such a relationship. For a detailed review of FLT-PET and its applications in cancer, the interested reader is referred to [40].


## 2.2.3. Hypoxia Imaging via PET

It is well known that hypoxia in malignant tumors can affect the outcome of anticancer treatments, as malignant tumors are more resistant to chemotherapy and irradiative therapy because of their lack of oxygen, which is a potent radiosensitizer [41]. Thus, there has been much effort in developing imaging methods for characterizing tissue oxygen status [42]; in particular, there have been several approaches employing derivatives of misonidazole. Upon entering a cell misonidazole is reduced and retained in hypoxic cells but freely exits the cell under normoxic conditions. Using 14C-labeled misonidazole,  Chapman was the first to demonstrate the potential of radiolabeled nitroimidazoles for the detection of hypoxia [43]. The clinical detection of hypoxia became available with the introduction of 18F-Fluoromisonidazole [44]. Several studies have since demonstrated the clinical prognostic value of FMISO [45, 46]. FMISO is also being used to assess hypoxia and oxygen-enhancing therapeutics in experimental cancer models [47, 48]. For a detailed review of FMISO-PET and its applications in cancer, the interested reader is referred to [49].Another common PET-based method for studying at hypoxia is based on copper(II)-diacetyl-bis(N4-methylthiosemicarbazone), Cu–ATSM, which can be labeled with a positron emitting isotope of copper (xCu, where x=60,61,62, or 64) [50].  Owing to its high membrane permeability and low redox potential, xCu-ATSM tracers can access the intracellular space and are stable there under normoxia. However, when entering hypoxic tissue the xCu(II) in xCu-ATSM is reduced to xCu(I) and is released from the ATSM chelate, thereby becoming trapped intracellularly and, therefore, xCu-ATSM is preferentially accumulated in hypoxic tissues [51, 52]. These agents have been used to study, for example, the ability to predict treatment response [53, 54] and determine prognosis [55, 56]. For a detailed review of xCu-ATSM-PET and its applications in cancer, the interested reader is referred to [57].


## 2.3. Multimodal Imaging Data

As should be apparent from the above discussions, many of the imaging modalities are complimentary in nature, and this realization has spurred the development of multimodal assessments of cancer. For example, combining the measurements available from DCE-MRI and FMISO-PET can provide insight into the relationship between vascular status and tissue oxygenation. Another obvious paring is DW-MRI and FLT-PET to explore the (temporal and spatial) relationships between cell proliferation and cellularity. To date the overwhelming majority of multimodal studies have been correlative in nature and there have been only minimal efforts at integrating such data into an appropriate biomathematical or biophysical model. (We review such investigations below in Section3.2) However, as multimodality imaging of cancer continues to develop [58–60], and clinical PET-MRI scanners come online [61–63], the ability to acquire quantitative, multimodality imaging data becomes increasingly accessible, and consequently, the ability to initialize and constrain mathematical models of tumor growth becomes increasingly realistic. This is an area for which there is great opportunity for investigations of real clinical impact.


## 3. Current Efforts at Integrating Imaging and Mathematical Models

We divide our review of the literature on integrating modeling and imaging into subsections in which a single imaging technique is employed, and a section where multiple imaging techniques are employed. We note that the order of the manuscripts reviewed below was determined by their publication date.


### 3.1. Modeling Studies Employing One Imaging Technique

As it is well known that glioma cells preferentially migrate along white matter fibers of the brain, it is natural to incorporate the structurala priori information on brain architecture provided by diffusion tensor imaging to guide how tumor cells will move from their source point.  This is the approach taken by Jbabdi et al. [64]. The authors began with a reaction-diffusion model of how tumor cells proliferate and migrate as follows:
(1)∂c∂t=D∇2c(x,t)+ρc(x,t),
where c(x,t) is the concentration of tumor cells at position x and time t, D is the diffusion coefficient associated with random tumor cell motion, and ρ is the proliferation rat of the tumor cells.  The first term on the right hand side of (1) is the diffusion component, while the “reaction” portion is characterized by the exponential growth term (the second term on the right hand side). As it is written, (1) does not account for the heterogeneous structure of the brain, and one way to incorporate such knowledge is to make D an explicit function of space dependent upon information provided by the DTI data. In this case, we have
(2)∂c∂t=∇·(D(x)∇c(x,t))+ρc(x,t),
where D(x) is the diffusion tensor that describes anisotropic tumor cell diffusion and therefore takes into account both location and direction as reported by the DTI data. By combining an average value for ρ (obtained from previous investigations [65, 66]) with appropriate initial (i.e., c(x,0)=co(x)) and boundary (D(x)∇c(x,t)t)·n=0 for all x on the sulcal and ventricular boundary, where n is the surface normal) conditions, the model can then be run forward to predict tumor shape and distribution which can then be compared to experimental results. In their initial effort, the authors used DTI data from a healthy volunteer, as it was not available on the patient data sets. The authors stated that their goal was to simulate tumor growth in a healthy brain to best fit the tumor size and shape that was experimentally observed in patients with brain tumors low-grade tumors. By “seeding” a tumor in the healthy brain data set and “growing” it forward in time, they were able to compare to their predicted tumors to those observed in real patients. The results indicated that when the tumors were allowed to grown anisotropically, as guided by the DTI data, the spatial extension and distribution of the tumor more closely resembled what was observed in the patient, then when the tumor cells were allowed to migrate isotropically without guidance from the DTI data. This preliminary study showed the importance of incorporating structural information available from imaging data in order to account for more realistic tumor cell growth.A similar approach to linking tumor growth and restricted movement was investigated by Clatz et al., who simulated the growth of glioblastoma multiforma (GBM) by incorporating the mechanical restrictions (based on classical continuum mechanics) presented by certain brain structures [67]. Their overall goal was to build a patient-specific simulator of GBM which includes brain deformation induced by the invading tumor. To conceptualize this, the tumor is broken down into two components: (1) the expansive, noninfiltrative component that produces the increase in tumor volume and is responsible for the mass effect and (2) the component associated with tumor diffusion that is responsible for the infiltrative component. The components are coupled in the sense that the mass effect is directly related to tumor cell density, but the tumor cell density is not influenced by the mass effect. To implement this approach, the authors first segment the images into an expansive component and a diffusion component, and the images are then registered to a brain (anatomical) atlas. Then a finite-element mesh is built based on the segmentations, and tissue properties are assigned to each mesh point including the tumor cell density based on the two segmented components. Based on this initialization, the reaction diffusion equation is run forward to generate a virtual GBM whose results are then compared to the patient data; the simulated tumor isodensity contours and brain deformation are compared to what is actually measured on MR images by registering the predicted tumor to that acquired approximately six months after the initial measurements. The algorithm was applied to two patients, and the predicted tumor sizes and shapes showed good qualitative agreement (no quantitative comparison was performed) to what was actually measured on anatomical T2-weighted MR data. These initial results were very encouraging and, importantly, showed the feasibility of using limited clinical imaging data to initialize a patient specific model of tumor growth. The authors concluded by noting that this success justified its further evaluation which then did in a follow-up paper.Bondiau et al. [68] combined the main features of the previous two paragraphs—DTI measurements of white matter tracks and mechanically restricted tumor growth—to model GBM in an excellent paper in 2008.  The DTI data is used to estimate the preferred direction of migration for the diffusion component of the reaction diffusion (2) which is then coupled to a mechanical component.  The authors again used classical linear elasticity theory to describe how the brain reacts to the tumor. In particular, they assume a linear relation between stress and strain as follows:
(3)σ=kε,
where σ is the stress tensor, k is the elasticity of the brain, and ε is the linearized Lagrange strain tensor. Since the deformation is assumed to be very slow the authors use the static equilibrium equation which states that the stress induced by tumor growth is balanced by the external force, Fext, due to tissue deformation:
(4)∇·σ+Fext=0.
The model parameters can then be tuned to match that of an MRI data obtained from a patient; that is, the different brain structures have different mechanical properties, and these are assigned in a straightforward fashion by image segmentation and literature values [69]. These assignments, together with the model can then be used to simulate tumor status at a future time point. The simulation is then compared to what is actually observed at that future time point. In this contribution, the authors were able to perform a set of quantitative comparisons including the average and odd ratio of the differences between the simulated and observed images, as well as displacement of reference points. The authors were able to show good correlation between the simulated and measured tumor masses; in particular, they stated that the inclusion of DTI data to “guide” tumor migration improves the “realism” of the diffusion component of (2). More recently, the authors have extended their approach to compare predicted growth patterns with irradiation margins in which they showed that microscopic invasion occurs past the irradiation margins and that some healthy tissue was within the radiation margins [68]. (We return to this important topic below.) Clearly, this is an exciting direction of investigation.Also using MRI data to drive the reaction-diffusion equation in GBM patients, Wang et al. employed serial pretreatment MRI data to generate patient specific rates of net glioma cell proliferation and invasion in 32 patients [70].  The results in this paper were based on previous efforts by this group in which an imaging-based method was developed to compute the parameters ρ and D from serial pre-treatment MRI data [71]. Briefly, this is done by using the serial T1-weighted (postcontrast) images in which the volumes of the tumor mass visible on these image sets were computed, and then the radius of a sphere of an equivalent volume is computed. The difference between the measured radii, along with the time interval between subsequent scans, allows for the calculation of the velocity of radial expansion, v. This velocity is then substituted into Fisher’s approximation [72]: (5)v=2Dρ.
In addition to the T1-weighted images, T2-weighted images are collected which allows for calculating D/ρ (the so-called, “index of invisibility” [71]. Given this system of two equations in two unknowns it is straight forward to compute D and ρ. Using this approach, the authors showed that their estimates of proliferation and invasion proved to be significantly related to prognosis; in particular, the ratio ρ/D was a significant predictor of survival (P=0.01). Additionally, the authors compute a new measure, termed the therapeutic response index (TRI), which is the ratio of the actual patient survival time to the predicted survival time achieved without therapeutic intervention. Thus, the TRI provides a quantitative assessment of the efficacy of a particular treatment paradigm. The authors concluded that their approach of developing patient-specific virtual controls for application in both treatment design as well as monitoring the effects of that treatment.Rockne et al. extended (2) by adding a term describing the linear-quadratic radiobiological model to describe the effects of radiation therapy [73, 74]. This model considers the surviving fraction, S, of cells as
(6)S=exp[-(αDT+βDT)],
where DT is total dose administered, and α and β are constants from experiment. In practice, DT is fractionated over n fractions, and the fractionated dose, d, is a function of both space and time, that is, d≡d(x,t). Then the biologically effective dose can be shown to be [74]
(7)BED=nd(1+dα/β).
With these assignments, equation (2) is thus updated to
(8)∂c∂t=∇·(D(x)∇c(x,t))+ρc(x,t)-R(α,d(x,t),t)c(x,t),
where
(9)R(α,d(x,t))={0,fortnotduringtherapy1-S(α,d(x,t)),fortduringtherapy}.
Using this model, the investigators could simulate how a tumor treated with radiation therapy would compare to an untreated tumor, as well as estimate a “typical response window.” The group then moved to apply this methodology to patient data. Using pretreatment MRIs to estimate patient specific ρ and D parameters as described above, Rockne et al. made predictions about what the response to radiation would be and correlated these results to clinical data [75]. In particular, they showed the proliferation rate was highly correlated to the radiation response parameter (r=0.89;P=0.0007) which allowed for the construction of a predictive relationship that was then tested in a “leave-one-out” approach. This result let the authors to conclude that their model can generate in silico tumors, with characteristics similar to the in vivo (i.e., patient) tumors and that this framework can then be used to predict individual response.Our own group has recently contributed a simple approach to integrating imaging and modeling, whereby serial DW-MRI data obtained before and early after the initiation of therapy are used to populate the logistic growth model [76–79] and predict tumor cellularity at the conclusion of therapy. The goal is to use the measurement of the ADC to first convert to estimates of cell number [16], and then use changes in the cell number to convert to a proliferation value. More specifically, we convert measurements of ADC to tumor cell number via (10) as follows:
(10)N(t)=θ(ADCw-ADC(t)ADCw-ADCmin),
where ADCw and ADCmin are the ADC of free water and the voxel with the minimum ADC, respectively; the voxel with the minimum ADC (ADCmin) contains the maximum number of cells that can fit into that voxel (i.e., the carrying capacity, θ). The voxel with the ADC of free water (ADCw) is taken to have zero tumor cells. Thus, given ADC data it is straightforward to convert to tumor cell number. If this approach is applied at multiple time points, it is possible to computer the proliferation rate, k, of the tumor on a region of interest or voxel level by fitting the N(t) data to, say, (11):
(11)N(t)=θN(t1)N(t1)+(θ-N(t1))e-kt.
(We note that the parameter k is actually a bulk (averaged) measure of both proliferation and death and can therefore take on either positive or negative values.) Our first effort [77] using this simple approach was in a small animal model of cancer in which 13 rats were given 9L gliosarcoma brain tumors; eight of the rats were treated with the chemotherapeutic drug 1,3-bis(2-chloroethyl)-1-nitrosourea, and five rats were untreated controls. All animals underwent DW-MRI immediately before, one day and three days after treatment. The ADC data from the first imaging time points were used to estimate k as just described which was then used to predict the number of tumor cells at data three, and this value was correlated with the day five experimental data.  On a region of interest basis, the Pearson’s correlation coefficients between the predicted and measured tumor cellularity was 0.88. Thus, this simple approach has some ability to capture the gross (i.e., at the ROI level) cellularity changes.We also applied this approach to breast cancer data acquired during neoadjuvant chemotherapy. Six patients received DW-MRI before, after one cycle, and after all cycles of neoadjuvant chemotherapy. Again, the proliferation rates were estimated using the ADC data from the first two time points and then used with (10) and (11) to predict cellularity after therapy. When comparing the predicted to the experimentally observed data, the Pearson’s correlation coefficient for the ROI analysis yielded 0.95 (P=0.004), and, encouragingly, after applying a conservative 3 × 3 mean filter to the ADC data, the voxel-by-voxel analysis yielded a Pearsons correlation coefficient of 0.70±0.10 (P<0.05). More recently, we have had some preliminary success in using this approach to separate responders from nonresponders [80].Another effort that has incorporated serially acquired DW-MRI data into a mathematical model was contributed by Ellingson et al. [81]. As there is an established negative correlation between ADC and cellularity [16], the ADC can be directly substituted into (2) to capture
(12)ddtADC=∇·(D∇ADC)+ρADC,
and, subsequently,
(13)ddtADC=D∇2ADC+∇D·∇ADC+ρADC.
Algebraically solving (13) for ρ at time point n-1 allows for an estimate of ρ; that is,
(14)ρn-1=(1ADCn-1)×(ddtADCn-1-D∇2ADCn-1-∇D·∇ADCn-1),
where d/dt of ADCn-1 defines the rate of change of ADC from the previous time point n-1 to the first time point n-2. After ρ is fixed, it is substituted back into (12) to yield a PDE that can be solved to estimate D. The authors applied this approach to obtain estimates of ρ, D, and v in 53 patients. In nine of the patients, magnetic resonance spectroscopy was performed to measured Cho/NAA (i.e., choline to N-acetylaspartate ratio); the Pearson’s correlation coefficient between ρ and Cho/NAA was 0.97 (P<0.0001). Furthermore, using a “hot-spot analysis” (based on abnormal anatomical MRI regions), the authors were able to use ρ and D to statistically separate (World Health Organization, WHO) grade II from grades III and IV. The authors then extended this effort to predict progression free and overall survival (OS) in 26 recurrent glioblastoma patients treated with the antiangiogenic drug bevacizumab [82]. Again, the authors were able to find strong correlations between ρ and progression-free survival (P<0.0001) and overall survival (P=0.0093).


### 3.2.  Modeling Studies Employing Multiple Imaging Techniques

One of the first studies to include multimodality imaging in a mathematical model of tumor growth was put forth by Titz and Jeraj in 2008 [83]. In this study, the authors employed FDG-PET/CT data to investigate the effects of oxygenate status on the response to radiation therapy. The CT data are used to extract the bulk tumor size and shape. The molecular imaging data they incorporated into the modeling was FLT-PET and 61Cu-ATSM-PET [84] data of head and neck cancer. This data is incorporated into an extension of the standard linear-quadratic model to describe the effects of radiation therapy. The authors’ goal was to explicitly incorporate the effects of oxygen status on the linear-quadratic model so they include an addition “oxygen-dependent modification factor,” or OMF, into the relation. Thus, the survival probability for cells being irradiated is given by
(15)P(D,pO2)=exp[-αD·OMF(pO2)·(1+D·OMF(pO2)α/β)],
where α and β are the linear and quadratic contributions, respectively, of the absorbed dose D to the induced damage. The OMF is then given by the oxygen enhancement ratio (OER) at the current (relative) partial pressure of oxygen, pO2:
(16)OER(pO2)=m·pO2+KpO2+K,
where m is the maximum OER, and K is the pO2 at the midpoint between 1 and m. The 61Cu-ATSM data is then used to determine the pO2 and the SUV of the FLT-PET data are used to initialize the tumor cell proliferation for each voxel. The main result of the paper is to employ pretreatment FLT-PET and 61Cu-ATSM data to initialize their model with the nonuniform pO2 map obtained from the 61Cu-ATSM data, a uniform pO2 of 2.5 mmHg, a uniform pO2 of 5.0 mmHg, and a uniform pO2 of 9.0 mmHg. The simulation showed that the greater the hypoxic fraction is (tissue is generally considered hypoxic if pO2 dips below 2.5 mmHg) the less efficacious the treatment is. In particular, using the inhomogeneous 61Cu-ATSM scan resulted in inhomogeneous tumor response with the least treatment benefit coming in regions of hypoxia. The authors note that this suggests that uniform dose distributions may result in an increased tumor cell survival in hypoxic regions. The authors concluded that their multiscale approach can successfully integrate molecular imaging data into a tumor growth and radiation therapy response tumor simulation.The first study to merge modeling with clinical MRI and PET data was contributed by Szeto et al. [85] who investigated the link between tumor growth kinetics and hypoxic tumor burden. The growth kinetics were examined using the approach summarized by (2), only this time the proliferation term (second term on the right hand side of (2)) was modeled by logistic growth; that is, there operational equation was
(17)∂c∂t=∇·(D(x)∇c(x,t))+ρc(x,t)(1-c(x,t)k),
where k is the carrying capacity of the section of tissue under investigation. The investigators studied the correlation between the relative hypoxic fraction, (1) cell proliferation, and (2) the modeling assessment of “biological aggressiveness.” More specifically, the relative hypoxic (RH) fraction was defined as the ratio of the hypoxic volume, as assessed by FMISO-PET, to the tumor volume defined by anatomical T2-weighted MRI. Cell proliferation is quantified by ρ, while biological aggressiveness is defined as ρ/D; both ρ and D are determined as described previously. Performing correlation analysis, they were able to show a strong (R2=0.88) and significant (P=0.000021) relationship between ρ/D and RH. Additionally, a strong (R2=0.74) and significant (P=0.00069) relationship was also found between ρ/D and the mean FMISO activity of the tumor scaled to the blood FMISO activity. The authors concluded that their data indicates that biological aggressiveness of brain tumors (as quantified by ρ/D) is linked to hypoxic burden.Building on these efforts, this group has more recently extended the model described by (17) to include the effects of hypoxia, necrosis, and angiogenesis [86]. This new model divides the population of tumor cells into three subpopulations: normoxic, hypoxic, and dead. Thus, (17) is modified into a system of PDEs with one each for normoxic, hypoxic, dead cells, tumor angiogenic factors, and vasculature. Cells are allowed to convert from one population to the other and that is determined, in part, by vascular and hypoxic status. In this preliminary effort, the authors used this approach combined with two presurgical MRI scans to predict what an FMISO image would look like just prior to surgery. The simulated image is then compared to the experimentally measured one, and it is found that several of the key features of the real image are captured by the simulation.


### 3.3. Two Points of Interest

Before leaving this section, we stress two points of interest. The first is that the overwhelming majority of the studies linking mathematical modeling and tumor growth or treatment response  have occurred in the brain. While the reasons for this are debatable (though the rigid structure of the skull and the presence of many internal fiducial markers facilitates image registration and, therefore, the application of several of the models described above), it is clearly a limitation in the current literature. Moving forward, the mathematical and biophysical modeling communities will need to address the current paucity of application of their techniques at disease sites outside of the brain. The second point worth noting is that the studies described above all have, at their core, the reaction-diffusion equation. As the excellent efforts described above convincingly indicate, this approach has much merit. But integrating the myriad of data available from quantitative imaging will most likely require new modeling approaches.


## 3.1. Modeling Studies Employing One Imaging Technique

As it is well known that glioma cells preferentially migrate along white matter fibers of the brain, it is natural to incorporate the structurala priori information on brain architecture provided by diffusion tensor imaging to guide how tumor cells will move from their source point.  This is the approach taken by Jbabdi et al. [64]. The authors began with a reaction-diffusion model of how tumor cells proliferate and migrate as follows:
(1)∂c∂t=D∇2c(x,t)+ρc(x,t),
where c(x,t) is the concentration of tumor cells at position x and time t, D is the diffusion coefficient associated with random tumor cell motion, and ρ is the proliferation rat of the tumor cells.  The first term on the right hand side of (1) is the diffusion component, while the “reaction” portion is characterized by the exponential growth term (the second term on the right hand side). As it is written, (1) does not account for the heterogeneous structure of the brain, and one way to incorporate such knowledge is to make D an explicit function of space dependent upon information provided by the DTI data. In this case, we have
(2)∂c∂t=∇·(D(x)∇c(x,t))+ρc(x,t),
where D(x) is the diffusion tensor that describes anisotropic tumor cell diffusion and therefore takes into account both location and direction as reported by the DTI data. By combining an average value for ρ (obtained from previous investigations [65, 66]) with appropriate initial (i.e., c(x,0)=co(x)) and boundary (D(x)∇c(x,t)t)·n=0 for all x on the sulcal and ventricular boundary, where n is the surface normal) conditions, the model can then be run forward to predict tumor shape and distribution which can then be compared to experimental results. In their initial effort, the authors used DTI data from a healthy volunteer, as it was not available on the patient data sets. The authors stated that their goal was to simulate tumor growth in a healthy brain to best fit the tumor size and shape that was experimentally observed in patients with brain tumors low-grade tumors. By “seeding” a tumor in the healthy brain data set and “growing” it forward in time, they were able to compare to their predicted tumors to those observed in real patients. The results indicated that when the tumors were allowed to grown anisotropically, as guided by the DTI data, the spatial extension and distribution of the tumor more closely resembled what was observed in the patient, then when the tumor cells were allowed to migrate isotropically without guidance from the DTI data. This preliminary study showed the importance of incorporating structural information available from imaging data in order to account for more realistic tumor cell growth.A similar approach to linking tumor growth and restricted movement was investigated by Clatz et al., who simulated the growth of glioblastoma multiforma (GBM) by incorporating the mechanical restrictions (based on classical continuum mechanics) presented by certain brain structures [67]. Their overall goal was to build a patient-specific simulator of GBM which includes brain deformation induced by the invading tumor. To conceptualize this, the tumor is broken down into two components: (1) the expansive, noninfiltrative component that produces the increase in tumor volume and is responsible for the mass effect and (2) the component associated with tumor diffusion that is responsible for the infiltrative component. The components are coupled in the sense that the mass effect is directly related to tumor cell density, but the tumor cell density is not influenced by the mass effect. To implement this approach, the authors first segment the images into an expansive component and a diffusion component, and the images are then registered to a brain (anatomical) atlas. Then a finite-element mesh is built based on the segmentations, and tissue properties are assigned to each mesh point including the tumor cell density based on the two segmented components. Based on this initialization, the reaction diffusion equation is run forward to generate a virtual GBM whose results are then compared to the patient data; the simulated tumor isodensity contours and brain deformation are compared to what is actually measured on MR images by registering the predicted tumor to that acquired approximately six months after the initial measurements. The algorithm was applied to two patients, and the predicted tumor sizes and shapes showed good qualitative agreement (no quantitative comparison was performed) to what was actually measured on anatomical T2-weighted MR data. These initial results were very encouraging and, importantly, showed the feasibility of using limited clinical imaging data to initialize a patient specific model of tumor growth. The authors concluded by noting that this success justified its further evaluation which then did in a follow-up paper.Bondiau et al. [68] combined the main features of the previous two paragraphs—DTI measurements of white matter tracks and mechanically restricted tumor growth—to model GBM in an excellent paper in 2008.  The DTI data is used to estimate the preferred direction of migration for the diffusion component of the reaction diffusion (2) which is then coupled to a mechanical component.  The authors again used classical linear elasticity theory to describe how the brain reacts to the tumor. In particular, they assume a linear relation between stress and strain as follows:
(3)σ=kε,
where σ is the stress tensor, k is the elasticity of the brain, and ε is the linearized Lagrange strain tensor. Since the deformation is assumed to be very slow the authors use the static equilibrium equation which states that the stress induced by tumor growth is balanced by the external force, Fext, due to tissue deformation:
(4)∇·σ+Fext=0.
The model parameters can then be tuned to match that of an MRI data obtained from a patient; that is, the different brain structures have different mechanical properties, and these are assigned in a straightforward fashion by image segmentation and literature values [69]. These assignments, together with the model can then be used to simulate tumor status at a future time point. The simulation is then compared to what is actually observed at that future time point. In this contribution, the authors were able to perform a set of quantitative comparisons including the average and odd ratio of the differences between the simulated and observed images, as well as displacement of reference points. The authors were able to show good correlation between the simulated and measured tumor masses; in particular, they stated that the inclusion of DTI data to “guide” tumor migration improves the “realism” of the diffusion component of (2). More recently, the authors have extended their approach to compare predicted growth patterns with irradiation margins in which they showed that microscopic invasion occurs past the irradiation margins and that some healthy tissue was within the radiation margins [68]. (We return to this important topic below.) Clearly, this is an exciting direction of investigation.Also using MRI data to drive the reaction-diffusion equation in GBM patients, Wang et al. employed serial pretreatment MRI data to generate patient specific rates of net glioma cell proliferation and invasion in 32 patients [70].  The results in this paper were based on previous efforts by this group in which an imaging-based method was developed to compute the parameters ρ and D from serial pre-treatment MRI data [71]. Briefly, this is done by using the serial T1-weighted (postcontrast) images in which the volumes of the tumor mass visible on these image sets were computed, and then the radius of a sphere of an equivalent volume is computed. The difference between the measured radii, along with the time interval between subsequent scans, allows for the calculation of the velocity of radial expansion, v. This velocity is then substituted into Fisher’s approximation [72]: (5)v=2Dρ.
In addition to the T1-weighted images, T2-weighted images are collected which allows for calculating D/ρ (the so-called, “index of invisibility” [71]. Given this system of two equations in two unknowns it is straight forward to compute D and ρ. Using this approach, the authors showed that their estimates of proliferation and invasion proved to be significantly related to prognosis; in particular, the ratio ρ/D was a significant predictor of survival (P=0.01). Additionally, the authors compute a new measure, termed the therapeutic response index (TRI), which is the ratio of the actual patient survival time to the predicted survival time achieved without therapeutic intervention. Thus, the TRI provides a quantitative assessment of the efficacy of a particular treatment paradigm. The authors concluded that their approach of developing patient-specific virtual controls for application in both treatment design as well as monitoring the effects of that treatment.Rockne et al. extended (2) by adding a term describing the linear-quadratic radiobiological model to describe the effects of radiation therapy [73, 74]. This model considers the surviving fraction, S, of cells as
(6)S=exp[-(αDT+βDT)],
where DT is total dose administered, and α and β are constants from experiment. In practice, DT is fractionated over n fractions, and the fractionated dose, d, is a function of both space and time, that is, d≡d(x,t). Then the biologically effective dose can be shown to be [74]
(7)BED=nd(1+dα/β).
With these assignments, equation (2) is thus updated to
(8)∂c∂t=∇·(D(x)∇c(x,t))+ρc(x,t)-R(α,d(x,t),t)c(x,t),
where
(9)R(α,d(x,t))={0,fortnotduringtherapy1-S(α,d(x,t)),fortduringtherapy}.
Using this model, the investigators could simulate how a tumor treated with radiation therapy would compare to an untreated tumor, as well as estimate a “typical response window.” The group then moved to apply this methodology to patient data. Using pretreatment MRIs to estimate patient specific ρ and D parameters as described above, Rockne et al. made predictions about what the response to radiation would be and correlated these results to clinical data [75]. In particular, they showed the proliferation rate was highly correlated to the radiation response parameter (r=0.89;P=0.0007) which allowed for the construction of a predictive relationship that was then tested in a “leave-one-out” approach. This result let the authors to conclude that their model can generate in silico tumors, with characteristics similar to the in vivo (i.e., patient) tumors and that this framework can then be used to predict individual response.Our own group has recently contributed a simple approach to integrating imaging and modeling, whereby serial DW-MRI data obtained before and early after the initiation of therapy are used to populate the logistic growth model [76–79] and predict tumor cellularity at the conclusion of therapy. The goal is to use the measurement of the ADC to first convert to estimates of cell number [16], and then use changes in the cell number to convert to a proliferation value. More specifically, we convert measurements of ADC to tumor cell number via (10) as follows:
(10)N(t)=θ(ADCw-ADC(t)ADCw-ADCmin),
where ADCw and ADCmin are the ADC of free water and the voxel with the minimum ADC, respectively; the voxel with the minimum ADC (ADCmin) contains the maximum number of cells that can fit into that voxel (i.e., the carrying capacity, θ). The voxel with the ADC of free water (ADCw) is taken to have zero tumor cells. Thus, given ADC data it is straightforward to convert to tumor cell number. If this approach is applied at multiple time points, it is possible to computer the proliferation rate, k, of the tumor on a region of interest or voxel level by fitting the N(t) data to, say, (11):
(11)N(t)=θN(t1)N(t1)+(θ-N(t1))e-kt.
(We note that the parameter k is actually a bulk (averaged) measure of both proliferation and death and can therefore take on either positive or negative values.) Our first effort [77] using this simple approach was in a small animal model of cancer in which 13 rats were given 9L gliosarcoma brain tumors; eight of the rats were treated with the chemotherapeutic drug 1,3-bis(2-chloroethyl)-1-nitrosourea, and five rats were untreated controls. All animals underwent DW-MRI immediately before, one day and three days after treatment. The ADC data from the first imaging time points were used to estimate k as just described which was then used to predict the number of tumor cells at data three, and this value was correlated with the day five experimental data.  On a region of interest basis, the Pearson’s correlation coefficients between the predicted and measured tumor cellularity was 0.88. Thus, this simple approach has some ability to capture the gross (i.e., at the ROI level) cellularity changes.We also applied this approach to breast cancer data acquired during neoadjuvant chemotherapy. Six patients received DW-MRI before, after one cycle, and after all cycles of neoadjuvant chemotherapy. Again, the proliferation rates were estimated using the ADC data from the first two time points and then used with (10) and (11) to predict cellularity after therapy. When comparing the predicted to the experimentally observed data, the Pearson’s correlation coefficient for the ROI analysis yielded 0.95 (P=0.004), and, encouragingly, after applying a conservative 3 × 3 mean filter to the ADC data, the voxel-by-voxel analysis yielded a Pearsons correlation coefficient of 0.70±0.10 (P<0.05). More recently, we have had some preliminary success in using this approach to separate responders from nonresponders [80].Another effort that has incorporated serially acquired DW-MRI data into a mathematical model was contributed by Ellingson et al. [81]. As there is an established negative correlation between ADC and cellularity [16], the ADC can be directly substituted into (2) to capture
(12)ddtADC=∇·(D∇ADC)+ρADC,
and, subsequently,
(13)ddtADC=D∇2ADC+∇D·∇ADC+ρADC.
Algebraically solving (13) for ρ at time point n-1 allows for an estimate of ρ; that is,
(14)ρn-1=(1ADCn-1)×(ddtADCn-1-D∇2ADCn-1-∇D·∇ADCn-1),
where d/dt of ADCn-1 defines the rate of change of ADC from the previous time point n-1 to the first time point n-2. After ρ is fixed, it is substituted back into (12) to yield a PDE that can be solved to estimate D. The authors applied this approach to obtain estimates of ρ, D, and v in 53 patients. In nine of the patients, magnetic resonance spectroscopy was performed to measured Cho/NAA (i.e., choline to N-acetylaspartate ratio); the Pearson’s correlation coefficient between ρ and Cho/NAA was 0.97 (P<0.0001). Furthermore, using a “hot-spot analysis” (based on abnormal anatomical MRI regions), the authors were able to use ρ and D to statistically separate (World Health Organization, WHO) grade II from grades III and IV. The authors then extended this effort to predict progression free and overall survival (OS) in 26 recurrent glioblastoma patients treated with the antiangiogenic drug bevacizumab [82]. Again, the authors were able to find strong correlations between ρ and progression-free survival (P<0.0001) and overall survival (P=0.0093).


## 3.2.  Modeling Studies Employing Multiple Imaging Techniques

One of the first studies to include multimodality imaging in a mathematical model of tumor growth was put forth by Titz and Jeraj in 2008 [83]. In this study, the authors employed FDG-PET/CT data to investigate the effects of oxygenate status on the response to radiation therapy. The CT data are used to extract the bulk tumor size and shape. The molecular imaging data they incorporated into the modeling was FLT-PET and 61Cu-ATSM-PET [84] data of head and neck cancer. This data is incorporated into an extension of the standard linear-quadratic model to describe the effects of radiation therapy. The authors’ goal was to explicitly incorporate the effects of oxygen status on the linear-quadratic model so they include an addition “oxygen-dependent modification factor,” or OMF, into the relation. Thus, the survival probability for cells being irradiated is given by
(15)P(D,pO2)=exp[-αD·OMF(pO2)·(1+D·OMF(pO2)α/β)],
where α and β are the linear and quadratic contributions, respectively, of the absorbed dose D to the induced damage. The OMF is then given by the oxygen enhancement ratio (OER) at the current (relative) partial pressure of oxygen, pO2:
(16)OER(pO2)=m·pO2+KpO2+K,
where m is the maximum OER, and K is the pO2 at the midpoint between 1 and m. The 61Cu-ATSM data is then used to determine the pO2 and the SUV of the FLT-PET data are used to initialize the tumor cell proliferation for each voxel. The main result of the paper is to employ pretreatment FLT-PET and 61Cu-ATSM data to initialize their model with the nonuniform pO2 map obtained from the 61Cu-ATSM data, a uniform pO2 of 2.5 mmHg, a uniform pO2 of 5.0 mmHg, and a uniform pO2 of 9.0 mmHg. The simulation showed that the greater the hypoxic fraction is (tissue is generally considered hypoxic if pO2 dips below 2.5 mmHg) the less efficacious the treatment is. In particular, using the inhomogeneous 61Cu-ATSM scan resulted in inhomogeneous tumor response with the least treatment benefit coming in regions of hypoxia. The authors note that this suggests that uniform dose distributions may result in an increased tumor cell survival in hypoxic regions. The authors concluded that their multiscale approach can successfully integrate molecular imaging data into a tumor growth and radiation therapy response tumor simulation.The first study to merge modeling with clinical MRI and PET data was contributed by Szeto et al. [85] who investigated the link between tumor growth kinetics and hypoxic tumor burden. The growth kinetics were examined using the approach summarized by (2), only this time the proliferation term (second term on the right hand side of (2)) was modeled by logistic growth; that is, there operational equation was
(17)∂c∂t=∇·(D(x)∇c(x,t))+ρc(x,t)(1-c(x,t)k),
where k is the carrying capacity of the section of tissue under investigation. The investigators studied the correlation between the relative hypoxic fraction, (1) cell proliferation, and (2) the modeling assessment of “biological aggressiveness.” More specifically, the relative hypoxic (RH) fraction was defined as the ratio of the hypoxic volume, as assessed by FMISO-PET, to the tumor volume defined by anatomical T2-weighted MRI. Cell proliferation is quantified by ρ, while biological aggressiveness is defined as ρ/D; both ρ and D are determined as described previously. Performing correlation analysis, they were able to show a strong (R2=0.88) and significant (P=0.000021) relationship between ρ/D and RH. Additionally, a strong (R2=0.74) and significant (P=0.00069) relationship was also found between ρ/D and the mean FMISO activity of the tumor scaled to the blood FMISO activity. The authors concluded that their data indicates that biological aggressiveness of brain tumors (as quantified by ρ/D) is linked to hypoxic burden.Building on these efforts, this group has more recently extended the model described by (17) to include the effects of hypoxia, necrosis, and angiogenesis [86]. This new model divides the population of tumor cells into three subpopulations: normoxic, hypoxic, and dead. Thus, (17) is modified into a system of PDEs with one each for normoxic, hypoxic, dead cells, tumor angiogenic factors, and vasculature. Cells are allowed to convert from one population to the other and that is determined, in part, by vascular and hypoxic status. In this preliminary effort, the authors used this approach combined with two presurgical MRI scans to predict what an FMISO image would look like just prior to surgery. The simulated image is then compared to the experimentally measured one, and it is found that several of the key features of the real image are captured by the simulation.


## 3.3. Two Points of Interest

Before leaving this section, we stress two points of interest. The first is that the overwhelming majority of the studies linking mathematical modeling and tumor growth or treatment response  have occurred in the brain. While the reasons for this are debatable (though the rigid structure of the skull and the presence of many internal fiducial markers facilitates image registration and, therefore, the application of several of the models described above), it is clearly a limitation in the current literature. Moving forward, the mathematical and biophysical modeling communities will need to address the current paucity of application of their techniques at disease sites outside of the brain. The second point worth noting is that the studies described above all have, at their core, the reaction-diffusion equation. As the excellent efforts described above convincingly indicate, this approach has much merit. But integrating the myriad of data available from quantitative imaging will most likely require new modeling approaches.


## 4. Future Directions

There are three key areas (two theoretical and one experimental) that are primary goals for investigators working at the interface of biomedical imaging and tumor modeling: (1) the need to construct models amenable to the incorporation of imaging data, (2) the construction of a framework, informed by clinical and cancer biology, that allows forin silico, individualized clinical, trials, and (3) the design and execution of appropriately designed experiments to test the theoretical developments of areas (1) and (2).The majority of existing biomathematical and biophysical models of tumor growth and treatment response are not of the form that are readily amenable to incorporating imaging data; more specifically, the models are not of the kind that can readily be populated by data that can be measured in an intact system at multiple time points with reasonable spatial resolution in 3D. Frequently, the models have a wide variety of parameters that must be estimated or assumed based on measurements from disparate sources. Conversely, the utility of imaging data is that it can be acquired on a patient specific basis so that models can be initialized and constrained by an individual patient’s characteristics. As the majority of current models are not constructed with this in mind, this is a fundamental—though artificial—barrier to progress on integrating imaging data in tumor growth models. Thus, going forward, models need to be constructed with the types of measurements that are available from imaging data in mind. In the author’s view, this should be one of the primary goals for investigators working in this field.Once tumor growth models are built that are readily populated by the imaging data described above in Section3, it becomes immediately possible to perform in silico experiments to predict, given a set of possible treatment regiments, an optimal treatment strategy. For example, in a well-vascularized tumor it would be important to determine the order and timing of (say) an antiangiogenic drug and a cytoxic drug (for an excellent study on this point, the interested reader is referred to [87]). With an appropriate model in hand, simulations could be run in which the order, dosing, and timing of the antiangiogenic and cytoxic drugs were varied, and this would, presumably, lead to different outcomes (i.e., predictions). Thus, the order, dosing, and timing of the drugs could be varied and explored as parameters to be optimized to yield a prediction of what treatment protocol would be most effective. Of course, such studies would require intimate involvement from experts in cancer biology as well as clinical oncologists with expertise in more standard clinical trials. But if such studies were carried out, they would make concrete predictions that are eminently testable in both the preclinical and clinical settings. In this way, the synthesis of imaging and modeling would truly enable the development of a field of “theoretical oncology.”As imaging data is acquired with minimal interaction with the system under investigation, measurements made early in an experiment can lead to predictions that can then be tested later in the experiment. More specifically, in the context of pre-clinical rodent tumor studies, imaging measurements could be made at baseline (i.e., after an implanted tumor reaches a particular size but before therapy commences), then at a second time point early in the course of treatment. These two data sets could then be used to initialize the model(s) to determine the optimal treatment approach as outlined in the previous two paragraphs. The optimized therapy could then be given to the (say) mouse and the disease progression/regression tracked. In this way the predictions of the imaging based model can be directly compared to experiment to determine theirin vivo utility and validity. There is no other currently available approach that enables this. Indeed, it is hard to imagine a methodology such as this that does not involve imaging since it is all based on making spatially resolved, noninvasive measurements; any method that can achieve this is some brand of imaging, almost by definition. Furthermore, there is, of course, no reason the same paradigm cannot be applied to human clinical trials. Thus, an appropriate model, initialized by patient specific imaging data, could be used to perform a wide range of in silico clinical trials to determine an optimized, personalized treatment approach for patients. Such a development would hasten the arrival of personalized medicine for oncology patients.


## 5. Summary and Conclusion

After reviewing a representative sampling of the emerging, quantitativein vivo imaging data that is readily available, we proceeded to show such data can potentially be used to initialize and constrain biomathematical and biophysical models of tumor growth and treatment response. As quantitative, imaging data can be acquired noninvasively, in 3D, and at multiple time points, it is a natural and powerful method for acquiring data that can directly populate appropriate models. Furthermore, precisely since the data are acquired noninvasively, the model predictions can then be directly compared to what actually happens in the subject from which the data were taken. Thus, this approach lends itself to make patient specific predictions that can be experimentally tested. It is the author’s hope that this paper will motivate others to work at the exciting and clinically relevant, interface of quantitative cancer imaging and modeling of tumor growth and treatment response. If such a line of investigation proves successful, then it has the ability to fundamentally shift existing paradigms of therapy selection and monitoring in cancer.



---

*Source: 287394-2012-12-26.xml*