|Home | About | Journals | Submit | Contact Us | Français|
Quantitative imaging data obtained from multiple modalities may be integrated into mathematical models of tumor growth and treatment response to achieve additional insights of practical predictive value. We show how this approach can describe the development of tumors that appear realistic in terms of producing proliferating tumor rims and necrotic cores. Two established models (the logistic model with and without the effects of treatment) and one novel model built a priori from available imaging data have been studied. We modify the logistic model to predict the spatial expansion of a tumor driven by tumor cell migration after a voxel’s carrying capacity has been reached. Depending on the efficacy of a simulated cytoxic treatment, we show that the tumor may either continue to expand, or contract. The novel model includes hypoxia as a driver of tumor cell movement. The starting conditions for these models are based on imaging data related to the tumor cell number (as estimated from diffusion-weighted MRI), apoptosis (from 99mTc-Annexin-V SPECT), cell proliferation and hypoxia (from PET). We conclude that integrating multi-modality imaging data into mathematical models of tumor growth is a promising combination that can capture the salient features of tumor growth and treatment response and this indicates the direction for additional research.
In recent years there have been dramatic increases in the range and quality of information available from non-invasive imaging methods, so that several potentially valuable imaging measurements are now available to quantitatively measure tumor growth, assess tumor status and predict treatment response. Magnetic resonance imaging (MRI), single photon emission computed tomography (SPECT) and positron emission tomography (PET) have matured to the point where they offer unique measures of tumor status at the physiological, cellular and molecular levels. Imaging techniques can provide non-invasive, quantitative and 3D characterizations of, for example, blood flow, vessel permeability, cellularity, cell membrane turnover, pH and pO2 (MRI); cell surface receptor expression, apoptosis (SPECT); and metabolism, cell proliferation and hypoxia (PET). These methods are commonly used in clinical trials and are slowly being incorporated into clinical practice; for example, fluordeoxyglucose PET (FDG-PET) has now been included in the most recent version of the response evaluation criteria in solid tumors (RECIST 1.1 (Eisenhauer et al 2009)) for assessing treatments.
Independent of these developments have been the major advances in the mathematical modeling of tumor growth. Such models have a wide range of complexity ranging from exponential growth to mechanical models of tumor growth, to systems of partial differential equations describing invasion and angiogenesis (Mantzaris et al 2004). However, current models are frequently limited in their practical applicability as they require input data measuring, for example, chemotaxis, haptotaxis or growth factor gradients—data that are very difficult to obtain in an intact organism with any reasonable spatial resolution. Consequently, there has been very little application of such models to clinical data sets and virtually no incorporation in clinical trials. In light of this realization a few groups are exploring the concept of using imaging data to drive mathematical models. The motivation for integrating imaging data into mathematical models of tumor growth is that imaging can provide information non-invasively, in 3D and at multiple time points. Thus, it may be possible to model the entire system without disturbance and to check the predictions of the model in individual animals or humans. Measurements can be made at, for example, the time of diagnosis and then early in the course of treatment and these changes could be modeled to predict response at the end of therapy. In this way, imaging allows models to be initialized with patient- or animal-specific data, as opposed to models that are driven from data obtained from assays or groups of animals. Furthermore, since the systems are left intact, the models can be updated with more measurements taken at later time points.
A major difficulty is being unable to parameterize a reasonable, realistic model of tumor growth by data that can be obtained from non-invasive imaging data. Typical models of tumor growth do not include parameters that are (currently) obtained by practical imaging, though there have been several recent efforts to accomplish this. For example, Swanson et al have extended the reaction-diffusion model of cell movement to model brain tumor growth as measured by MRI and have used this approach to compute a net invasion rate (Swanson et al 2000). In a more recent effort, they correlated such data to hypoxia measurements obtained from PET (Szeto et al 2009). Another example was offered by Hogea et al who presented a framework for modeling glioma growth and their mechanical impact on surrounding tissue (Hogea et al 2007). Their efforts are based on a reaction diffusion model for tumor growth and an elastic material for the background. They showed the development of a model with predictive capabilities based on model parameters estimated for an individual patient. While other efforts have been reported (Bondiau et al 2008, Zhang et al 2007), this is clearly a young and exciting field that shows much promise.
Here we explore in simulations, how quantitative imaging data obtained from multiple modalities may serve as initial conditions for simple mathematical models of tumor growth and treatment response. We incorporate measurements of cellularity via apparent diffusion coefficient (ADC) mapping by diffusion-weighted MRI (DW-MRI), cell proliferation via fluorodeoxythymidine PET (FLT-PET), apoptosis via 99mTc-Annexin-V SPECT and hypoxia imaging via FMISO-PET into three models of tumor growth. Detailed descriptions of these techniques can be found elsewhere (Bihan et al 2001, Charles-Edwards and deSouza 2006, Shields et al 1998, Muzi et al 2005, Blankenberg et al 1999, Takei et al 2004). We show how the readily available in vivo imaging data can be combined with simple mathematical models to generate tumors that display realistic growth patterns that capture several salient features of tumor growth and treatment response. We do this by studying modified versions of the logistic model (with and without the effects of treatment) and a novel model that has been built a priori from imaging data.
The logistic growth model incorporates exponential growth of the number of tumor cells early in time, which then asymptotically approaches the limiting cellular carrying capacity. This model is explained in detail in Byrne (2003) and the relevant differential equation is
where N(r, t) is the number of cells per voxel at position r and time t, k(r) is the cell proliferative rate (in units of 1/day) at position r and θ is the carrying capacity of the population. The solution is
where N(r, 0) is the number of cells present at t = 0. It is straightforward to include treatment within the logistic growth law (Byrne 2003). We merely amend equation (1) by adding a term that characterizes the concentration of a chemotherapeutic drug present at time t and position r, A(r, t):
where μ is the rate at which the drug kills tumor cells. A(r, t) is assumed to be a periodic function that is constant in the time interval when the drug is present in the tumor.
If we make the following assignments:
it is clear that equation (4) is an ordinary differential equation of the Bernoulli type which can be linearized by setting
Solving equation (8) yields
and replacing z with N returns
Thus, the solution to equation (3) is therefore given by
We note that equation (13) is valid only when A(r, t) is piecewise constant, periodic function in the time intervals in which the drug is present at the tumor site. Below, we show how equations (2) and (13), with some modification, can be used to model tumor growth and treatment response.
We consider the cells within a voxel as proliferating, undergoing apoptosis or migrating into or out of the voxel. Mathematically, we write
where Ni(t) is the number of tumor cells in the ith voxel at time t, μp,i and μd,i are the proliferation and death rates of the tumor cells in the ith voxel, respectively, kim represents the net ‘transfer constant’ of cells from the ith voxel to the mth voxel. This is the transfer between the ith voxel and the voxel that produces the maximum absolute hypoxia difference between itself and the ith voxel in a two-dimensional space. Here we assume that kim is proportional to the difference in pO2 as estimated by FMISO-PET, though it could be easily amended to indicate other factors that determine cell movement. In terms of the accumulation of a PET tracer as estimated by the standard uptake value (SUV (Eschmann et al 2007)), we have
where the terms on the right-hand side of equation (15) represent the SUV of FMISO in the ith and mth voxels as given by
and we have assumed that cells will tend to move in the direction of increased pO2. The kim term determines the values of N(t): if kim < 0, then N(t) = Ni(t) and if kim > 0, then N(t) = Nm(t). In this model we can use imaging to provide estimates of cell proliferation (μp,i via FLT-PET), apoptosis (μd,i via 99mTc-Annexin-V) and hypoxia (kim via FMISO-PET). If this is combined with an initial estimate of the cell number (Ni via DW-MRI), then the model can make predictions for both the number and spatial distribution of tumor cells as a function of time.
We aim to use imaging data acquired early in a tumor’s life cycle to predict the shape and distribution of tumor cells at a later time point. In this way the imaging data serve as initial conditions in which to initialize the models described in section 2.1. We ‘grow’ in silico tumors in a candidate rat brain represented by a 3D MRI data set that is naturally gridded by voxels. We use those discrete elements to grow a tumor that begins as a small number of voxels and grows outward. We take our voxels to be 100 × 100 × 100 μm3. The first step in this process is to make a reasonable correspondence between the terms in equations (2), (13) and (14), and the data obtained from imaging experiments. For example, in equation (2), we need to make assignments for N(r, 0), k(r) and θ. Similarly, in equation (13), we need to make assignments for N(r, 0), k(r), θ, μ(r) and A(r, t). In equation (14), we must make assignments for Ni, Nm, μp,i, μd,i and kim.
The carrying capacity, θ, can be interpreted as the maximum number of cells that can be contained within a given volume. For imaging, the relevant volume scale is a voxel and therefore θ = (voxel volume)/(cell volume). The voxel volume and the individual cell volume (for a given cell species) are relatively well-defined quantities, so θ can be calculated. In the simulations below, we take the cell volume to be 1150 μm3 (Das et al 2005), the voxel size to be 1 × 10−6 μm3; thus, the carrying capacity is 869 cells per voxel.
In controlled situations, it has been shown that ADC values vary inversely with the cell volume fraction. Here we use the mean tumor ADC values reported by Eis et al (1995) to describe the initial distribution of cell density within the tumor. We require one ADC value for each imaging voxel (i.e. an ADC map) so we construct a Gaussian noise vector of length equal to the number of voxels present in the tumor at time zero and add the mean ADC value to each entry in the vector. The resulting values are then distributed throughout the tumor to establish the initial ADC map. The standard deviation of the Gaussian noise was assigned as the standard deviation reported by Eis et al (1995). From a variety of studies in cells and tissues, we can also write
where ADCw is the ADC of free water and λ is a proportionality constant that varies with the cell type. This relationship has been shown to hold in well-controlled settings (Anderson et al 2000) and can be re-arranged to give N(r, t) in terms of ADC(r, t), ADCw, which can be directly measured from DW-MRI experiments. In order to convert from ADC to cell number, we need to determine the calibration constant, λ. This is done by taking the ADC of free water to be 3.0 × 10−3 mm2 s−1 (Hagmann et al 2005) and relating this to a cellular volume fraction of zero. Then the ADC value at a cellular volume fraction of 1 is taken to correspond to the carrying capacity θ which is the maximum number of cells that can fit into a voxel; this ADC was assigned to be 1.1 × 10−3 mm2 s−1 which was the maximum value reported by Anderson et al (2000). The ADC map is then converted to a cell number by assuming equation (17) with λ = 2.19 × 10−6 (mm2/s*cells). This calibration is used to assign N (r, 0), Ni and Nm in equations (2), (13) and (14).
Since FLT-PET reports on cell proliferation (in specific circumstances; see, e.g. Salskov et al (2007)), we can obtain estimates of k(r) and μp,i from FLT-PET images. By performing compartmental modeling on FLT kinetics, it is also possible to extract various parameters that relate to FLT delivery and retention. A parameter of central interest in these analyses is the metabolic flux of fluorothymidine, denoted by the parameter KFLT in units of mL min−1 g−1. In the simulations below, the proliferation rate k(r) was estimated from the mean KFLT values reported by Bradbury et al (2008) in their study of the DF-1 cells in the TV-a/Ink4 a−/−Art−/− mouse which is a glioma tumor model. Since we require one value of KFLT for each imaging voxel (i.e. a KFLT map), we construct a Gaussian noise vector of length equal to the number of voxels present in the tumor at time zero and added the mean KFLT value to each entry in the vector. The resulting values are then distributed throughout the tumor to establish the initial KFLT map. The standard deviation of the Gaussian noise was assigned as the standard deviation reported by Bradbury et al (2008). We take the KFLT value to be proportional to the proliferation rate of the tumor and this is used to generate the k(r) map. The proliferation rate values, k(r), were taken from the range reported by Kiran et al (2009). The proliferation rate that gave reasonable results from the range of proliferation rates was chosen. The k(r) chosen was 1/day with a standard deviation of 10% for the logistic model, 0.3/day for the logistic model with drug treatment for the responding tumor treatment, 0.6/day for the logistic with drug treatment for the non-responding tumor and 0.5/day for the proliferation-apoptosis-migration model. The standard deviation of 10% was chosen to match the percentage standard deviation of the KFLT values. Specifically, the individual k(r) values for each voxel are determined by linearly mapping the maximum KFLT value to a maximum k(r) value and the minimum KFLT value to a minimum k(r) value. The maximum and minimum k(r) values are defined as the mean k(r) value plus or minus three times the standard deviation of the k(r) values, respectively. These two points are then used to determine a line through which all other KFLT values are mapped to k(r) values. This is summarized by
The same approach is employed to assign μp,i in equation (14).
We simulate a case of tumor shrinkage by delivering a drug to our simulated tumor. We assume the tumor to be most well vascularized in the periphery and less so toward the core so that the drug concentration within the tumor has approximate spherical symmetry with a maximum concentration on the periphery. The total number of tumor cells killed by the drug depends on a combination of its kill rate and total concentration so cells on the periphery of the tumor will suffer the greatest destruction. A linear relationship between the drug dosage and the cell death rate is assumed whereby an increase in the dosage of the drug leads to an increase in the death rate (Powathil et al 2007). Therefore, the drug delivered A(r, t) in equation (13) is dimensionless and equal to 1 and the death rate for that dosage is μ(r) in equation (13).
The death rates μ(r) in equation (13) and μd,i in equation (14) were estimated from the SUVAnV of 99mTc-Annexin V in a rat brain tumor treated after a single dose of single 150 mg kg−1 dose of cyclophosphamide (Takei et al 2004). We can assume that μ(r) is proportional to the maximum accumulation of the 150 mg kg−1 of cyclophosphamide dose and concentrations lead to a fraction of μ. In generating the SUVAnV map, we proceed as above for the FLT and ADC maps; namely, the mean reported SUVAnV value for the tumor is summed with a Gaussian vector of length equal to the tumor size and with a standard deviation equal to that reported by Takei et al (2004). These values are then distributed within the tumor voxels. The mean death rate chosen was 2.1/day (Sinek et al 2004) with a standard deviation of 10%. This was chosen to match the percentage standard deviation of the SUVAnV. The SUVAnV is assumed to be proportional to the tumor cell death rate and this is used to generate the μ(r) and μd,i maps. The SUVAnV map was then converted to a death rate map by making the maximum death rate (mean death rate + 3 × standard deviation of the noise) correspond to maximum SUVAnV and minimum death rate (mean death rate − 3 × standard deviation of the noise) correspond to minimum SUVAnV. These two points are then used to determine a line through which all other SUVAnV values are mapped to μ(r) values. This is summarized by
The SUV of 18F-fluoromisonidazole (FMISO), SUVFMISO, was used to generate a hypoxia map. Similar to above approaches, the mean and standard deviation of the SUVFMISO values reported in the literature (Eschmann et al 2007) were used to generate a map of the hypoxia in different voxels. The net transfer of cells kim in equation (15) was calculated from the SUVFMISO of the voxel of interest and the adjacent voxels. The voxel that produces the highest absolute SUVFMISO is the mth voxel in equation (16), and kim is then calculated using this voxel as shown in equation (15). This supposes that the cells will migrate to the voxel with the lowest pO2 value. Again, if hypoxia does not significantly determine tumor cell migration (this is a topic of some debate), it can be changed to model other phenomena in a straightforward manner. The main point is that hypoxia (or its suitable replacement) can be imaged and therefore provide initial conditions into the tumor growth model.
We perform simulations on three models: the modified logistic model, the modified logistic model with treatment and the novel proliferation-apoptosis-migration model. Table 1 summarizes the parameters used for the models.
A voxel size of 100 μm was chosen as the spatial resolution obtainable by standard small animal MRI methods in a reasonable time period (i.e. approximately 1 h). The tumor cell volume was set at 1150 μm3 (Das et al 2005). Given the voxel size and the tumor cell volume, the carrying capacity, θ, is calculated. We initialize the 2D simulations by seeding a tumor 5 × 5 voxels in size, and initial conditions (i.e. the cell density and proliferation maps obtain from ADC and KFLT data) are assigned to the tumor as described in the previous sections. The simulation is then run for 80 days with a time step of 1 day. At each step, if the number of cells within a given voxel is greater than 0.90 θ, the cells are moved randomly to the adjacent neighbors. This is done by calculating the number of cells above 0.90 θ for each voxel and for each cell randomly selecting any of the three to eight neighboring voxels (depending on the position of the voxels in relation to the boundary) and moving the cell into that voxel and adjusting the number of cells in that voxel. This is done until no voxel has more than 0.90 θ cells. As the tumor spreads into a new voxel, a proliferation value is calculated for the new voxel with tumor cells by adding Gaussian noise to the mean proliferation value as described in section 2.2.3.
A flowchart outlining the entire simulation process is shown in figure 1. As the tumor cells proliferate and spread to nearby voxels, cells in a voxel become quiescent when they are more than 100 μm from the tumor periphery and are also surrounded by adjacent voxels with cell numbers greater than 0.90 θ. This was chosen because it has been shown that in four malignant glioma cell lines, Ki67 staining (which is an immunohistochemical measure of cell proliferation) revealed proliferation was restricted to an outer zone or periphery of 100 μm in width (Bell et al 2001). Tumor cells that are between 100 μm and 300 μm from the periphery are assumed to be quiescent and tumor cells become necrotic and exponentially die when they are at least 300 μm deep and surrounded by voxels with cell numbers greater than 0.90 θ. The necrotic cells die off with an exponential decay according to this rule:
where Nnecrotic(t) is the number of necrotic cells at time t, Nnecroticinitial refers to the number of cells at the time the voxel becomes necrotic and μnecrosis is the rate at which necrotic cells die. The death rate due to necrosis was set at 0.17/day, which is an average of four different tumor cell lines (Busini et al 2007). The dead cells are assumed to be removed by natural processes so that they no longer contribute to the total tumor cell number in each voxel. The quiescent and necrotic cells do not proliferate and therefore do not increase the total number of cells at each iteration.
For this model, we use the same method as above until the tumor has grown for approximately 80 days. The cell densities and proliferation rates at that time point then furnish the initial conditions for the drug treatment model. The concentration of the drug over time is assumed to be pulsed with the drug concentration equal to 1 (A is dimensionless) on the days the drug is delivered and zero on the days the drug is not delivered to the tumor cells. Also, since the time step of the simulation is in days and the clearance of the drug is in hours, the drug concentration can be assumed to be instantaneous so that the concentration of the drug is 1 on the day the drug is applied to the tumor. The concentration of the drug varies with radial symmetry, with the maximum and minimum concentrations occurring in the periphery and center of the tumor, respectively. The equation describing the concentration of the drug as a function of distance is as follows:
where dist(x, y) is the Euclidean distance between the tumor center and the point (x, y) and A(r, t) is the concentration of the drug at position r and time t as in equation (13). The tumor center is calculated by taking the maximum tumor length in the x and y directions and using these values to determine the tumor center; i.e., we simply divide the x and y lengths by 2.
We evaluated cyclophosphamide as the cytoxic drug and assume that it acts only on the proliferating cells. As the proliferating cells are killed according to equation (13) and are removed, the quiescent cells become less than 100 μm from the periphery and thus now have access to oxygen and nutrients to begin proliferating. The proliferating cells are then acted upon by the chemotherapy drugs. Also, note that the distribution of the drug changes with the shrinking tumor according to equation (21), though the radial symmetry is maintained. Figure 2 details the steps for the modified logistic model with drug treatment. By changing the timing of the application of the drug and the proliferation rate of the tumor cells, we can simulate effective and ineffective treatments.
The same voxel size, carrying capacity, time step and initial condition methods that are employed for the modified logistic models are used for the PAM model. Cells also become necrotic and quiescent depending on their position and the number of cells that surround them as described in section 2.4.1. Again, a tumor measuring 5 × 5 voxels initializes the simulation with initial condition maps generated as described above. For each cell in each voxel, we define and calculate a ‘hypoxic load’, H, which is a summary metric of oxygen present for each voxel. H is calculated by dividing the SUVFMISO in a voxel by the number of cells in that voxel. At each iteration, the cells grow, die, or are dispersed to other voxels as described in section 2.3.1 As cells move, their H value travels with them to the next voxel which is then used to update the total oxygenation status of the voxel. After the movement of the cells, the voxels that have a cell number greater than 0.90 θ are identified and the cells are then moved to the nearest adjacent cells based on the hypoxia of the adjacent cells and a probability; the cells have an 80% chance of migrating to the voxel with the lowest hypoxia, a 10% probability of going to the voxel with the next lowest hypoxia and a 10% probability of being dispersed to the rest of the adjacent cells. (These probability values, though reasonable, are arbitrary.) If hr (i, t) is the hypoxic load of the rth cell in voxel i at time t, Ni(t) is the total number of cells in voxel i at time t, and Hi(t + 1) is the hypoxic load of the cells in voxel i at the next iteration, then
where n is the total number of cells in the voxel. Figure 3 depicts this simulation process.
Figure 4 depicts a typical set of ‘imaging initial conditions’ which begin the simulations; again, by initial conditions we mean the set of imaging data that are available when the tumor is first realized. Panel A is a simulated map of the PET parameter KFLT which is a measure of the overall activity of TK1 and therefore represents an estimate of tumor cell proliferation within the tumor. TK1 is the abbreviation for ‘thymidine kinase 1’ which is an enzyme that phosphorylates FLT and thereby traps it within the intracellular space. We map this parameter into a ‘proliferation index’ as discussed above and the result is displayed in panel B. Similarly, C is a simulated ADC map which is mapped to a ‘cell number’ as discussed above and the result is shown in panel D. These maps are then input into equation (2) and allowed to run forward in time.
The bottom row of figure 5 shows the proliferation rate of the cells at each of the five time points (days 0, 20, 40, 60 and 80, respectively) studied in the simulation. Each panel represents the central ‘slice’ through the three-dimensional tumor. The top row presents the results of implementing the proliferation initial conditions to drive tumor growth and the color bar indicates the number of cells within each voxel at each time point. Note that as the tumor enlarges, some cells become quiescent and some voxels become necrotic yielding different tumor cell ‘layers’. The quiescent cells do not proliferate and the necrotic cells die off according to equation (20) and are removed from the tumor mass thereby resulting in a decreasing number of cells and a proliferation rate of zero in the core of the tumor. The overall effect of the model is to produce a growing tumor with an expanding, proliferating rim and a steadily necrosing central core—not unlike what is seen in many actual brain tumor models.
Next we simulate the effects of treatment on a tumor mass developed from the modified logistic model. Figure 6 presents the initial conditions used to initialize the modified logistic model with the incorporation of the effects of treatment. Panels A and B are the initial maps of KFLT and the corresponding proliferation rate, respectively. Panel C is the SUVAnV which can be interpreted as an overall measure of cell death due to apoptosis. This is transformed to the death rate of the tumor shown in panel D as described by equation (19). Panel E displays the distribution and concentration of normalized drug within the tumor. We assume that the highest concentration of the drug is achieved on the periphery with a decreasing amount toward the core of the tumor as described by equation (21). This distribution is then related to the μ · A(r, t) as described in section 2.2.4. We present two illustrative cases: (1) the treatment is effective and the tumor shrinks (figure 7) and (2) the treatment is ineffective and the tumor continues to grow (figure 8).
Figure 7 depicts the central ‘slice’ of the tumor as it responds to effective therapy (which is assumed to act on only actively proliferating cells) over a 40 day time period. The top row presents the normalized distribution of the cytoxic drug, A(r, t), that is administered daily, while the second and third rows present the proliferation and the death rate maps, respectively. As the proliferating cells die, they free up resources (oxygen/nutrients) and the quiescent cells can transform into proliferating cells that can then be acted upon by the drug. The bottom row depicts the corresponding number of tumor cells during the treatment window. It is clear that this combination of initial conditions and treatment leads to an effective reduction in tumor size.
Figure 8 again shows the central slice of a tumor grown via equation (2), but this time the treatment regimen is not effective and the tumor does not respond over a 40 day period. The top row presents the normalized distribution of the cytoxic drug, A(r, t), that is administered every 7 days, while the second and third rows present the proliferation and the death rate maps, respectively. The bottom row depicts the corresponding number of tumor cells during the treatment window. The cell death observed is mainly determined by necrosis and therefore there are fewer cells in the core than the periphery. The tumor does not shrink but actually grows larger as seen in columns B–E. Ineffective tumor treatment can be simulated by increasing the proliferation rate of the tumor and also by increasing the dosing interval. In the case below the proliferation rate was doubled from 0.3/day to 0.6/day and the dosing interval was lengthened to once every 7 days.
Proceeding as with the previous two models, figure 9 depicts a typical set of ‘initial conditions’ which initialize the PAM simulations. Similar to figures 4 and and6,6, panels A–G correspond to the initial condition data provided by a multi-modality imaging experiment. Panels A and B are the parametric maps of KFLT and the corresponding proliferation index, respectively. Panels C and D are the SUVAnV, and the corresponding death rate, respectively. Panels E and F are the ADC map, and the corresponding cell number map, respectively. Finally, panel G is a map of the SUVFMISO. These maps are then used to initialize equation (14) which is then allowed to run forward in time over a period of 16 days.
Figure 10 presents the time course of each parameter employed in the PAM model as a function of time. The first row displays the spatial distribution of the cell proliferation rate, while rows two and three present the death rate and hypoxia maps for the tumor, respectively, as a function of time. The last row in the figure is the resulting distribution of the cell number as a function of time given the conditions displayed in the other three rows. As the tumor enlarges, some cells become quiescent and some cells become necrotic as seen (first) in panel B. The quiescent cells do not proliferate but can be moved to other voxels based on the hypoxia difference. If conditions become favorable, the quiescent cells can then become proliferating cells. Due to the migration of the tumor cells, the PAM model produces a more aggressive tumor than the modified logistic model as manifested by the time scale. The proliferation rate of the PAM model is half the proliferation of the logistic model but after 14 days, the PAM model had a tumor volume that was approximately two times larger than the tumor volume of the logistic model.
We have introduced, and explored in simulations, the concept of using quantitative imaging data obtained from multiple modalities as initial conditions to drive simple models of tumor growth and treatment response. The inputs to the simulation include cellularity obtained via ADC mapping by diffusion-weighted MRI (DW-MRI), cell proliferation via fluorodeoxythymidine PET (FLT-PET), apoptosis via 99mTc-Annexin-V SPECT and hypoxia imaging via fluoromisonidazole PET (FMISO-PET) and these were used to simulate tumor growth and treatment response. We explored three tumor models: modified logistic model with and without treatment, and the proliferation-apoptosis-migration model.
In general, the modified logistic model (equation (2)) leads to the development of three different layers of tumor cells: a proliferating rim, a quiescent region just beneath the rim and a necrotic core. With this model, as the tumor develops, most of the tumor begins to consist of a necrotic region. This is not unlike what is seen in some real tumors (Barth and Kaur 2009). As is expected, the tumor morphology and rate of growth are determined by the initial condition maps; in particular, the greater the proliferation rate, the greater the resulting necrotic region of the tumor and the greater the overall size of the tumor. This is not unexpected, but the fundamental point is that different initial conditions lead to different tumor morphology; thus modeling individual tumor growth become possible using readily available imaging data.
When treatment is added to the modified logistic model (equation (13)), the effects of effective and ineffective cytoxic therapies can be studied. For example, tumors that are highly proliferating (increased proliferation rate) are more likely to be resistant to drug treatment than tumors that have a lower proliferation rate. The timing of the drug delivery is also important as increasing the duration between doses can lead to tumors that are less responsive to drugs. Again, these findings are not surprising, but it is encouraging that models can be built and driven by imaging data and they can generate results that are not unreasonable.
The last model considered was new and built from the ‘ground up’ strictly in terms of data that can be supplied by currently available imaging studies. The PAM model incorporates the effect of proliferation, apoptosis and tumor cell migration and leads to a more aggressive tumor phenotype than the modified logistic model. The change of the tumor shape with time is influenced by the pO2 gradient in the local tissue volume of the tumor and the model assumes that tumor cells are driven into regions of higher pO2. Additionally, since the quiescent cells are also allowed to move to a region with more oxygenation, they can become proliferating cells and therefore increase the number of tumor cells with time. While selecting hypoxia as a driver of tumor cell migration is not without a basis, we acknowledge that the picture is much more complicated than what we have assumed in this implementation of the model. The model can however be amended to fit several other phenomena that are also potential drivers of cell migration. If these phenomena are able to be imaged in vivo and non-invasively, then such data could be incorporated into an updated version of equation (14).
Overall, the models explored here can simulate the effects of different tumor characteristics as a function of space and time. As the tumor models lead to different tumor morphology depending on the parameters maps that are used to initialize the simulations, the models are subject specific and as such changes in the initial condition maps will produce tumors with different shapes and volumes.
The analyses presented here represent an initial step, and will most likely need to be adjusted to increase the models’ biological relevancy. A fundamental major limitation is that these models do not explicitly incorporate vascularity, which is a determining factor in tumor development (Fukumura and Jain 2008). However there are a number of imaging methods that can report on tumor vascularity (e.g. dynamic contrast-enhanced MRI (Yankeelov et al 2007, Yang et al 2003) or angiography (van Vliet et al 2005) so the possibility of incorporating such data into these (or other) models is possible. The practical utility of the approaches described above remains to be proven. In particular, the models require data input from nuclear imaging technologies (PET and SPECT) that have relatively poor spatial resolution (approximately a few millimeters (Frangioni, 2008)) and therefore the tumor has to be of an appreciable size before this approach can be tested. While this may not be a significant limitation for clinical data, it may be so for small animal imaging. Conversely, the multi-modality techniques are certainly more practical for small animal imaging, but advances in hybrid technology suggest that several may be translated to clinical application. Other limitations include uncertainty in the relationship between the different imaging parameters and those presented in these (or other) models.
Another limitation is the somewhat large number of assumptions that still need to be made to perform the modeling. In particular, calibrations between imaging data and the parameters actually imported into the models; these include the mapping between a voxel’s ADC value and the absolute cell number in that voxel, or a voxel’s KFLT value and the actual tumor cell proliferation rate. Additionally, knowledge of the concentration and distribution of the cytoxic drug under consideration is required to run the modified logistic model with treatment effects. The models also do not include any anatomical barriers or preferred directions of tumor growth. For example, in the brain it is known that tumor cells may preferentially migrate along fiber tracks (Bondiau et al 2008). In other regions, interstitial pressure gradients or other mechanical barriers may force tumors to grow in certain directions. Such details are not included in these early models, though they can be incorporated. While several assumptions are currently required, we do not feel that they are fundamental barriers to the approach outlined.
The approach described is promising because the models can make predictions based on individual parameter values obtained from individual patients/subjects. Future efforts must include experiments to validate or refute the models. For example, a series of imaging measurements could be obtained just after treatment and used to drive one of the models above to predict how the tumor will look at some time point in the future. Then when the patient/subject reaches that time point, another set of imaging data can be obtained to test the predictions of the model(s). In this way each individual patient/subject would generate their own hypotheses which could then be tested. Such an approach is not currently available to mathematical models of tumor growth that are driven by parameters obtained invasively or from population averages.
It is clear that biomedical imaging and mathematical modeling can, independently, yield fundamental insights into cancer biology. When the two are combined the possibility of performing clinically relevant modeling of tumor growth and treatment response is compelling. In this contribution we presented evidence that it is possible to reparameterize existing tumor growth models (the logistic equation with and without the effects of treatment) in terms of readily available imaging data and that it can yield realistic tumor cell growth and distribution both with and without treatment. We further introduced a novel model that is described entirely in terms of imaging data and showed that it too can predict realistic growth and distribution of tumor cells. Further studies will include the experimental counterparts to these simulations to test the predictions of these models.
We thank the National Institutes of Health for funding through NIBIB 1K25EB005936, NIBIB 5 R01 EB000214 and NCI P50 CA128323 (In vivo Cellular and Molecular Imaging of Cancer program). We thank Drs Michael Miga, PhD, James Nutaro, PhD, Natasha Deane, PhD, and Michelle Martin, PhD, for many helpful discussions.