Home | About | Journals | Submit | Contact Us | Français |

**|**PLoS Comput Biol**|**v.9(3); 2013 March**|**PMC3610621

Formats

Article sections

- Abstract
- Author Summary
- Introduction
- Materials and Methods
- Results
- Discussion
- Supporting Information
- References

Authors

Related links

PLoS Comput Biol. 2013 March; 9(3): e1003008.

Published online 2013 March 28. doi: 10.1371/journal.pcbi.1003008

PMCID: PMC3610621

Hermann B. Frieboes,^{1,}^{2,}^{3,}^{*} Bryan R. Smith,^{4} Yao-Li Chuang,^{3} Ken Ito,^{4} Allison M. Roettgers,^{1} Sanjiv S. Gambhir,^{4,}^{5} and Vittorio Cristini^{3,}^{6}

Mark S. Alber, Editor^{}

University of Notre Dame, United States of America

* E-mail: hbfrie01/at/louisville.edu

The authors have declared that no competing interests exist.

Conceived and designed the experiments: HBF BRS YLC KI SSG VC. Performed the experiments: BRS YLC KI. Analyzed the data: HBF BRS KI AMR. Contributed reagents/materials/analysis tools: SSG VC. Wrote the paper: HBF BRS YLC KI SSG VC.

Received 2012 July 5; Accepted 2013 February 13.

This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

This article has been cited by other articles in PMC.

Non-Hodgkin's lymphoma is a disseminated, highly malignant cancer, with resistance to drug treatment based on molecular- and tissue-scale characteristics that are intricately linked. A critical element of molecular resistance has been traced to the loss of functionality in proteins such as the tumor suppressor *p53*. We investigate the tissue-scale physiologic effects of this loss by integrating *in vivo* and immunohistological data with computational modeling to study the spatiotemporal physical dynamics of lymphoma growth. We compare between drug-sensitive *Eμ-myc Arf-/-* and drug-resistant *Eμ-myc p53-/-* lymphoma cell tumors grown in live mice. Initial values for the model parameters are obtained in part by extracting values from the cellular-scale from whole-tumor histological staining of the tumor-infiltrated inguinal lymph node *in vivo*. We compare model-predicted tumor growth with that observed from intravital microscopy and macroscopic imaging *in vivo*, finding that the model is able to accurately predict lymphoma growth. A critical physical mechanism underlying drug-resistant phenotypes may be that the *Eμ-myc p53-/-* cells seem to pack more closely within the tumor than the *Eμ-myc Arf-/-* cells, thus possibly exacerbating diffusion gradients of oxygen, leading to cell quiescence and hence resistance to cell-cycle specific drugs. Tighter cell packing could also maintain steeper gradients of drug and lead to insufficient toxicity. The transport phenomena within the lymphoma may thus contribute in nontrivial, complex ways to the difference in drug sensitivity between *Eμ-myc Arf-/-* and *Eμ-myc p53-/-* tumors, beyond what might be solely expected from loss of functionality at the molecular scale. We conclude that computational modeling tightly integrated with experimental data gives insight into the dynamics of Non-Hodgkin's lymphoma and provides a platform to generate confirmable predictions of tumor growth.

Non-Hodgkin's lymphoma is a cancer that develops from white blood cells called lymphocytes in the immune system, whose role is to fight disease throughout the body. This cancer can spread throughout the whole body and be very lethal – in the US, one third of patients will die from this disease within five years of diagnosis. Chemotherapy is a usual treatment for lymphoma, but the cancer can become highly resistant to it. One reason is that a critical gene called *p53* can become mutated and help the cancer to survive. In this work we investigate how cells with this mutation affect the cancer growth by performing experiments in mice and using a computer model. By inputting the model parameters based on data from the experiments, we are able to accurately predict the growth of the tumor as compared to tumor measurements in living mice. We conclude that computational modeling integrated with experimental data gives insight into the dynamics of Non-Hodgkin's lymphoma, and provides a platform to generate confirmable predictions of tumor growth.

Monoclonal antibodies and small molecule inhibitors of intracellular targets are being developed alongside a host of anti-non-Hodgkin's lymphoma therapeutic options [1]. Yet the tumor tissue-scale effects from these molecular-scale manipulations are not well-understood. With the ultimate goal to more rationally optimize lymphoma treatment, we integrate pre-clinical *in vivo* observations of lymphoma growth with computational modeling to create a platform that could lead to optimized therapy. As a first step towards this goal, we develop the capability for simulation in order to gain insight into the tissue-scale effect of molecular-scale mechanisms that drive lymphoma growth. We use the modeling to study these mechanisms and their association to cell proliferation, death, and physical transport barriers within the tumor tissue.

Tumor growth and treatment response have been modeled using mathematics and numerical simulation for the past several decades (see recent reviews [2]–[9]). Models are usually either discrete or continuum depending on how the tumor tissue is represented. *Discrete models* represent individual cells according to a specific set of bio-physical and -chemical rules, which is particularly useful for studying carcinogenesis, natural selection, genetic instability, and cell-cell and cell-microenvironment interaction (see reviews by [10]–[20]). *Continuum models* treat tumors as a collection of tissue, applying principles from continuum mechanics to describe cancer-related variables (e.g., cell volume fractions and concentrations of oxygen and nutrients) as continuous fields by means of partial differential and integro-differential equations [2]. A third modeling approach employs a *hybrid* combination of both continuum and discrete representations of tumor cells and microenvironment components, aiming to develop multiscale models where the discrete scale can be directly fitted to molecular and cell-scale data and then upscaled to inform the phenomenological parameters at the continuum scale (see recent work by [21]–[23]).

There is a paucity of mathematical oncology work applied to the study of non-Hodgkin's lymphoma, with some notable exceptions providing insight into the role of the tumor microenvironment heterogeneity in the treatment response [24], [25] and the disease origin [26]. Like many other cancers (solid tumors), two critical tissue-scale effects in lymphoma are hypoxia and angiogenesis, as observed in our studies and other work [27]. Supporting previous qualitative observations of physiological resistance, mathematical modeling and computational simulation have shown that the diffusion barrier alone can result in poor tumor response to chemotherapy due to diminished delivery of drug, oxygen, and cell nutrients [28], [29]. Local depletion of oxygen and cell nutrients may further promote survival to cell cycle-specific drugs through cell quiescence.

