Identification of Genes Regulating Breast Cancer Dormancy in 3D Bone Endosteal Niche Cultures Free

The progression of breast cancer at metastatic sites continues to be the largest contributor to patient mortality (1, 2). There has been renewed focus on identifying the mechanisms governing the establishment of dormancy in specific peripheral sites and the reawakening of dormant cells as major contributors of this cancer's lethal phenotype (3), especially because only a fraction of disseminated tumor cells, and even fewer circulating tumor cells, give rise to clinical macrometastases (4–6). Both estrogen receptor (ER)-positive and triple-negative breast cancers (TNBC: ER⁻/PR⁻/Her2⁻) metastasize to and enter dormancy in the bone, with ER⁺ tumors often exhibiting a longer-lasting dormancy (7). Importantly, dormant versus active growth depends on which bone niche is colonized (8). Colonization of the endosteal niche (EN), enriched in osteoblasts and marked by low oxygen and high calcium levels, results in dormancy, whereas colonization of the perivascular niche, enriched in hematopoietic stem cells, results in active proliferation and formation of macrometastases (7, 9). Signaling pathways in tumor cells define whether disseminated tumor cells will remain dormant or become proliferative: activation of the p38 MAPK pathway in the absence of ERK1/2 MAPK activity favors quiescence (10). In addition, long-term survival of dormant cells is likely to require adoption of stem-like quiescence properties through increased activity of lineage plasticity pathways. Although several groups identified genes and pathways differentially regulated in dormant versus proliferating breast cancer metastatic cells (8, 11–16), few studies have addressed possible causative roles for these genes, especially in the context of p38 control (16, 17), confounding attempts at therapeutic targeting.

Using a novel 3D-EN culture system developed by Marlow and colleagues (18), in which otherwise aggressive MDA-MB-231 human breast cancer cells become dormant for up to 60 days due to direct contact with EN environmental cells, we show here that the knockdown of p38-MAPK induced proliferation, confirming the notion that dormancy is p38-dependent. To identify genes that promote or maintain dormancy, we used a high-stringency genomic short hairpin RNA (shRNA) screen to identify gene knockdowns that induced proliferation of MDA-MB-231 in 3D-EN cultures. Identification of the top hits, BHLHE41, a known p38 target required for induction of quiescence (19), and LSP1, a suppressor of ERK1/2 activity (20) increased confidence in the screen's validity. Other gene hits grouped into several regulatory categories not previously identified as suppressors of dormancy. These include genes involved in neurogenesis, translation, and noncanonical WNT signaling, all of which play roles in regulating the maintenance of cancer stem cells. These data strengthen the notion that breast cancer dormancy in the EN is promoted by p38-MAPK–controlled stem cell pathways.

MDA-MB-231 (ATCC HTB-26), MDA-MB-231 BoM1833 (gift from Joan Massague, Memorial Sloan Kettering Cancer Center; ref. 21), HEK293T (ATCC CRL-11268), Phoenix 293T (ATCC CRL-3213), and HS-5 (ATCC CRL-11882) were cultured in DMEM supplemented with 10% FBS and 1% penicillin–streptomycin and incubated at 37°C and 5% CO₂. Human fetal osteoblasts (hFOB; ATCC CRL-11372) were cultured with DMEM/F12 (1:1) media without phenol red and supplemented with 15% FBS and G418 (0.3 mg/mL) and incubated at 32°C and 5% CO₂. Human umbilical vein endothelial cell (HUVEC; Lonza C2517A) were cultured to no more than five passages with EBM-2 media (Lonza CC-3156) supplemented with a SingleQuot pack (Lonza CC-4176) and grown at 37°C and 5% CO₂.

