Department of Microbiology and Molecular Genetics, University of California, Irvine, Irvine, CA, USAChao Family Comprehensive Cancer Center, University of California, Irvine, Irvine, CA, USACenter for Complex Biological Systems, University of California, Irvine, Irvine, CA, USA
Department of Mathematics, University of California, Irvine, Irvine, CA, USAChao Family Comprehensive Cancer Center, University of California, Irvine, Irvine, CA, USACenter for Complex Biological Systems, University of California, Irvine, Irvine, CA, USADepartment of Biomedical Engineering, University of California, Irvine, Irvine, CA, USA
Figure 1.SW480 xenograft tumors reveal a spotted pattern of metabolic heterogeneity
A, B SW480 cells lentivirally transduced with empty pCDH vector (mock) were subcutaneously injected into immunocompromised mice. The resulting tumors were stained for (A) phosphorylated pyruvate dehydrogenase (pPDH) and counterstained with hematoxylin or (B) lymphoid enhancer factor‐1 (LEF‐1). Scale bars indicate 100 μm in the series of 4×, 20×, and 40× images. The red curves denote spot contours and the blue curves denote convex hulls, which group together spots that are sufficiently close to one another (see Appendix A1).
C. Image analysis of spot size versus distance of spot to nearest neighbor, using analyzed 20× images (third panels of A and B). The outlined data points indicate the average distance and area for pPDH and LEF‐1 spots. Results show that quantifiable features of the spotted patterns in pPDH and LEF‐1 are similar.
D. Colorectal carcinoma patient samples (tumors 1, 2, and 3) stained for pPDH (top) and LEF‐1 (bottom) show spatial heterogeneity in expression levels. Scale bars are 100 μm (LEF‐1 samples from Uhlén et al, 2015).
E. Serial section of human colorectal carcinoma stained with pPDH and LEF‐1 antisera. Scale bars are 100 μm.
Figure EV1.Overlap of pPDH and LEF1 spots in xenograft tumors
Convex hull image analysis of serial sections of SW480 mock xenograft tumors stained for pPDH and LEF1, as shown in Fig 1A and B, second panels.
A, B Isolated contour maps with convex hull outlines for pPDH and LEF1.
C–E pPDH and LEF1 contour maps were overlaid on each other and overlapped regions highlighted in blue. Different thresholds for spot detection were tested; each threshold condition revealed between 65 and 77% overlap between pPDH and LEF1 spots.
Figure 2.A mathematical model for Wnt signaling regulation of metabolism
This set of reaction–diffusion equations describes the change over time of oxidative (Po) and glycolytic (Pg) cell populations, Wnt signaling activity (W), and Wnt inhibitor activity (WI).
The cells can diffuse, proliferate, and “switch” metabolism programs depending on Wnt signaling activity and nutrient levels and die from lack of nutrient (N).
Wnt and Wnt inhibitor activity equations are based on the Gierer–Meinhardt activator–inhibitor model. The Wnt signal diffuses short range relative to the longer‐range diffusion of the Wnt inhibitor. Wnt also auto‐upregulates its activity in glycolytic cells at a rate proportional to nutrient level, is inhibited by a Wnt inhibitor, is constitutively upregulated in both cell types, and decays (downregulation term). The Wnt inhibitor diffuses long range, is nonlinearly upregulated by Wnt, and decays. Equations for nutrient and dead cells (Pd) are not shown; their descriptions are in the main text.
Three‐dimensional numerical simulations that model the spatial distribution and level of glycolytic and oxidative cells, Wnt, and Wnt inhibitor reveal an emergent self‐organizing pattern of metabolic heterogeneity (spots). The simulations shown depict the heterogeneity in a 3D and 2D representation. The 3D representation includes a portion of the “tumor” removed to visualize the interior of the domain. The 2D representation is a horizontal slice of the respective 3D simulation in the center of the domain. Color bars refer to unitless concentrations.
Summary of parameter effects on the spotted pattern.
Figure 3.Decreasing Wnt signaling leads to changes in metabolic patterning in xenograft tumors
A, B SW480 cells were lentivirally transduced to express dominant negative LEF‐1 (dnLEF‐1), and transduced cells were injected subcutaneously into immunocompromised mice. Tumor sections were stained for phosphorylated PDH (A) and β‐catenin (B) and counterstained with hematoxylin. Compared to mock tumors, the spots are larger and more heterogeneous and the background staining is lighter, which reflects an overall decrease in Wnt signaling. Scale bars are 100 μm. The red curves denote spot contours and the blue curves denote convex hulls, which group spots that are sufficiently close to one another (see Appendix A1).
C. Numerical simulations that lower Wnt signaling activity in the model show an overall decrease in glycolysis and a change in the spotted pattern that closely mimics that observed in the dnLEF tumors. Color bars refer to unitless concentrations.
D. Image analysis of spot size versus distance of spot to nearest neighbor, using analyzed images. Averages for mock and dnLEF spot simulations are denoted in white outlined symbols (pPDH: spot sizes/inter‐spot distances: mock simulation average: 225 ± 278 μm2/29 ± 12 μm; mock xenograft tumor average: 309 ± 367 μm2/29 ± 10 μm; dnLEF‐1 simulation average: 423 ± 327 μm2/41 ± 14 μm; dnLEF‐1 xenograft tumor average: 1,139 ± 1,042 μm2/53 ± 15 μm). Results are also shown for dominant negative transcription factor 1 (dnTCF‐1) tumors (see Appendix A1.10–A1.12 and Appendix Figs S8–S10). The metabolic pattern in dnTCF‐1 tumors is consistent with that in dnLEF‐1 tumors. The analysis and model predict that the changes in the metabolic spotted pattern (larger spots, greater distance between spots) are due to an increase in the diffusion of Wnt and the Wnt inhibitor.
E. Comparison of mock β‐catenin spots to dnLEF‐1 and dnTCF‐1 β‐catenin spots from image analysis. Averages for mock and dnLEF‐1 spot simulations are denoted in white outlined symbols (β‐catenin spot sizes/inter‐spot distances: mock simulation average: 139 ± 145 μm2/26 ± 11 μm; mock xenograft tumor average: 97 ± 209 μm2/16 ± 4 μm; dnLEF‐1 simulation average: 342 ± 221 μm2/39 ± 14 μm; dnLEF‐1 xenograft tumor average: 312 ± 277 μm2/32 ± 9 μm; dnTCF‐1 xenograft tumor average: 603 ± 578 μm2/42 ± 14 μm). The analysis and model predict that Wnt signaling diffuses further with dominant negative LEF‐1 expression.
Figure 4.Model predictions revealed in xenograft tumors and human colorectal cancer
The model predicted that lowering Wnt signaling results in an increase in the expression of factors that increase the range of diffusion of Wnt and Wnt inhibitors.
A. Known regulators of Wnt ligand diffusion.
B, C Quantitative PCR of diffusion regulators in SW480 mock, dnLEF‐1, and dnTCF‐1 (B) transduced cells, and (C) xenograft tumors show human SPOCK2, GPC4, and SFRP5 mRNA are notably upregulated in xenograft tumors but not 2D in vitro culture. In vitro data represent an average of three sample sets (± SEM), and xenograft tumor data represent the average of five independent tumor sets (± SEM); * denotes P < 0.05. Statistical significance was determined using Student's two‐tailed t‐test.
D, E Gene expression data, from GEO dataset GDS3756, of 21 rectal cancer patient tissue with or without radio‐chemotherapy (Snipstad et al, 2010). Significant changes in expression levels of GPC1 (P = 0.00019), SFRP1 (P = 0.016), SFRP2 (P = 0.0006), and SFRP4 (P = 0.0006) are observed in post‐therapy tumor cells compared to before treatment. SD shown. Statistical significance was determined using the Mann–Whitney U‐test with Benjamini–Hochberg correction for multiple hypotheses (RStudio).
Figure EV2.Wnt ligand and glycolysis gene expression in rectal cancer patients post‐radiochemotherapy
Gene expression data, from GEO dataset GDS3756, of 21 rectal cancer patients before and after radio‐chemotherapy (Snipstad et al, 2010).
A, B Expression of Wnt ligands WNT5B, WNT8B, and WNT10B shows trends toward increased expression in rectal tumor tissue treated with radiochemotherapy, but these changes do not reach statistical significance when P <0.05 is used as a cutoff (specific P‐values are 0.10, 0.10, and 0.08, respectively). Statistical significance was determined using the Mann–Whitney U‐test with Benjamini–Hochberg correction for multiple hypothesis testing.
C, D Expression of the glycolytic enzyme ENO2 is specifically increased in tumor tissue after radiochemotherapy (P =0.008); HIF1A expression also shows a trend in increased expression (P =0.06). * denotes adjusted P‐value < 0.05; + denotes adjusted P‐value < 0.10. Statistical significance was determined using the Mann–Whitney U‐test with Benjamini–Hochberg correction for multiple hypothesis testing.
Figure 5.Simulations identify the glycolytic cell population as a sensitive drug target
We target either cells with more oxidative phosphorylation (Po; left) or cells with more glycolysis (Pg; right) selectively, starting from a metabolically patterned state, for 2.5, 5, or 7.5 (arbitrary) time units, with a treatment dose between 0.25 and 1. After therapy is halted, the cells are allowed to evolve according to the original model (Fig 2). The total cell populations, relative to the initial, starting cell population are shown. Corresponding populations of Po and Pg cells can be found in Fig EV3.
Figure EV3.Targeted therapy simulations for Pg and Po populations
Figure 5 gives the results for total tumor size after individually targeting either Pg or Po cells with given treatment doses and treatment times. Here are the effects for the individual Pg and Po populations in those targeted therapy simulations. Simulations suggest that the glycolytic cell population is a more sensitive drug target than the oxidative cell population. We target either Po (left) or Pg (right) cells selectively, starting from a metabolically patterned state, for 2.5, 5, or 7.5 (arbitrary) time units, with a death rate between 0.25 and 1. After therapy is stopped, the cells are allowed to evolve according to the original model (Fig 2). The Pg and Po cell populations, relative to their initial cell populations, are shown.
Figure 6.Therapies targeting metabolism and Wnt synergize for tumor death in mathematical simulations
A, B Modeled therapies, their targets, and the model parameters influenced by therapy.
C. Starting with a metabolically patterned state, treatment of tumors with dichloroacetic acid (DCA) and XAV939 combined leads to an effective crash in the system, as shown by the complete loss of cells (1e and 2e). The panels on the left show the cell arrangements for the oxidative (Po) and glycolytic (Pg) populations (metabolic patterning), and the three graphs on the right show the fractions of Po or Pg cells relative to their initial values, after applying the therapy for 50 (arbitrary) time units. The effects of therapy on the total cell population, relative to the initial cell population, for the same DW and 1/τgo values, are shown in the third graph. XAV939 treatment is modeled by decreasing SW and increasing DW and DWI linearly with respect to the decrease in SW; legend values are listed relative to mock SW. The dashed curves labeled 0.80 correspond to the case in which 0.80, but DW and DWI are unchanged and take the values used in the mock tumor simulations. The panels on the left correspond to the red curve in the graphs and show the effect on patterning for 1/τgo values 1, 4, 12, 17, and 18, respectively [denoted by labels (1a) through (2e)]. Color bar refers to unitless concentrations.
Combined Wnt signaling and glycolysis targeting therapies significantly decrease SW480 spheroid size in vitro.
SW480 cells were embedded in a fibrin gel using the method shown. Media containing a mock treatment, 0.5 mM DCA, 2 mM DCA, 10 μM of XAV939, or a combination of DCA and XAV939.
Representative 4× images of spheroids each condition, imaged 14 days after treatment, are shown.
Analysis of 75 spheroids per condition shows 2 mM DCA significantly increases SW480 spheroid size, while combined 2 mM DCA treatment with 10 μM XAV939 significantly decreases their size. Statistical significance was determined using Student's two‐tailed t‐test.
The effects of therapy on the total cell population, relative to the initial cell population, of combined XAV939 and DCA treatment were simulated using an in vitro version of the model (Appendix A7). As in Fig 6, DCA treatment was modeled by increasing the rate at which cells switched from a glycolytic metabolic phenotype to an OXPHOS phenotype (e.g., 1/τgo is increased) to reflect the tendency of cells to perform OXPHOS when PDK is inhibited. XAV939 treatment was modeled by decreasing the general Wnt signaling term SW and increasing the range of Wnt and its inhibitor (due to upregulation of Wnt and Wnt inhibitor diffusers), which we modeled by increasing DW and DWI proportionally. The dashed curves labeled 0.80 correspond to the case in which 0.80 but DW and DWI are unchanged.