In order to study these effects in lymphoma, we implement an integrated computational/experimental approach to quantitatively link the processes from the cell scale to the tumor tissue-scale behavior in order to gain insight into their cause and progression in time. We extend a version of our 3D continuum model [30]–[32], building upon extensive mathematical oncology work [2], [3], [33]–[35], and calibrate both parameters and equations, i.e., functional relationships that are not conservation laws, from detailed experimental data to produce a virtual lymphoma. We obtain the experimental data by very fine sectioning of both drug-sensitive and -resistant lymphomas, thus visualizing molecular, cellular, and tissue-scale parameter information across the whole tumor geometry. We further develop the protocols for calibration of parameters by building on recent work based on patient histopathology [36], [37]. We also use the data to derive the relationships between model parameters for apoptosis, proliferation, and vasculature. We verify the model results at the tumor-scale through tissue-scale observations *in vivo* of tumor size, morphology, and vasculature using intravital microscopy and macroscopic imaging of the inguinal lymph node. We note that comparison of model results to experimental data has been done to various extents for different cancers (see reviews above); here, we perform a tissue-scale comparison after extensive calibration of cell-scale parameters in order to validate the model results. We undertake simulations to study how the growth of drug-resistant Non-Hodgkin's lymphoma may be governed by the cellular phenotype, and use this information to better elucidate the links between physical drug resistance and molecular-scale phenotype by experimental and computational comparison to drug-sensitive tumors.

This process yields a lymphoma simulator as an initial step to study detailed tumor progression and provide further insight into drug resistance, and, ultimately, may provide a tool to design better personalized treatments for Non-Hodgkin's lymphoma. Since the cell-scale measurements used for calibration are different from those at the tissue-scale used for verification, this methodology enables the model to bridge from the cell to the tumor scale to calculate tumor growth and hypothesize associated mechanisms predictively, i.e., without resorting to fitting to the experimental data. This process quantitatively links the cellular phenotype to the tumor tissue-scale behavior, and may serve to highlight the importance of physical heterogeneity and interactions in the tumor microenvironment when evaluating chemotherapeutic agents in addition to consideration of chemo-protective effects such as cell-specific phenotypic properties and cell-cell and cell-ECM adhesion [38].

We choose an *Eμ-myc* murine orthotopic lymphoma experimental model because of its similarity to human Non-Hodgkin's Lymphoma [39], and select five parameters to measure based on their importance to lymphoma progression: viability, hypoxia, vascularization, proliferation, and apoptosis. In order to investigate the role of physical heterogeneity in the development of drug resistance, including the impediment of transport barriers, we focus on two types of lymphoma cells: *Eμ-myc Arf-/-* cells (Doxorubicin (DOX) and Cyclophosphamide (CTX) sensitive, with IC_{50}=3.5 nM and 16.0 µM, respectively; the IC_{50} is the amount of drug needed to kill 50% of a cell population), and *Eμ-myc p53-/-* cells (DOX and CTX resistant: IC_{50}=46.2 nM and 75.8 µM, respectively). The *Eμ-myc* transgenic mouse model expresses the *Myc* oncogene in the B cell compartment, resulting in mice with transplantable B cell lymphomas. We chose this *in vivo* model because it captures genetic and pathological features of the human disease and, given the appropriate genetic mutation, drug-resistant and drug-sensitive tumors can be directly compared [39],[40].

*Eμ-myc/Arf-/-* and *Eμ-myc/p53-/-* lymphoma cells, which harbor loss-of-function regions in the *Arf* and *p53* genes respectively, were previously derived by intercrossing *Eμ-myc* transgenic mice with *Arf*-null and *p53*-null mice, all in the C57BL/6 background as described previously [39]. *Eμ-myc/Arf-/-* lymphoma cells and *Eμ-myc/p53-/-* lymphoma cells were cultured in 45% Dulbecco's modified Eagle medium (DMEM) and 45% Iscove's Modified Dulbecco's Medium (IMDM) with 10% fetal bovine serum (FBS) and 1% penicillin G-streptomycin onto the feeder cells – Mouse Embryonic Fibroblasts (MEFs).

C57BL/6 mice were obtained from Charles River Laboratories (Wilmington, Massachusetts). All animal studies were approved by The Stanford University Institutional Animal Care and Use Committee. Lymphoma cells (1×10^{6}) *Eμ-myc/Arf-/-* and *Eμ-myc/p53-/-* were diluted with 200 µl of PBS and injected intravenously via the tail vein as described previously [39]. The intravital microscopy and macroscopic tumor observations were obtained for at least n=4 mice per tumor group.

We isolated both *Eμ-myc/Arf-/-* and *Eμ-myc/p53-/-* driven tumors at day 21 after tail-vein injection of lymphoma cells. Typical murine lymphomas were observed to range from about 4 to 6 mm in diameter prior to fixation. Lymph node tissues were fixed and paraffin-embedded. The tissues were used for immunohistochemical (IHC) identification of cell viability (H&E staining), hypoxia (HIF-1α), vascularization (CD31), proliferation (Ki-67), and apoptosis (Caspase-3). Five 2-µm thick sections were cut 5 µm apart from each other in order to stain for these markers (**Figure 1**). A total of five sets (S1 through S5) of five stained sections each was collected every 100 µm along the lymphoma, in order to section and stain the entire tumor for sequential microscopic scanning of the stained sections. Sections S1 and S5 were at the tumor top and bottom, respectively, while the other sections were towards the center with S3 being in the middle. Note that due to tissue processing and dehydration, the tumors as cut were smaller than measured when removed from the animal. All the sections were de-paraffinized and rehydrated in PBS. Then the sections in each set were incubated at 4°C with the primary antibody overnight: rabbit anti-mouse HIF-1 antibody (Abcam, Santa Cruz, CA), rabbit anti-mouse Ki-67 antibody (Labvision, Fremont, CA), rabbit anti-mouse Caspase-3 antibody (Cell Signaling Technology, Beverly, CA), and rat anti-mouse CD31 antibody (BD Pharmingen, San Diego, CA), and incubated for 1 hour at room temperature with a peroxidase-conjugated secondary antibody. The samples were fully scanned and stitched together using a digital pathology BioImagene instrument (Ventana Medical Systems, Tucson AZ) at ×20 magnification.