Three-dimensional (3D) cultures recapitulating growth in the EN or control 3D collagen matrix growth were produced as described previously (18). Briefly, GELFOAM (Pfizer) discs were cut (5-mm diameter and 3-mm thickness) using a sterile hole punch and scalpel and then UV treated for 2 hours (1.2 × 10⁵ μ J/CM²) using a UV cross linker. GELFOAM discs were placed into a low attachment 96-well plate (Eppendorf Microplate 96/F-PP, Sigma-Aldrich #C150179G) using sterile forceps. Discs were incubated in 200 μL of 1× Dulbecco PBS at 37°C/5% CO₂ for 20 minutes. EN cultures were produced by seeding HS-5, hFOB, and HUVECs onto GELFOAM discs (4 × 10⁴ cells/5 μL each) via capillary action, incubating for 2 hours, and then topped off with 200-μL niche media. After 24 hours, media were removed and MDA-MB-231 cells (transduced with RFP-expressing shRNA libraries) were seeded onto EN cultures (10⁴ cells/10 μL) and grown for 7 days, with niche media (100 μL) replaced daily. MDA-MB-231 cells grown on tissue culture plates (2D) or in GELFOAM alone (3D) served as negative controls. Bone-tropic MDA-MB-231 (BoM1833) cells, infected with pGIPZ (GFP-expressing) lentivirus vector, served as a positive control for 3D-EN growth. MDA-MB-231 growth was monitored every other day via fluorescence microscopy using a Nikon Eclipse TS100 inverted microscope and SPOT Insight Fire Wire Camera and SPOT 5.2 software. Cell numbers were quantified using ImageJ software (NIH) from five different fields containing >100 cells/field.

Lentivirus shRNA clones (pGIPZ-based; Supplementary Table S2), Cellecta Human DECIPHER RFP screening libraries (Table 1), and DNA expression vectors (ORFeome 8.1) were provided by the Roswell Park Comprehensive Cancer Center Gene Modulation Core (I.H. Gelman, Director). Plasmids were propagated in Stbl3 bacteria (Thermo Fisher Scientific) in LB media supplemented with 100 μg/mL ampicillin for 18 hours at 37°C at a constant speed of 200 rpm. Plasmids were extracted using QIAprep Spin Miniprep Kit (Qiagen catalog no. 27104) according to the manufacturer's protocol. Relative plasmid concentrations were quantified using a Nanodrop 2000 (Thermo Fisher Scientific).

HEK293T cells were transfected with LipoD293/DNA mixtures as described previously (22).

Polytropic lentiviruses were packaged in HEK293T cells using psPAX2 and pMD2.G packaging constructs as described previously (23). For infection of target cells plated the previous day at confluency in 6-well dishes, 150 μL of lentiviral supernatant and 4 μg/mL polybrene (Sigma) was added to the cells along with 1–2 mL DMEM (10% FBS, 1% penicillin/streptomycin) until cells were completely covered. Cells were then incubated for 30 minutes at 37°C with 5% CO₂, centrifuged at 1,800 rpm for 45 minutes, and then incubated overnight again at 37°C/5% CO₂. The cultures were split at a 1:5 ratio in DMEM containing 10% FBS, 1% penicillin/streptomycin, 2 mg/mL puromycin.

RNA was extracted using TRIzol (Life Technologies) according to the manufacturer's protocol. One microgram of total RNA/reaction was used in qRT-PCR reactions (50 μL total volume) containing High Capacity cDNA Reverse Transcriptase Kit reagents (Life Technologies/Applied Biosystems #4368814). qPCR was performed using Power SYBR green PCR Mastermix (Life Technologies #4367659) on a Step One Plus thermocycler (Applied Biosystems). GAPDH housekeeping gene was used as the loading control. qRT-PCR reactions were performed in triplicate and relative concentrations were determined using the 2^−ΔΔC_t method (StepOne software).

Immunoblot analysis was performed, as described previously (24), using the following antibodies: primary—V5 tag (Thermo Fisher Scientific, #37-7500), GAPDH (Santa Cruz Biotechnology, #sc-25778), α-tubulin (Santa Cruz Biotechnology, #sc-5286), p38-MAPK (Cell Signaling Technology, #9212), p38^poT180/Y182 (Cell Signaling Technology, #9211), ERK1/2 (Cell Signaling Technology, #9102S), ERK1/2^poT202/Y204 (Cell Signaling Technology, #9101S); secondary—Alexa Fluor700 anti-mouse Ig (1:1,000) or Alexa Fluor800 anti-rabbit Ig (1:10,000).

DNA was isolated from 7-day 2D and 3D-EN cultures using phenol/chloroform/isoamyl (25:24:1). Cells in the 3D culture were isolated by removing GELFOAM discs with cells and incubating on a rotator for 1 hour at 37°C in 5 mL 1× collagenase/hyaluronidase solution (STEMCELL Technologies). After vigorous mixing, cells were pelleted and washed three times with PBS. First-round and nested PCR were performed according to the Cellecta shRNA library manual (http://www.cellecta.com/resources/product-manuals-and-certificates/), as we described previously (23) using primers described in Supplementary Table S3. All experiments were conducted in triplicate. The PCR products were cleaned using QIAquick Gel Extraction Kits, and then subjected to single-end rapid mode sequencing on an Illumina HiSeq2500 as performed by the Roswell Park Comprehensive Cancer Center Genomics Shared Resource (Sean Glenn, Director). Using FASTQ sequencing data files, barcodes were trimmed from flanking sequence using the ShortRead package from Bioconductor (25). The isolated barcode sequences were aligned to a reference file matching shRNA clones to gene targets using the DECIPHER BarCode Deconvoluter program (Cellecta), which allows for up to two incorrect base changes for accurate barcode identification. Individual sequence read counts were normalized by total reads sequenced, and top hits were filtered on the basis of a threshold determined by luciferase shRNA-negative controls (16 clones). An analysis of row sums was performed to identify genes targeted by multiple shRNA clones and across replicates.

Statistical analysis was performed on the fold-change between the cell counts from day 1 to 7 using the Student two-tailed t test. Error bars indicate SEM. Significant differences between experimental groups had a P value lower than 0.05.

Using a novel 3D model of dormancy for bone-metastatic breast cancer (18), we endeavored to identify genes that suppress tumor cell quiescence in a cultured microenvironment recapitulating bone EN. In this model, the human TNBC cell line MDA-MB-231 proliferates in a GELFOAM biomatrix, whereas it is growth-arrested in EN conditions (hFOBSs, HUVECs, and HS-5 diploid fibroblasts in GELFOAM; Fig. 1A). Importantly, the inclusion of bone marrow origin fibroblasts (HS-5) and human endothelial cells (HUVEC) promoted the long-term survival of hFOB osteoblasts even after these cells reached initial confluence after 24 hours of growth. This EN culture condition was previously shown to induce growth arrest of ER-positive (MCF7, T47D, ZR75-1, and BT474) and ER-negative (SUM149, SUM159, MDA-MB-231, and MDA-MB-453) human breast cancer cell lines, whereas these lines could proliferate in either GELFOAM alone, or in GELFOAM seeded with primary human bone marrow stem cells, representing a perivascular niche (18). In contrast, the bone-metastatic MDA-MB-231 variant, BoM1833, which was selected in vivo for increased bone growth (26), proliferates in either niche (Fig. 1B). Consistent with the notion that activated p38 MAPK in the absence of MEK-ERK activation favors dormancy, we showed that the knockdown of p38 by shRNA (shRNA clones #15 and #18) also induced MDA-MB-231 proliferation in the EN (Fig. 1C), consistent with previous data (18) using the p38 kinase inhibitor, SB203580.

To identify the suppressors of tumor cell proliferation in a bone niche, MDA-MB-231 cells were transduced with a genomic shRNA library (Cellecta DECIPHER library covering 15,377 human genes with 82,500 independent shRNA clones, divided into 3 modules; Table 1) and clones that proliferated in EN cultures were enriched. Genes that are potentially required for MDA-MB-231 dormancy within the EN were identified by performing next-generation sequencing of shRNA clone barcodes from DNA taken from triplicate screen aliquots of freshly infected cells (24 hours) and from infected cells incubated for 7 days in 3D-EN. The barcode sequences were trimmed from flanking sequences and shRNA-targeted genes then identified using Cellecta's BarCode Deconvoluter software. We selected gene targets (shRNA barcodes) that were found in ≥2 of three independent screens, identified by ≥3 independent shRNA clones/gene, each at >1.5-fold increase over background (normalized against the relative abundance of each clone in the library; Table 1). This analysis identified 416 potential “dormancy-reactivation suppressor” genes (DRSG) in the 3 shRNA clone modules (Table 1).

One of the ways we established statistical significance for potential DSRG candidates was to compare the relative frequency of shRNA clones with the 16 luciferase shRNA controls (shLuc) contained within each module. For example, 2 module-2 genes, DOLK and MICALL2, had at least 2 independent shRNA clones with fold-change sequence reads over the 1.55-fold shLuc threshold (Supplementary Fig. S1A), yet when compared between screening replicates, only the MICALL2 clones showed consistent statistical significance over the luciferase clones (Supplementary Fig. S1B) in ≥2 replicates. Indeed, the knockdown of DOLK using independent shRNAs failed to induce increase MDA-MB-231 proliferation in 3D-EN cultures (Table 2).

Many studies have shown that metastatic dormancy is controlled by the simultaneous upregulation of p38MAPK and downregulation of ERK activation (27), yet little is known about mediators of dormancy downstream of p38MAPK. Thus, DSRG were subjected to Ingenuity Pathway Analysis and PubMed search to bin them based on a predicted or known relationship to p38MAPK signaling and metastasis. An example is shown in Table 3 for the 65 potential DSRGs from module-1, relative to their relationship to p38MAPK signaling and metastasis. Of these candidates, 14 genes were associated with both p38MAPK signaling and metastasis (ADRB2, BHLHE41, CASR, CD63, CDC2L1, FLT1, HBP1, KEAP1, LSP1, NOB1, NRG1, P11, TTF1, WNT3), 12 genes were associated with metastasis but not with p38MAPK signaling (BIK, BRD4, CLDN2, EIF4A1, FUBP1, HSPD1, KIF11, NES, NOX1, RELN, SERPIN1, WNT8A), and 5 genes were associated with p38MAPK signaling alone (CACNB3, EHF, NNAT, OASL, TAC4). As a whole, there was a selection for DSRG candidates in module-1 that are involved in the regulation of neurogenesis or stem cell biology (Table 3), and/or protein translation (Tables 3 and 4). Loss of differentiation (neurogenesis) or stem cell induction genes resulting in active breast cancer cell growth in the 3D-EN is consistent with the notion that disseminated tumor cells exhibit stem cell–like properties (28, 29). Differential expression of genes controlling ribosome biogenesis are known to control stem cell homeostasis (30), and indeed, the antagonism of this process was shown to inhibit tumor formation induced by CD44⁺/CD24⁻ human breast cancer stem cells (31).

In contrast, several gene candidates had roles that might directly control known dormancy functions. For example, BHLHE41 has been reported to play a role in p38MAPK-mediated dormancy (19), ADRB2 suppresses prostate cancer proliferation in bones by downregulating osteoblast-expressed GAS6 (32), LSP1 negatively controls ERK signaling (20), and P11 (S110A10) controls breast cancer adhesion to endothelial cells in the metastatic niche (33). Finally, gene set enrichment analysis of all the module-1 genes showed that 18 of 65 genes (27.7%) were likely MYC targets (Fig. 2A). Although a role of MYC as a driver of dormancy reawakening has not been addressed, MYC amplification is associated with high-grade breast cancer and worse prognosis (34), and in several non-breast cancer models, the forced reexpression of MYC rescues proliferation in dormant tumor cells (35).

We then sought to independently confirm that the downregulation of several DSRG candidates leads to MDA-MB-231 proliferation in 3D-EN cultures. Thus, MDA-MB-231 cells were transduced with 2 independent shRNA clones/gene, and following confirmation of gene knockdown by either qRT-PCR or immunoblotting, the cells were assessed for proliferation (vs. scrambled shRNA controls) in 3D-EN cultures as in Fig. 1. For this analysis, we chose three predicted DSRG from module-1, BHLHE41, HBP1, and WNT3, which were both p38- and metastasis-associated (Table 3), and 1 gene, NES, not known to be p38-regulated. As well, we chose two negative controls (not predicted to be DSRG): DOLK, which was neither p38- nor metastasis-associated (Supplementary Table S1) and that was likely not significant due to lack of replicate hits (Supplementary Fig. S1B), and HTATIP2, a module-2 gene that failed to make the cut because it had only 2 shRNA hits in one of three replicates. BHLHE41, HBP1, and WNT3 were validated as DSRG, that is, their knockdown resulted in significantly increased proliferation in EN over controls (BHLHE41: Table 2, HBP1 and WNT3: Fig. 2C and D), whereas DOLK, NES, and HTATIP2 knockdown failed to induce MDA-MB-231 proliferation in the 3D-EN (Table 3).

WNT signaling largely has been linked to metastatic progression, especially in models of TNBC (36). However, recent data suggest that specific WNT family members, such as WNT5A, might promote either metastatic progression or dormancy, depending on whether signaling is through canonical or noncanonical pathways (37). Therefore, we sought to follow-up our finding that WNT3 knockdown induced MDA-MB-231 proliferation in our 3D-EN assay. First, we confirmed that WNT deficiency caused proliferation (Fig. 2D), using two independent WNT3-specific shRNAs, which knocked down the WNT3 expression in MDA-MB-231 cells approximately 2.5-fold over scrambled controls (shCont; Fig. 2E). Interestingly, WNT3 levels were relatively decreased in BoM1833 cells, consistent with the notion that WNT3 loss facilitates proliferation in the 3D-EN (Fig. 2F). However, WNT3 knockdown in MDA-MB-231 cells had no effect on relative p38^poT180/Y182 levels, indicating that the increased proliferation of MDA-MB-231 cells after WNT3 deficiency was not a result of loss of p38 activation. Although this finding would be consistent with WNT3 being a downstream mediator of p38 signaling, we cannot rule out that the WNT3 effect is p38-independent. In contrast, the BoM1833 variant, which failed to growth-arrest in the 3D-EN cultures (Fig. 1B), exhibited decreased relative p38 activation levels.

Secreted forms of frizzled-related proteins (SFRP) are thought to antagonize WNT/β-catenin signaling by directly binding WNT members (38), and more specifically, SFRP2 and 4, but not SFRP3, are known to bind WNT3 at high affinity (39). We therefore transduced MDA-MB-231 cells with V5-epitope tagged SFRP2, 3, or 4 expressing vectors, confirmed ectopic expression versus an empty vector control (Fig. 2H), and tested these cells for proliferation in 3D-EN cultures. Figure 2G shows that SFRP2 and 4, but not 3, could release MDA-MB-231 cells from dormancy. Similarly, the overexpression of WNT3 in BoM1833 suppressed proliferation in 3D-EN, but not in 3D-control cultures (Fig. 2G). Taken together, these data strongly suggest that WNT3 promotes dormancy in our 3D-EN model. This is consistent with a report showing that increased Hedgehog-mediated Sfrp1 expression in liver stroma increased the metastatic potential of human Capan-1 pancreatic tumor cells through the suppression of WNT3 signaling, and that the overexpression of WNT3A in Capan-1 cells decreased experimental metastasis formation (40).

To address this in a clinical context, we compared the relative WNT3 expression levels in primary-site versus metastatic breast cancers in The Cancer Genome Atlas (TCGA) and Radvanyi Oncomine datasets (41). These data show lower levels of WNT3 in clinical macrometastases compared with primary-site breast cancer (Fig. 3A), suggesting that WNT3 deficiency promotes active metastatic progression. One drawback, however, is that both studies have very few metastatic cases (three each), with none derived from bone, confounding the determination of statistical significance. We then analyzed how the loss of WNT3 correlates with either progression-free or overall survival using the TCGA Breast and METABRIC (42, 43) datasets in cBioPortal (http://www.cbioportal.org). We noted that WNT3 loss, either due to gene deletion or transcriptional downregulation, occurred in 31% and 21% of all breast cancer cases in TCGA and METABRIC, respectively, and that these cases showed strong cooccurrence with the loss/mutation of TP53 (P values of 1.36e-21 and 2.98e-58, respectively). Indeed, clinical cases of TNBC are marked by specific TP53 mutations (44). We then sought to determine whether the loss of WNT3 and/or TP53 predicted poorer clinical survival, or whether any survival correlation associated with WNT3 loss was potentiated in a background of ER, progesterone receptor (PR) and HER2 loss (encoded by ESR1, PGR, and ERBB2, respectively), reflecting TNBC. Indeed, a large portion of the WNT3-deficient cases in the TCGA breast dataset had coincident losses of ESR1, PGR, and ERBB2 (Fig. 3B). Figure 3C shows that in the TCGA dataset, poorer progression-free survival was only detected in cases with combined losses of WNT3, TP53, ESR1, PGR, and ERBB2; the loss of any of these genes alone or the combined loss of ESR1, PGR, and ERBB2 did not affect survival (WNT3 loss alone: P = 0.652; TP53 loss alone: P = 0.125; WNT3 plus TP53 loss: P = 0.0855; ESR1, PGR, and ERBB2 combined loss: P = 0.300). It is noteworthy that BRCA1 mutational status, which represents fewer than 3% of all the WNT3-deficient cases in the TGCA Breast database, has no effect on survival (WNT3 loss + BRCA1 mutation: P = 0.693). In contrast, loss of only WNT3 and TP53 in the METABRIC dataset showed poorer survival (Fig. 3D). Taken together, these data show that WNT3 loss contributes to poorer survival, especially in the context of TP53 loss. The superior powering of the METABRIC dataset, and the fact that it includes many more cases of disease recurrence/progression associated with metastasis, allows for the conclusion that WNT3/TP53 loss is sufficient for poorer survival, a value that worsens in the context of what are likely triple-negative cases (ESR1, PGR, and ERBB2 loss). These data correlate with the fact that TNBC dormancy in the bones is shorter in duration than that of ER-positive breast cancer (7), suggesting that the combined loss of ER, PR, and HER2 might sensitize dormant breast cancer cells toward WNT3 loss. It is important to note that we only studied the role of WNT3 signaling in tumor cells in the context of a 3D-EN microenvironment. Although secreted factors are likely to play an important role, it should be noted that cell–cell contact was required for MDA-MB-231 dormancy in this 3D-EN (18). The exact roles played by each EN niche cell type, whether in its direct interaction to breast cancer cells or through its secretome, remains to be elucidated. In addition, the role of WNT3 in controlling tumor dormancy in the bone may be cancer type–dependent because Nandana and colleagues (45) showed that prostate cancer invasiveness and bone colonization required TBX2-regulated WNT3A expression.

We performed a similar analysis on the two other validated DSRG, which indicated that BHLHE41, but not HBP1, exhibited decreased expression in breast cancer metastases compared to levels in primary tumors (Fig. 4A; Supplementary Fig. S1C). In the TCGA Breast dataset, 11% of cases exhibited partial or full loss of BHLHE41 (Fig. 4B). As with WNT3 loss, the loss of BHLHE41 had a statistical cooccurrence with mutation/loss of TP53 (q = 0.0242) in this dataset. Moreover, the combined loss of BHLHE41 and TP53 correlated with decreased progression-free or overall survival, respectively, in the TCGA Breast and METABRIC datasets (TCGA: P = 0.0379; METABRIC: P = 6.669e-5). In contrast, loss of BHLHE41 alone did not correlate with decreased survival in TCGA Breast (P = 0.568). Interestingly, the combined loss of WNT3, BHLHE41, and TP53 did significantly change the rate of progression-free survival in TCGA Breast cases, strongly suggesting that the WNT3/TP53 and BHLHE41/TP53 loss cohorts were independent groups, and that either loss of WNT3 or BHLHE41 was individually capable of initiating reawakening in the context of TP53 loss. This also suggests that TP53 loss is the main driver of the poorer prognosis. In regards to possible mechanisms underlying BHLHE41 as a DSRG, the gene product, BHLHE41 (also known as DEC2 and SHARP1), functions as a transcriptional repressor of epithelial-to-mesenchymal transition and invasion factors, SNAI1, SNAI2, and TWIST (46). In addition, Adorno and colleagues (47) showed that breast cancer cases with higher levels of BHLHE41 and CCNG2, two p63-induced genes, correlated with lower metastatic risk. Interestingly, specific sets of p53 mutations abrogate p63 activity, likely leading to BHLHE41 loss (48).

In conclusion, this study marks a novel method to identify and validate potential DSRG based on an in vitro 3D-EN dormancy model for breast cancer. Our data suggest several therapeutic avenues, but these would likely be divided into treatments that either secure dormancy, that is, antagonize reawakening, or that induce large-scale reawakening in a neoadjuvant setting, linked to standard chemotherapies prescribed for TNBC. Examples of reawakening suppressors might include small-molecule inhibitors of MYC (49) or Nutlin-3a to normalize mutant p53 function (50), whereas inducers of reawakening might include inhibitors of WNT3 signaling (51, 52) or p38 kinase activity (53).

I.H. Gelman is a consultant/advisory board member for Gerson Lehman Group and Kinex LLC (now Athenex Inc.). No potential conflicts of interest were disclosed by the other authors.

Conception and design: I.H. Gelman

Development of methodology: J.E. McGrath, R. Ransom

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): J.E. McGrath, L. Panzica, H.G. Withers

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J.E. McGrath, L. Panzica, H.G. Withers

Writing, review, and/or revision of the manuscript: L. Panzica, H.G. Withers, I.H. Gelman

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): L. Panzica, H.G. Withers

Study supervision: I.H. Gelman

The authors thank G. Dontu for critical discussion regarding the 3D bone growth models. This work was supported by the Roswell Park Alliance Foundation and by NCI grant P30-CA016056 involving the use of Roswell Park Comprehensive Cancer Center's Genomics, Bioinformatics, and Gene Modulation Shared Resource. This work was also supported by grant CA94108 (NIH/NCI) and by an Alliance Foundation grant (to I.H. Gelman), and in part, through NCI Comprehensive Cancer funds (P30-CA016056) involving the use of Roswell Park Comprehensive Cancer Center's Genomics and Gene Modulation Shared Resources.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.