The model treats tissue as a mixture of various cell species, water, and ECM; each component is subject to physical conservation laws described by diffusion-taxis-reaction equations (see below). Briefly, the tissue microstructure is modeled through the proper choice of parameter values and through biologically-justified functional relationships between these parameters, e.g., cellular transitions from quiescence to proliferation depend upon oxygen concentration [41]. The model simulates non-symmetric tumor evolution in 2D and 3D, and dynamically couples heterogeneous growth, vascularization, and tissue biomechanics (**Figure 2**). In [36] we calibrated models using cell-scale data to predict tissue scale parameters such as size and growth rate. These models are predictive because they are not calibrated with the same data used for model validation, which avoids data fitting. While in [36] we focused on the final predicted tumor sizes, here we focus on the growth rate as an essential first step; in follow-up work, we will evaluate the complex problem of drug response. Our approach to constrain the computational model involves both cell- and tumor-scale approaches as described in **Figure 3**.

We approximate the healthy lymph node as a sphere to represent the experiments in the mouse model (**Figure 4**). To simulate node expansion and deformation of surrounding tissue to accommodate the growing tumor, as a first step we delineate the tumor boundary by decreasing the value of the cell mobility parameter beyond the sphere diameter (see below). For the multigrid algorithm, we pick a computational domain that is a 6.4 mm× 6.4 mm× 6.4 mm box, with finest mesh grid size=100 microns; this grid size provides adequate resolution to resolve the tumor boundaries without incurring excessive computational cost.

We assume that the tumor is a mixture of cells, interstitial fluid, and extracellular matrix (ECM). The temporal rate of change in viable and dead tumor tissue at any location within the tumor equals the amount of mass that is pushed, transported, and pulled due to cell motion, adhesion, and tissue pressure, plus the net result of production and destruction of mass due to cell proliferation and death:

(1)

The rate of change in the volume fraction *ρ _{i}* of cell species

Tumor angiogenesis is driven by excessive accumulation of cancerous cells, leading to a chronic under-supply of oxygen and cell nutrients (generically here labeled “nutrients”) in tumor regions farther removed from pre-existing vessels [44]. Hypoxic cells in lymphoma release a net balance of pro-angiogenic factors such as VEGF-A, bFGF, PDGF and VEGF-C, which promote neo-vascularization mainly through sprouting angiogenesis of mature resident endothelial cells and, to a lesser extent, through vasculogenesis from recruitment of bone marrow-derived progenitor cells [45]. Accordingly, the model incorporates angiogenesis into the lymphoma by coupling with a multiscale representation of tumor vessel growth, branching, and anastomosis based on earlier work [46]–[48] (further details in **Text S1**).

The vasculature releases oxygen and nutrients *n* that diffuse through the tissue and are uptaken by cells during metabolism, while tumor cells secrete VEGF (*n _{V}*) in response to hypoxia [32]. The oxygen and nutrients are non-dimensionalized by the maximum inside vessels, hence their levels are ≤1 and are assumed to be stationary. The transport can be described as:

(2)

where *D _{n}* and are diffusion constants (1×10

The tumor species viable (*V*) volume fraction *ρ _{V}* is assumed to increase through proliferation and decrease through apoptosis and necrosis. We assume that normal host cells (

(3)

where *λ*_{M,i}, *λ*_{A,i}, and *λ*_{N,i} are mitosis, apoptosis, and necrosis rates, *λ*_{D} is the cell degradation rate (varies due to the differences between apoptosis and necrosis), and *H*(*x*) is the Heaviside “switch” function.

The movement of a cell species is determined by the balance of proliferation-generated oncotic pressure, cell-cell and cell-ECM adhesion, as well as chemotaxis (due to substrate gradients), and haptotaxis (due to gradients in the ECM density). We model the motion of cells and interstitial fluid through the ECM as a viscous, inertialess flow through a porous medium. Therefore, no distinction between interstitial fluid hydrostatic pressure and mechanical pressure due to cell-cell interactions is made. Cell velocity is a function of cell mobility and tissue oncotic (solid) pressure (Darcy's law); cell-cell adhesion is modeled using an energy approach from continuum thermodynamics (see **Text S1**). For simplicity, the interstitial fluid is modeled as moving freely through the ECM (i.e., at a faster time scale than the cells).

(4)

The variational derivative δ*E*/δ*ρ _{i}* of the cell-cell interaction potential, combined with the remaining contributions to the flux

We used the IHC staining to estimate the number and spatial localization of cells that were viable (from H&E), proliferating (from Ki-67), apoptotic (from Caspase-3), hypoxic (from HIF-1α), and with vascular endothelial characteristics (from CD31). These estimates were calculated for both *Eμ-myc Arf-/-* and *Eμ-myc p53-/-* cells for each set of five sections obtained every 100 µm across the lymphoma (**Figures 5**
**and**
**6****6**).

A comparison of viable *Eμ-myc p53-/-* to *Eμ-myc Arf-/-* cells along the lymphoma (**Figure 5**) indicates that the viability is higher for the drug-resistant tumors in the middle of the tumor (Section S3) compared to the drug-sensitive tumors, with a corresponding statistically significant increase in cell density (p=0.024; Student's t-test with α=0.05). In contrast to the *Eμ-myc p53-/-* tumors, the *Eμ-myc Arf-/-* seemed to be more dense in the peripheral regions (p=0.002 on one end (Section S1) and p=0.009 on the other end (Section S5)), whereas they were about the same for both tumor types in the intermediate sections S2 and S4. Tumors with drug-resistant cells have a 4-fold increase in endothelial cells in the core of the tumor (Section S3) compared to drug-sensitive tumors (**Figure 6A**). Hypoxia is higher in the peripheral regions for the *Eμ-myc p53-/-* (**Figure 6B**) even though for both tumor types the peripheral regions seem to be equally vascularized (based on the endothelial cell density). This could be due to the vasculature on the periphery not being fully functional, with a potential difference in vascular function between the two tumor types leading to a more hypoxic phenotype for the *Eμ-myc p53-/-*. Although the core proportionally holds almost twice the number of proliferating cells for the drug-resistant tumors as compared to the drug-sensitive case (**Figure 6C**), a correlation between proliferation and vascularization/hypoxia is precluded. Interestingly, the number of apoptotic cells is consistently higher for *Eμ-myc p53-/-* (**Figure 6D**), suggesting non-hypoxia driven apoptosis for these tumors.

By analyzing each IHC section longitudinally along the tumor, a range of baseline values can be calculated from the experimental data for key model parameters (**Table S1**), inspired by recent methods in mathematical pathology [36]: cell viability, necrosis, and spatial distribution pattern (from H&E), cell proliferation (from Ki-67), cell apoptosis (from Caspase-3), oxygen diffusion distance (from HIF-1α), and blood vessel density (from CD31). These values are obtained for both *Eμ-myc Arf-/-* and *Eμ-myc p53-/-* tumors for each of the five sections obtained longitudinally along the tumor, with values sampled from the middle (core) and the edge (periphery) of each section. The measured values are not resolved in space but averaged over each section, thus yielding information averaged over space. The periphery was defined as the region approximately within 200 µm of the tumor boundary.

**Figure S1** shows an example of this calibration process for proliferation at the periphery and middle from two histology sections in the center of the tumor (Section S3). Taking an average proliferation cycle of 20 hours that we observed for the lymphoma cells in culture, the proliferation calculation in units of day^{−1} is λ_{M}**<n>*=[(stained/(stained+unstained))/20 hours/prolif.] * 24 hours/day. The average nutrient *<n>* indicates that this proliferation rate depends on the model diffusion of cell substrates such as glucose and oxygen in the 3D space (**Eq. 2**). Similarly, since the apoptosis cycle was detectable up to 5 hours, the apoptosis calculation in units of day^{−1} is λ_{A}=[(stained/(stained+unstained))/5 hours/apoptosis] * 24 hours/day.

We calculate the average nutrient from the blood vessel density by assuming a uniform nutrient delivery rate from the blood to the tissue adjacent to the vessels (**Eq. 2**). Estimating blood vessel area versus surrounding tissue provides a measure of the magnitude of cell substrates transferred into the tumor. Thus, we calculate the fraction of cells supported per endothelial cell in a unit volume to be (number unstained/(number stained+unstained))^{3/2}. When the viable cell fraction in the simulations matches what is directly observed from microscopy, this implies that the vascular and nutrient distributions have been correctly represented in the model (**Figure 3A****, middle**). Similarly, we calculate the hypoxic cell fraction per unit volume as (number stained/(number stained+unstained))^{3/2}.

The node is represented by the computational model initially as a spherical capsule in 3D with a membrane boundary separating it from the surrounding tissue (**Figure 4**) (see **Text S1**). Lymphoma cells are assumed to enter the lymph node through the afferent lymph vessels. As they accumulate in the node during tumor progression in time, they compete for cell substrates such as oxygen and nutrients with the normal lymphocytes. These substrates are assumed to diffuse radially outward toward the node periphery from the pre-existing vasculature, situated mainly in the core of the node (see **Figure 4A** and **Figure 3B****, left**, at the intersection of three large blood vessels). Once a tumor has begun to form in the core of the node, this diffusion process presents a transport barrier for oxygen and nutrients to the lymphoma cells incoming through the afferent vessels into the node.

We investigate the effect of initially available oxygen and cell substrates needed for cell proliferation, since lymphoma growth is hypothesized to depend on access to these through the vasculature. Preliminary calculations suggested that the initially available nutrient level has a significant effect on the growth phase of the tumor but not on its terminal size, which according to a theoretical analysis of the model (**Text S1**) depends mainly on the ratio of apoptosis to proliferation [51]. A further investigation revealed that the initial guess of parameter values results in a mismatch between the ratio of hypoxic cells and the average apoptosis rate: where the range of hypoxic ratio matches the experiments, the apoptosis rate range in the model is too low.

Accordingly, we calibrated the cell necrosis rate so that the key parameter values remain invariant when the initial nutrient is set to a threshold of 0.5. With this set of parameters, a necrosis rate from 5 to 7 (non-dimensional units) would satisfy the experimentally observed ranges of both the hypoxic fractions and the average apoptosis rate (**Figure S2A–B**). We then varied the initial nutrient threshold while maintaining the necrosis rate invariant to confirm that the fraction of hypoxic cells and average apoptosis rate would remain within the experimentally observed range of values (**Figure S2C–D**). This calibration suggests that the initially available nutrient still affects the growth phase of lymphoma. In this model, the lymphoma tumor and the lymph node greatly outgrow the original lymph node size, which we consistently observed *in vivo* in addition to the distortion of the lymph node geometry (we are currently implementing the Diffuse Domain Method [52] to better represent this geometry).

After using the IHC data to perform a cell-scale calibration of the lymphoma model, we verify the simulated tissue-scale lymphoma size from *in vivo* macroscopic observations and intravital imaging at the tissue scale. Recently, it has been discovered with bioluminescence imaging by Gambhir and co-workers that lymphoma cells coming from the spleen and bone marrow seed the inguinal lymph node around Day 9 *in vivo*
[53]. Using this seeding as the initial condition for the simulations, the model predicts the tumor diameter to be ~5.2±0.5 mm by Day 21 (**Figure 7**). This figure also shows the gross tumor size from our caliper measurements in time, indicating that the model-predicted tumor diameter for the maximum possible value of initial nutrient falls within the range of the measurements *in vivo* (the experiments show that there is no statistical difference in the tumor growth between the two cell types, **Figure 3B****, right**). The model simulations are based on an oxygen diffusion distance from the vessels estimated to be directly proportional to the distance at which hypoxia is detected away from blood vessels, measured experimentally from the HIF-1α staining to be 80±20 µm. The variation in this measurement leads to the variation in the simulated diameter. If the lymphoma is begun at sites within the lymph node other than the center (**Figure 4**), similar growth curves are computationally obtained as the whole node volume is eventually taken over by the proliferating tumor cells (results not shown). We note that since there is a distributed source of vessels in the tumor, the proliferation is relatively weakly sensitive to additional outside sources.

The tumor growth from the model calibrated from the cell-scale can be validated through theoretical analysis of the model based on previous mathematical and computational work [51], [54]–[56] (see **Text S1**). Assuming that the lymph node geometry is approximated by a 3D sphere, the model can be used to predict the tumor radius in time based on the ratio *A* of the rates of apoptosis to proliferation calculated from the experimental IHC data. The average ratio *A*=λ_{A}/λ_{M} ~0.4 for both drug-sensitive and drug-resistant cells. In comparison with the simulations based on the cell-scale calibration, this analysis predicts that the tumor would reach a diameter of ~6 mm. Both the theoretical analysis and the tumor growth obtained through the simulations agree with the similar diameters observed experimentally *in vivo* (~5 to 6 mm) (**Figure 7**).

In the model, simulations of the vasculature were qualitatively compared to independent intravital microscopy observations *in vivo* of a *Eμ-myc p53-/-* tumor in the same animal over time (**Figure 8**). The density of simulated viable tumor tissue (**Figure 3A****, right**) as a function of the vascularization at day 21 qualitatively matches the density of the tissue observed experimentally (fraction of simulated viable cells in the 2D plane, >90% per mm^{2} in inset in **Figure 8**, vs. the average fraction of viable cells measured from H&E staining, 87%±6% per mm^{2}), indicating that the overall vasculature function was modeled properly. The density of simulated endothelial tissue is also highest in the tumor core, as observed from histology. The increase in the lymphoma cell population disturbs the homogeneous distribution of cell substrates (such as oxygen and cell nutrients), leading to diffusion gradients of these substances that in turn affect the lymphoma cell viability. If the cell viability is established heterogeneously within the tumor, e.g., as observed experimentally in IHC with the *Eμ-myc Arf-/-* cells near the tumor periphery, the model predicts that the diffusion gradients would not be as pronounced. If the cell viability is higher near the center of the tumor, which is observed in IHC with the *Eμ-myc p53-/-* cells (**Figure 5**), then the gradients are predicted to be steeper and more uniform [28].

We integrate *in vivo* lymphoma data with computational modeling to develop a basic model of Non-Hodgkin's lymphoma. Through this work we seek a deeper quantitative understanding of the dynamics of lymphoma growth in the inguinal lymph node and associated physical transport barriers to effective treatment. We obtain histology data by very fine sectioning across whole lymph node tumors, thus providing detailed three-dimensional lymphoma information. We develop a computational model that is calibrated from these cell-scale data and show that the model can independently predict the tissue-scale tumor size observed *in vivo* without fitting to the data. We further show that this approach can shed insight into the tumor progression within the node, particularly regarding the physical reasons why some tumors might be resistant to drug treatment – a critical consideration when attempting to quantify and predict the treatment response. We envision that the modeling and functional relationships derived in this study could contribute with further development to patient-specific predictors of lymphoma growth and drug response.

Although the number of mice used for the experimental *in vivo* validation is limited, the model results are consistent with previous work. For example, a well-studied mechanism of physiological resistance is the dependence of cancer cell sensitivity to many chemotherapeutic agents on the proliferative state of the cell [28]. This physical mechanism is likely important in the difference in drug-sensitivity between the tumors formed from the two cell lines and will be explored in further studies. We found that the *Eμ-myc Arf-/-* cells tend to congregate at the periphery of the tumor (**Figure 5**), even though there are vessels in the interior of the tumor. This suggests the hypothesis that the more drug-sensitive *Eμ-myc Arf-/-* cells maintain better oxygenation at the expense of higher drug sensitivity by growing less compactly in the interior of the tumor – where there would be stronger competition for oxygen and cell nutrients – whereas the *Eμ-myc p53-/-* lymphoma cells may enhance their survival by closer packing in the core of the tumor. Cell packing density may present a barrier to effective drug penetration [57], which we have also modeled previously [28]. Closer packing could further increase the number of cells that would be quiescent due to depletion of oxygen and nutrients, as we specify in the model (**Materials and Methods**) and as we have simulated in previous work [28]. However, the proportion of chemoresistance inherent with *Eμ-myc p53-/-* that can be attributed to resistance at the genetic level compared to what can be attributed to suboptimal drug delivery and quiescence is unclear. In follow-up work we plan to measure drug amounts near various cells in order to begin answering this question, and to perform sensitivity analyses of the IC50 of each cell line with the computational model. This would provide a (model-generated) measure of how much of an effect suboptimal delivery could be attributed to resistance as compared to genetic effects (as measured by IC50).

Lymphoma cells are known to retain cell-cell adhesion, with strength associated with the lymphoma's originating cell type (B- or T-cell) [58]. Mechanisms of cell packing related to drug resistance may include weaker cell adhesion in *Eμ-myc Arf-/-* than in *Eμ-myc p53-/-* leading to higher cell density as well as a denser extra-cellular matrix in the latter [57]. Loss of ARF has been linked to increased cancer cell migration and invasion, and hence weaker cell-cell adhesion [59], associated with the binding of ARF to the transcriptional corepressor CtBP2 and promoting CtBP2 degradation [60]–[62].

Perhaps surprisingly, the experimental data indicate minimal presence of hypoxia within the tumor (**Figure 6B**). This may be due to the fact that lymphoma cells may associate with other cells including stromal cells in the tumor, and the consequent cytokine stimulation (e.g., IL-7) may also trigger proliferation [63]. We note that the oxygen diffusion length estimate is subject to variation, as calculated to be directly proportional to the hypoxic distances observed from the IHC; this may be improved by directly measuring the diffusing substances, e.g., oxygen. The simulated elastic tumor boundary may also introduce some variation into the size calculation. Nevertheless, even taking these variations into account, the model-calculated average ratio of apoptosis to proliferation, established from cell-scale measurements, implies that the tumor sizes fall within the range of the sizes estimated from the diameter measured with calipers *in vivo*. The hypothesis we test with the model by successful comparison to the experimental data is that the growth and eventual slowdown of these tumors is the balance of proliferation and death, which we have also previously observed for ductal carcinoma in situ [38]. Experimental evidence using bioluminescence imaging of living mice [53] demonstrates that lymphoma cells seed the tumor in the inguinal lymph node from other sites (e.g., spleen and bone marrow) in the mouse body at earlier times during the tumor growth. The model results are robust, however, because the tumor size by Day 21 predicted by the theory is independent of the earlier times; any influx of cells only provides an initial (transient) condition.

The staining also shows that apoptosis seems highest for drug-sensitive cells at the periphery of the tumor (Sections S1 and S5) compared to the center (Section S3) (both p-values=0.04 using a Student's t-test with α=0.05), and for drug-resistant cells it is highest in the more central regions (**Figure 6**). In accordance with biological observations [64], [65], [41], the model hypothesizes that increased hypoxia may lead to higher cell quiescence and hence drug resistance. In the experiments, angiogenesis is higher in the central regions, and is more pronounced for drug-resistant cells, suggesting that these cells are in a more angiogenic environment as a result of ongoing hypoxic stimulus. Higher tumor cell density around blood vessels suggests a functional relationship of cell viability as a function of nutrients, as we have implemented in the model (see **Materials and Methods**). However, apoptosis may not necessarily be driven solely by hypoxia, since lymphoma cells are known to have a cellular turnover rate that is on the order of days [66], [67]. We further note that angiogenesis is not necessarily triggered only by hypoxia. Lymphoma as well as stromal cells (such as tumor associated macrophages) may produce factors promoting angiogenesis (e.g., vascular endothelial growth factor or VEGF) under otherwise normoxic conditions.

The present work calibrates a computational model of lymphoma with experimental data from drug-sensitive and drug-resistant tumors. This data was derived from detailed IHC analysis of whole tumors, and validation of the model was performed via intravital microscopy measurements. The results suggest that differences in spatial localization of cells and vasculature, as well as in the transport phenomena in the tumor microenvironment may play a nontrivial role in the tumor behavior. This suggests that the genetic differences (*Eμ-myc Arf-/-* and *Eμ-myc p53-/-*) may provide a substantial compensation mechanism for these phenomena at the tissue scale in addition to the molecular as it relates to their drug resistance. We plan to verify this hypothesis in the future by assessing model predictions for therapeutic response of drug-sensitive and drug-resistant tumors in terms of cellular parameters such as proliferation, apoptosis, and hypoxia via both IHC and intravital microscopy.

**Example of calibration process of model parameters from the Ki-67 IHC data.** The proliferation parameter is calculated for both *Eμ-myc Arf-/-* (drug-sensitive) and *Eμ-myc p53-/-* (drug-resistant) lymphoma cells. This sample (from Set S3 in the center of the tumor) shows measurements obtained at the edge (periphery) and middle (center) of the section. Positive staining shown in the panels A–D is converted to red and negative staining to green in panels E–H to obtain a quantitative measure of proliferative activity, as calculated in the text. Results are shown in bottom right insets in panels E–H.

(TIF)

Click here for additional data file.^{(8.3M, tif)}

**Determination of optimal necrotic rate threshold for cell viability.** The necrosis rate is varied while the initial nutrient threshold is fixed at 0.5 to determine a range for which both the hypoxic fractions (A) and average apoptosis rate (B) match what is observed experimentally, finding that this range is from 5 to 7 (non-dimensionalized). We then varied the initial nutrient threshold while maintaining the necrosis rate invariant to confirm that the fraction of hypoxic cells (C) and average apoptosis rate (D) would remain within the experimentally observed ranges.

(TIF)

Click here for additional data file.^{(1.2M, tif)}

**Range of key parameter values and corresponding baseline values for the computational model.** (M) values were calculated from the cell-scale immunohistochemistry data, (C) values were calibrated using these data, and (ND) are non-dimensionalized values.

(TIF)

Click here for additional data file.^{(923K, tif)}

We are grateful to John Lowengrub (Mathematics, UCI) for useful discussions and advice, and to Fang Jin (Pathology, UNM) for enhancements to the tumor angiogenesis model. We wish to thank the reviewers for their valuable contribution.

This work was supported in part by NCI PSOC MC-START U54CA143907 (VC, SSG, HF), NCI ICBP 1U54CA151668 (VC), NCI ICMIC P50CA114747 (SSG), NCI RO1 CA082214 (SSG), NCI CCNE-TR U54 CA119367 (SSG), CCNE-T U54 U54CA151459 (SSG), and Canary Foundation (SSG), and K99 CA160764 (BRS). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

1. Mahadevan D, Fisher RI (2011) Novel therapeutics for aggressive non-Hodgkin's lymphoma. J Clin Oncol 29: 1876–1884. [PMC free article] [PubMed]

2. Lowengrub JS, Frieboes HB, Jin F, Chuang YL, Li X, et al. (2010) Nonlinear modelling of cancer: bridging the gap between cells and tumours. Nonlinearity 23: R1–R9. [PMC free article] [PubMed]

3. Frieboes HB, Chaplain MA, Thompson AM, Bearer EL, Lowengrub JS, et al. (2011) Physical oncology: a bench-to-bedside quantitative and predictive approach. Cancer Res 71: 298–302. [PMC free article] [PubMed]

4. Byrne HM (2010) Dissecting cancer through mathematics: from the cell to the animal model. Nat Rev Cancer 10: 221–230. [PubMed]

5. Astanin S, Preziosi L (2007) Multiphase Models of Tumour Growth. In: Bellomo N, Chaplain M, DeAngelis E, editors. Selected Topics on Cancer Modelling: Genesis - Evolution - Immune Competition – Therapy. Boston: Birkhäuser. pp. 1–31.

6. Harpold HL, Alvord Jr EC, Swanson KR (2007) The evolution of mathematical modeling of glioma proliferation and invasion. J Neuropathol Exp Neurol 66: 1–9. [PubMed]

7. Anderson ARA, Quaranta V (2008) Integrative mathematical oncology. Nat Rev Cancer 8: 227–244. [PubMed]

8. Deisboeck TS, Zhang L, Yoon J, Costa J (2009) In silico cancer modeling: is it ready for prime time? Nature Clin Practice Oncol 6: 34–42. [PMC free article] [PubMed]

9. Ventura AC, Jackson TL, Merajver SD (2009) On the role of cell signaling models in cancer research. Cancer Res 69: 400–402. [PubMed]

10. Quaranta V, Weaver AM, Cummings PT, Anderson ARA (2005) Mathematical modeling of cancer: the future of prognosis and treatment. Clin Chim Acta 357: 173–179. [PubMed]

11. Nagy JD (2005) The ecology and evolutionary biology of cancer: A review of mathematical models of necrosis and tumor cell diversity. Math Biosci Eng 2: 381–418. [PubMed]

12. Abbott RG, Forrest S, Pienta KJ (2006) Simulating the hallmarks of cancer. Art Life 12: 617–34. [PubMed]

13. Byrne HM, Alarcón T, Owen MR, Webb SW, Maini PK (2006) Modeling aspects of cancer dynamics: A review. Phi Trans R Soc A 364: 1563–1578. [PubMed]

14. Fasano A, Bertuzzi A, Gandolfi A (2006) Mathematical modelling of tumour growth and treatment. In: Complex Systems in Biomedicine. Milan: Springer. Pp. 71–108.

15. Galle J, Aust G, Schaller G, Beyer T, Drasdo D (2006) Individual cell-based models of the spatial temporal organization of multicellular systems-achievements and limitations. Cytometry A 69: 704–710. [PubMed]

16. Drasdo D, Höhme S (2007) On the role of physics in the growth and pattern of multicellular systems: What we learn from individual-cell based models? J Stat Phys 128: 287–345.

17. Thorne BC, Bailey AM, Peirce SM (2007) Combining experiments with multi-cell agent-based modeling to study biological tissue patterning. Briefings Bioinformatics 8: 245–257. [PubMed]

18. Anderson ARA, Chaplain M, Rejniak K, Fozard J (2008) Single-cell based models in biology and medicine. Math Med Biol 25: 185–186.

19. Quaranta V, Rejniak K, Gerlee P, Anderson ARA (2008) Invasion emerges from cancer cell adaptation to competitive microenvironments: quantitative predictions from multiscale mathematical models. Sem Cancer Biol 18: 338–348. [PubMed]

20. Zhang L, Wang Z, Sagotsky JA, Deisboeck TS (2009) Multiscale agent-based cancer modeling. J Math Biol 58: 545–559. [PubMed]

21. Kim Y, Stolarska MA, Othmer HG (2007) A hybrid model for tumor spheroid growth in vitro: I. Theoretical development and early results. Math Methods Appl Sci 17: 1773–1798.

22. Stolarska MA, Kim Y, Othmer HG (2009) Multiscale models of cell and tissue dynamics. Phil Trans R Soc A 367: 3525–3553. [PMC free article] [PubMed]

23. Bearer EL, Lowengrub JS, Frieboes HB, Chuang YL, Jin F, et al. (2009) Multiparameter computational modeling of tumor invasion. Cancer Res 69: 4493–4501. [PMC free article] [PubMed]

24. Ribba B, Marron K, Agur Z, Alarcón T, Maini PK (2005) A mathematical model of Doxorubicin treatment efficacy for non-Hodgkin's lymphoma: investigation of the current protocol through theoretical modelling results. Bull Math Biol 67: 79–99. [PubMed]

25. Alarcón T, Marches R, Page KM (2006) Mathematical models of the fate of lymphoma B cells after antigen receptor ligation with specific antibodies. J Theor Biol 240: 54–71. [PubMed]

26. Meyer-Hermann ME (2007) Are T cells at the origin of B cell lymphomas? J Theor Biol 244: 656–669. [PubMed]

27. Evens AM, Schumacker PT, Helenowski IB, Singh AT, Dokic D, et al. (2008) Hypoxia inducible factor-alpha activation in lymphoma and relationship to the thioredoxin family. Br J Haematol 141: 676–680. [PMC free article] [PubMed]

28. Frieboes HB, Edgerton ME, Fruehauf JP, Rose FR, Worrall LK, et al. (2009) Prediction of drug response in breast cancer using integrative experimental/computational modeling. Cancer Res 69: 4484–4492. [PMC free article] [PubMed]

29. Sinek JP, Sanga S, Zheng X, Frieboes HB, Ferrari M, et al. (2009) Predicting drug pharmacokinetics and effect in vascularized tumors using computer simulation. J Math Biol 58: 485–510. [PMC free article] [PubMed]

30. Frieboes HB, Lowengrub JS, Wise S, Zheng X, Macklin P, et al. (2007) Computer simulation of glioma growth and morphology. NeuroImage 37Suppl 1: S59–70. [PMC free article] [PubMed]

31. Wise SM, Lowengrub J, Frieboes HB, Cristini V (2008) Three-dimensional multispecies nonlinear tumor growth–I Model and numerical method. J Theor Biol 253: 524–543. [PMC free article] [PubMed]

32. Frieboes HB, Jin F, Chuang YL, Wise SM, Lowengrub JS, et al. (2010) Three-dimensional multispecies tumor growth-II: Tumor invasion and angiogenesis. J Theor Biol 264: 1254–1278. [PMC free article] [PubMed]

33. Chauviére AH, Hatzikirou H, Lowengrub JS, Frieboes HB, Thompson AM, et al. (2010) Mathematical Oncology: How Are the Mathematical and Physical Sciences Contributing to the War on Breast Cancer? Curr Breast Cancer Rep 2: 121–129. [PMC free article] [PubMed]

34. Deisboeck TS, Wang Z, Macklin P, Cristini V (2011) Multiscale cancer modeling. Annu Rev Biomed Eng 13: 127–155. [PubMed]

35. Hatzikirou H, Chauviere A, Bauer AL, Leier A, Lewis MT, et al. (2012) Integrative physical oncology. Wiley Interdiscip Rev Syst Biol Med 4: 1–14. [PMC free article] [PubMed]

36. Edgerton M, Chuang Y, Macklin P, Yang W, Bearer EL, et al. (2011) A Novel, Patient-Specific Mathematical Pathology Approach For Assessment Of Surgical Volume: Application To Ductal Carcinoma In Situ Of The Breast. Anal Cell Pathol 34: 1–17. [PMC free article] [PubMed]

37. Macklin P, Edgerton ME, Thompson AM, Cristini V (2012) Patient-calibrated agent-based modeling of ductal carcinoma in situ (DCIS): From microscopic measurements to macroscopic predictions of clinical progression. J Theor Biol 301: 122–140. [PMC free article] [PubMed]

38. Morin PJ (2003) Drug resistance and the microenvironment: nature and nurture. Drug Resist Updates 6: 169–172. [PubMed]

39. Schmitt CA, Fridman JS, Yang M, Lee S, Baranov E, et al. (2002) A senescence program controlled by p53 and p16INK4a contributes to the outcome of cancer therapy. Cell 109: 335–346. [PubMed]

40. Lowe SW, Sherr CJ (2003) Tumor suppression by Ink4a–Arf: progress and puzzles. Curr Opin Genet Dev 13: 77–83. [PubMed]

41. Tomida A, Tsuruo T (1999) Drug resistance mediated by cellular stress response to the microenvironment of solid tumors. Anti-Cancer Drug Design 14: 169–177. [PubMed]

42. Sobocinski GP, Toy K, Bobrowski WF, Shaw S, Anderson AO, et al. (2010) Ultrastructural localization of extracellular matrix proteins of the lymph node cortex: evidence supporting the reticular network as a pathway for lymphocyte migration. BMC Immunology 11: 42. [PMC free article] [PubMed]

43. Ma B, Jablonska J, Lindenmaier W, Dittmar KEJ (2007) Immunohistochemical study of the reticular and vascular network of mouse lymph node using vibratome sections. Acta Histochemica 109: 15–28. [PubMed]

44. Koster A, Raemaekers JM (2005) Angiogenesis in malignant lymphoma. Curr Opin Oncol 17: 611–616. [PubMed]

45. Ruan J, Hajjar K, Rafii S, Leonard JP (2009) Angiogenesis and antiangiogenic therapy in non-Hodgkin's lymphoma. Annals of Oncology 20: 413–424. [PMC free article] [PubMed]

46. Anderson AR, Chaplain MA (1998) Continuous and discrete mathematical models of tumor-induced angiogenesis. Bull Math Biol 60: 857–899. [PubMed]

47. Plank MJ, Sleeman BD (2003) A reinforced random walk model of tumour angiogenesis and anti-angiogenic strategies. Math Med Biol 20: 135–181. [PubMed]

48. Plank MJ, Sleeman BD (2004) Lattice and non-lattice models of tumour angiogenesis. Bull Math Biol 66: 1785–1819. [PubMed]

49. Nugent LJ, Jain RK (1984) Extravascular diffusion in normal and neoplastic tissues. Cancer Res 44: 238–244. [PubMed]

50. Sherratt JA, Murray JD (1990) Models of epidermal wound healing. Proc Roy Soc Lond B241: 29–36. [PubMed]

51. Cristini V, Lowengrub J, Nie Q (2003) Nonlinear simulation of tumor growth. J Math Biol 46: 191–224. [PubMed]

52. Li X, Lowengub J, Rätz A, Voigt A (2009) Solving PDEs in complex geometries: a diffuse domain approach. Commun Math Sci 7: 81–107. [PMC free article] [PubMed]

53. Ito K, Smith BR, Parashurama N, Yoon J-K, Song SY, et al. (2012) Unexpected dissemination patterns in lymphoma progression revealed by serial imaging within a murine lymph node. Cancer Res 72: 6111–6118. [PubMed]

54. Frieboes HB, Zheng X, Sun CH, Tromberg B, Gatenby R, et al. (2006) An Integrated Computational/Experimental Model of Tumor Invasion. Cancer Res 66: 1597–1604. [PubMed]

55. Li X, Cristini V, Nie Q, Lowengrub J (2007) Nonlinear three-dimensional simulation of solid tumor growth. Discrete Cont Dyn Sys B 7: 581–604.

56. Cristini V, Lowengrub JS (2010) Multiscale Modeling of Cancer: An Integrated Experimental and Mathematical Modeling Approach. Cambridge University Press .

57. Grantab R, Sivananthan S, Tannock IF (2006) The penetration of anticancer drugs through tumor tissue as a function of cellular adhesion and packing density of tumor cells. Cancer Res 66: 1033–1039. [PubMed]

58. Drillenburg P, Pals ST (2000) Cell adhesion receptors in lymphoma dissemination. Blood 95: 1900–1910. [PubMed]

59. Muniz VP, Barnes JM, Paliwal S, Zhang X, Tan X, et al. (2011) The ARF Tumor Suppressor Inhibits Tumor Cell Colonization Independent of p53 in a Novel Mouse Model of Pancreatic Ductal Adenocarcinoma Metastasis. Mol Cancer Res 9: 867–877. [PMC free article] [PubMed]

60. Paliwal S, Pande S, Kovi RC, Sharpless NE, Bardeesy N, et al. (2006) Targeting of C-terminal binding protein (CtBP) by ARF results in p53-independent apoptosis. Mol Cell Biol 26: 2360–2372. [PMC free article] [PubMed]

61. Paliwal S, Kovi RC, Nath B, Chen YW, Lewis BC, et al. (2007) The alternative reading frame tumor suppressor antagonizes hypoxia-induced cancer cell migration via interaction with the COOH-terminal binding protein corepressor. Cancer Res 67: 9322–9329. [PubMed]

62. Chen YW, Paliwal S, Draheim K, Grossman SR, Lewis BC (2008) p19Arf inhibits the invasion of hepatocellular carcinoma cells by binding to C-terminal binding protein. Cancer Res 68: 476–482. [PMC free article] [PubMed]

63. Takemoto S, Mulloy JC, Cereseto A, Migone TS, Patel BK, et al. (1997) Proliferation of adult T cell leukemiaylymphoma cells is associated with the constitutive activation of JAK/STAT proteins. Proc Natl Acad Sci USA 94: 13897–13902. [PubMed]

64. Graeber TG, Osmanian C, Jacks T, Housman DE, Koch CJ, et al. (1996) Hypoxia-mediated selection of cells with diminished apoptotic potential in solid tumours. Nature 379: 88–91. [PubMed]

65. Semenza GL (2000) Hypoxia, clonal selection, and the role of HIF-1 in tumor progression. Crit Rev Biochem Mol Biol 35: 71–103. [PubMed]

66. Freitas AA, Rocha B, Coutinho AA (1986) Life span of B lymphocytes: the experimental basis for conflicting results. J Immunol 136: 470–476. [PubMed]

67. Young AJ, Hay JB, Miyasaka M (1995) Rapid turnover of the recirculating lymphocyte pool in vivo. Int Immunol 7: 1607–1615. [PubMed]

Articles from PLoS Computational Biology are provided here courtesy of **Public Library of Science**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's Canada Institute for Scientific and Technical Information in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |