Evolutionary dynamics at the tumor edge reveal metabolic imaging biomarkers

Juan Jiménez-Sánchez, Jesús J. Bosque, Germán A. Jiménez Londoño, David Molina-García, Álvaro Martínez, Julián Pérez-Beteta, Carmen Ortega-Sabater, Antonio F. Honguero Martínez, Ana M. García Vicente, Gabriel F. Calvo, Víctor M. Pérez-García

Edited by Michael E. Phelps, University of California, Los Angeles School of Medicine, Los Angeles, CA, and approved January 4, 2021 (received for review August 26, 2020)

Author contributions: G.F.C. and V.M.P.-G. designed research; J.J.-S., J.J.B., D.M.-G., Á.M., C.O.-S., and G.F.C. performed research; J.J.-S., J.J.B., G.A.J.L., D.M.-G., Á.M., J.P.-B., A.F.H.M., A.M.G.V., G.F.C., and V.M.P.-G. analyzed data; J.J.-S., J.J.B., G.F.C., and V.M.P.-G. wrote the paper; and G.A.J.L., J.P.-B., A.F.H.M., A.M.G.V., and V.M.P.-G. collected and processed the data and analyzed the medical implications.

¹J.J.-S. and J.J.B. contributed equally to this work.

²G.F.C. and V.M.P.-G. contributed equally to this work.

https://doi.org/10.1073/pnas.2018110118, Volume: 118, Issue: 6, Pages: 1-9

Article Type: Research Article Article History

Publisher: National Academy of Sciences

- Facebook
- Twitter
- Linkedin
- Whatsapp
Altmetric

Table of Contents

Results
Discussion and Conclusion
Materials and Methods
The 18F-FDG PET Image Analysis and Computation of the Relevant Metrics
Supplementary Material
Data Availability.
Change History

Abstract

Through the use of different in silico modeling approaches capturing tumor heterogeneity, we postulated that areas of high metabolic activity would shift toward the periphery as tumors become more aggressive. To confirm the hypothesis and provide clinical value for the finding, we collected ¹⁸F-FDG PET images of breast cancers and non–small-cell lung cancers, where we measured the distance from the point of maximum activity to the tumor centroid, normalizing it by a surrogate of the volume. The metric, NHOC, showed higher prognostic value than other classical PET-based metabolic biomarkers used in oncology, evidencing that the shift of the hotspot of activity from the center of the tumor to its periphery correlates with a poor prognosis.

Human cancers are biologically and morphologically heterogeneous. A variety of clonal populations emerge within these neoplasms and their interaction leads to complex spatiotemporal dynamics during tumor growth. We studied the reshaping of metabolic activity in human cancers by means of continuous and discrete mathematical models and matched the results to positron emission tomography (PET) imaging data. Our models revealed that the location of increasingly active proliferative cellular spots progressively drifted from the center of the tumor to the periphery, as a result of the competition between gradually more aggressive phenotypes. This computational finding led to the development of a metric, normalized distance from ¹⁸F-fluorodeoxyglucose (¹⁸F-FDG) hotspot to centroid (NHOC), based on the separation from the location of the activity (proliferation) hotspot to the tumor centroid. The NHOC metric can be computed for patients using ¹⁸F-FDG PET–computed tomography (PET/CT) images where the voxel of maximum uptake (standardized uptake value [SUV]_max) is taken as the activity hotspot. Two datasets of ¹⁸F-FDG PET/CT images were collected, one from 61 breast cancer patients and another from 161 non–small-cell lung cancer patients. In both cohorts, survival analyses were carried out for the NHOC and for other classical PET/CT-based biomarkers, finding that the former had a high prognostic value, outperforming the latter. In summary, our work offers additional insights into the evolutionary mechanisms behind tumor progression, provides a different PET/CT-based biomarker, and reveals that an activity hotspot closer to the tumor periphery is associated to a worst patient outcome.

Keywords

Cancer, 18F-FDG PET /CT, evolutionary dynamics, prognostic biomarker

License: CC BY

Jiménez-Sánchez,Bosque,Jiménez Londoño,Molina-García,Martínez,Pérez-Beteta,Ortega-Sabater,Honguero Martínez,García Vicente,Calvo,and Pérez-García: Evolutionary dynamics at the tumor edge reveal metabolic imaging biomarkers

Human cancers are genetically and morphologically heterogeneous (1, 2). This is generally attributed to the evolutionary dynamics of different clonal cell populations coexisting in the tumor ecosystem and undergoing stochastic branching processes over time (3 4–5). Successively acquired driver mutations, somatic alterations, and nongenetic modifications may confer increased fitness on certain cancer cell phenotypes, which subsequently outcompete those that do not experience such selection benefits within their microenvironment (4, 6, 7). Cells with specific advantageous traits may not show uniform spatial distribution across the tumor, particularly in large tumors. In fact, trade-offs exist that preclude the occurrence of optimal phenotypes, as exemplified by the hallmarks of cancer (8), and thus only local selection is expected to take place. This produces the spatial phenotypic diversity found in primary tumors and distant metastases (9).

Sustained metabolic reorganization during tumor progression, due to bioenergetically very demanding processes such as rapid proliferation, is a major hallmark of cancer (8, 10). This gives rise to a global metabolic plasticity and fitness optimization that confers evolutionary advantages under specific selective pressures, such as hypoxia (11). Positron emission tomography (PET) has been proposed as a way to assess macroscopic tumor heterogeneity in human patients (12). The technique is used in clinical practice with the radiotracer ¹⁸F-fluorodeoxyglucose (¹⁸F-FDG) (13), which is an analog of glucose and thus a marker of glycolysis (14). The altered tumor metabolism leads to an up-regulation of glycolysis and an increase in glucose consumption (15). This happens even in the presence of oxygen and is referred to as the Warburg effect. Even though this process is energetically inefficient (16), cancer cells may find it beneficial to satisfy the biomass demands required by their high proliferation rates (17). This is confirmed by studies that relate the uptake of ¹⁸F-FDG in PET images to proliferation markers (18, 19). Therefore, the spatial map of glucose consumption provided by ¹⁸F-FDG PET images, as measured at each voxel by the standardized uptake value (SUV), is of great utility in portraying the spatial distribution of proliferation within the tumor.

The degree and impact of intertumor diversity and intratumor heterogeneity in patients has driven the need for quantitative frameworks to account for this variability (20). We considered how the metabolic activity might be distributed inside the tumor and how that information could be related to ¹⁸F-FDG PET images. Specifically, we looked at whether the location of prominent proliferation hotspots, as measured by the voxel of maximum radiotracer uptake (SUV_max), could convey information about patient prognosis. We did this by analyzing how these spots changed over time and space within the tumor in silico using two mathematical models of different levels of complexity. The implications of these results, embodied in the definition of a prognostic biomarker named NHOC, were tested on datasets of breast and lung cancer patients.

Results

Phenotype Variability Supports a Drift of the Highest Metabolic Activity toward the Tumor Boundary.

To describe the emergence of metabolic heterogeneity, we studied in silico a simple biological scenario assuming the tumor to be composed of a clonal population of cells that can migrate, proliferate until the physical space is full, and die. To account for phenotypic heterogeneity, a transition probability that a cell proliferating at a rate $ρ$ could increase or decrease its rate was introduced. The mathematical model used was a continuous nonlocal Fisher–Kolmogorov-type equation (21) which considered the tumor cell population to be structured both by a spatial position vector $x \in Ω \subset R^{3}$ , inside a domain $Ω$ , and a proliferation rate $ρ \in [0, ρ_{m}]$ , where $ρ_{m}$ is a maximum proliferation rate. Let $u = u (x, ρ, t)$ denote the cell density function, such that $u (x, ρ, t) d^{3} x d ρ$ represents the number of tumor cells that, at time $t$ , have a proliferation rate $ρ$ at point $x$ . We modeled the dynamics of $u (x, ρ, t)$ via the following migration–proliferation integro-differential equation:

\begin{align} \frac{\partial u}{\partial t} & = D_{c} \nabla^{2} u + D_{ρ} \frac{\partial^{2} u}{\partial ρ^{2}} \\ + (ρ - μ) (1 - \frac{1}{K} \int_{0}^{ρ_{m}} u (x, ρ^{'}, t) d ρ^{'}) u (x, ρ, t) . \end{align}

[1]

The first term accounts for cell migration with a diffusion constant

D_{c} > 0

. The second term captures the effect of nongenetic instability, mediated by fluctuations in the proliferation phenotype occurring with a diffusion constant

D_{ρ} > 0

. Note that the proliferation phenotype is a hallmark in tumors resulting from alterations in growth regulation (). The third term comprises two main factors. The first one includes the proliferation rate

ρ

minus a constant death rate

μ > 0

; those cells having a larger factor

ρ - μ

will tend to display a fitness advantage unless exogenous mechanisms (e.g., cytotoxic drugs targeting actively dividing cells) exert a negative selection effect on the phenotype. The second factor consists of a nonlocal logistic form with a carrying capacity

K > 0

. This factor represents the interplay between intratumor subpopulations with different proliferations competing for the available space.

A number of quantities are useful for summarizing the information contained in Eq. 1. The first one is the marginal cell density $n (x, t) = \int_{0}^{ρ_{m}} u (x, ρ, t) d ρ$ , with typical radially symmetric profiles as shown in Fig. 1A. The second one is the proliferation density $M (x, t)$ (Eq. 9 in Materials and Methods), which gives the spatiotemporal proliferation map and allows the tumor regions with high metabolic activity to be identified. Fig. 1B depicts $M (x, t) / K$ and shows how the location of the highest activity shifts from the tumor centroid toward the boundary as it grows in silico. This observed displacement, which was found to be linear with time, was quantified using two metrics. The first one was the distance from the highest activity, corresponding to the point of maximum proliferation, to the tumor centroid. We named this metric the distance from the metabolic hotspot to the tumor centroid (HOC). To make HOC independent of size, so that it can be compared among different tumors, we conceived a second metric, normalized HOC (NHOC), defined as the ratio between HOC and the mean metabolic radius of the tumor, $R_{m e t}$ (Fig. 1C and Eq. 10 in Materials and Methods). Simulations of Eq. 1 showed that, during the early stages of the natural history of the tumor, the metric HOC was found to be zero or very small. However, as the inner regions were filled with cells, HOC increased linearly with time (Fig. 1D). NHOC changed steadily from zero to one, since the maximum proliferation spot can only occur between the tumor center and its edge, indicating that this spot will move from the central regions of the tumor to its boundaries (Fig. 1 D, Inset).

Fig. 1.

Nonlocal Fisher–Kolmogorov model 1 predicts a drift of the highest metabolic activity from the tumor centroid to the periphery with time. (A) Normalized cell density $n (x, t_{j}) / K$ at $t_{j} =$ 6, 12, 18, 24, 30, and 36 mo (from left to right) for a radially symmetric tumor. (B) Pseudocolor plots of the normalized spatiotemporal proliferation density $M (x, t)$ and profiles (Inset) of $M (x, t_{j})$ calculated at $t_{j} =$ 6, 12, 18, 24, 30, and 36 mo. (C) Mean metabolic radius $R_{m e t} (t)$ and (Inset) average proliferation rate $ρ_{a} (t)$ . (D) Variation over time of the distance from the tumor centroid to the hotspot of proliferation (HOC) and (Inset) normalized HOC by the mean metabolic radius (NHOC). Simulation parameters are listed in Materials and Methods.

The simulations of Eq. 1 revealed other noteworthy effects. First, the amplitude of the maximum activity in $M (x, t)$ grew with time, meaning that tumors at later stages of their evolutionary history show increasing hotspot activity values. This was in line with the frequently observed association between SUV_max and prognosis for different tumor histologies (22, 23). Second, the distribution of the proliferation rates displayed sustained growth toward higher values of $ρ$ , reflected in the average proliferation rate $ρ_{a}$ (Fig. 1 C, Inset). The growth of the tumor proliferation rate with time (size) has been experimentally observed in other studies (24).

Genotype Evolutionary Dynamics Support the Drift of Tumor Highest Metabolic Activity toward the Boundary.

We next resorted to a more complex and realistic biological scenario accounting for genotypic alterations. We did this by considering a stochastic discrete model based on fundamental cell features. At the cellular level, cancer cells can be characterized by four deregulated processes: proliferation, migration, mutation, and death. These processes can be easily implemented as rules in a discrete mathematical model to mimic the main characteristics of the real system, with the drawback of facing high computational cost, especially when simulating clinically relevant volumes (25). To overcome this problem, we developed a hybrid stochastic mesoscale model of tumor growth that allowed clinically relevant tumor sizes to be simulated while retaining the basic cancer hallmarks (24).

The model was parameterized for two of the most prominent cancer types, namely breast and lung cancer (non–small-cell lung carcinoma [NSCLC]). A summary of our available data can be seen in Table 1. Mutational landscapes were constructed based on a simplification of their known mutational spectra. Alterations in EGFR and ALK, which are strongly associated with nonsquamous lung adenocarcinoma, were considered to model NSCLC, while driver mutations in PIK3CA and TP53 were considered for breast cancer (26 27 28–29). Therefore, the mutational tree in both types of tumors simulated had two possible altered genes, leading to four possible combinations or “genotypes” that define four different clonal populations. Basal rates associated a characteristic time to each basic cell process, and mutation weights determined how these basal rates were affected once a given alteration was acquired. Mutation weights were taken to contribute equally for all alterations, and their effect was cumulative, so that a cell carrying two alterations simultaneously would perform basic processes with a double advantage. Thus, the stochastic mesoscopic model provided a richer scenario to explore intratumoral heterogeneity during tumor growth.

Table 1.

Stochastic model parameters

Parameter	Meaning	Value	Ref.
$L$	No. of voxels per side	80	—
$Δ x$	Voxel side length (mm)	1	(62)
$Δ t$	Time-step length (h)	24	—
$K$	Carrying capacity per voxel (cells)	$2 \cdot 1 0^{5}$	(63)
$N_{t}$	Threshold cell number	0.2 $\cdot$ K
$N_{0}$	Initial population (cells)	1	—
$V_{e n d}$	Maximum reachable tumor volume ( ${c m}^{3}$ )	50 (breast)
		120 (lung)	(64)
$V_{diag}$	Tumor volume at diagnosis ( ${c m}^{3}$ )	0.3 to 5 (breast)	(65)
		0.2 to 15 (lung)	(66)

We ran 100 simulations of breast cancer and 100 simulations of NSCLC with random parameters uniformly sampled from the ranges in Tables 1 and 2 (Materials and Methods). Cell number, activity (number of newborn cells), and most abundant clonal population were calculated for each voxel and time step (Fig. 2). Tumor volumes were measured from the number of voxels containing more than a threshold number of cells $N_{t}$ (Table 1), and the mean spherical radii (MSR) were computed from these.

Table 2.

Basal rates and mutation weights

Processes	Proliferation	Migration	Death	Mutation
Breast basal rates ( $d^{- 1}$ )	0.0133 to 0.04	0.02 to 0.0303	0.01 to 0.02	0.01 to 0.02
TP53/PIK3CA (%)	20 to 40	25 to 45	(−30) to (−10)	25 to 40
Lung basal rates ( $d^{- 1}$ )	0.04 to 0.2	0.03 to 0.2	0.03 to 0.1205	0.03 to 0.1205
TP53/KRAS (%)	20 to 40	25 to 45	(−30) to (−10)	25 to 40

Fig. 2.

Hybrid stochastic mesoscale model shows that competition among progressively more aggressive phenotypes is pushed to the edge. (A) Three-dimensional volume renderings at different time frames (from left to right: 78, 85, and 92% of simulation) of a simulation of breast cancer growth depicting clonal populations within the tumor. Color of cell populations ranges from green (less aggressive) to red (more aggressive). Rates are proliferation 0.0315 $d^{- 1}$ , death 0.0157 $d^{- 1}$ , mutation 0.0160 $d^{- 1}$ , and migration 0.0235 $d^{- 1}$ . (B) Central section for the same simulation and time frames as in A showing the most abundant clonal population per voxel. (C) Central section of tumor activity for the same time frames as in A. (D and E) HOC progression for every simulation of breast cancer (D) and NSCLC (E). (F and G) Longitudinal NHOC dynamics for simulations of breast cancer (F) and NSCLC (G) growth, with individual runs colored in reddish orange and all-simulation averaged NHOC in blue; crosses depict the time points at which each simulation ended.

As cells mutated in silico, new clonal populations emerged with higher, more advantageous migration and proliferation rates. These new clones increased their relative abundance in the tumor, eventually becoming fixed in the system. As the tumors grew larger, cell division occurred preferentially at the tumor periphery. This was as expected, since inner voxels became progressively filled with cells that prevented them from proliferating. Voxels where the most aggressive clonal population was more abundant were associated with hotspots of maximum proliferation. Therefore, evolution was pushed toward the tumor edge: Cells with higher fitness (especially those having higher proliferation rates) appeared farther from the tumor center as they grew. At each time step, the maximum proliferation spot was identified as the voxel with the largest number of cell births, and its distance to the tumor centroid (HOC) was calculated. Fig. 2 D and E shows a monotonic increase of HOC with time for both histologies. Normalizing with respect to the MSR to get the NHOC showed that the point of maximum proliferation was displaced toward the boundary (Fig. 2 F and G) in all of the simulations performed. The only difference between simulations was the time that the maximum proliferation spot took to reach the edge. Thus, NHOC was predicted to be a robust property related to the evolutionary state of the disease.

PET Imaging Data Confirm Evolutionary Dynamics of the Maximum Metabolic Activity Validating Related Biomarkers.

The computational results suggested that NHOC might contain meaningful prognostic information. In the clinical setting, the metabolic activity distribution of the tumor can be evaluated by means of ¹⁸F-FDG PET, which reflects the biological processes taking place at a lower level (30) and is frequently used on newly diagnosed breast cancer and NSCLC patients. To confirm or refute the theoretical predictions, we performed a study on our patient cohort (Materials and Methods). For each patient, the tumor was delineated in the images and the locations of centroids and SUV_max were obtained computationally from the segmented distribution, as detailed in Materials and Methods. The metabolic tumor volume (MTV), total lesion glycolysis (TLG) (integral of the SUV distribution over the volume), and NHOC metrics were calculated for all tumors for both histologies. Two typical examples of ¹⁸F-FDG PET images from breast cancer patients are shown in Fig. 3 A and B, respectively. Small values of NHOC, with SUV_max close to the tumor centroid as in Fig. 3 B and I, were expected to correspond to less developed disease, in accordance with the previous theoretical framework. In contrast, the cases shown in Fig. 3 C and J with SUV_max displaced in relation to the centroid would correspond to tumors with a poorer prognosis.

Analysis of NHOC in 18F-FDG PET images reveals well-behaved distributions with prognostic potential. (A and H) Examples of PET images for breast cancer (A) and NSCLC (H) patients in our dataset. (B, C, I, and J) Two-dimensional slices from patients with small (B and I) and large (C and J) NHOC values for breast cancer (B and C) and NSCLC (I and J) patients. The centroid of each segmented lesion and the voxel of SUVmax are marked with a white cross and a green dot, respectively. (Scale bars, 1 cm.) (D–G and K–N) Histograms showing the distributions of metabolic tumor volume (D and K), total lesion glycolysis (E and L), SUVmax (F and M), and NHOC (G and N) for breast cancer (D–G) and NSCLC (K–N) patients in our datasets.

Fig. 3.

Analysis of NHOC in ¹⁸F-FDG PET images reveals well-behaved distributions with prognostic potential. (A and H) Examples of PET images for breast cancer (A) and NSCLC (H) patients in our dataset. (B, C, I, and J) Two-dimensional slices from patients with small (B and I) and large (C and J) NHOC values for breast cancer (B and C) and NSCLC (I and J) patients. The centroid of each segmented lesion and the voxel of SUV_max are marked with a white cross and a green dot, respectively. (Scale bars, 1 cm.) (D–G and K–N) Histograms showing the distributions of metabolic tumor volume (D and K), total lesion glycolysis (E and L), SUV_max (F and M), and NHOC (G and N) for breast cancer (D–G) and NSCLC (K–N) patients in our datasets.

The histograms in Fig. 3 D–G and K–N depict the distributions of MTV, TLG, SUV_max, and NHOC for both histologies. It is noteworthy that the NHOC has a more regular distribution than the other PET-based measures, with definite values between 0 and 1 and a centered mean (breast cancer, $0.51 \pm 0.18$ , median 0.50; NSCLC, $0.43 \pm 0.2$ , median 0.39). It is clear from Fig. 3 G and N that at the time of diagnosis the point of maximum uptake is typically located away from the geometrical center of the tumor. To inspect the relationships between the PET variables, we calculated a correlation matrix analysis (SI Appendix, Fig. S1). The results show that NHOC does not strongly correlate with the conventional measures and can therefore be considered as an independent metric.

The classical measures (MTV, TLG, SUV_max) are known to be prognostic biomarkers in breast cancer and NSCLC (22, 23). For these variables we performed Kaplan–Meier analyses on overall survival (OS) and disease-free survival (DFS) (SI Appendix, Figs. S3 and S4). All of the variables had prognostic value in the breast cancer cohort, but only MTV returned significant results (P value < 0.05) in the NSCLC cohort.

We then tested the prognostic value of NHOC by Kaplan–Meier analyses with OS and DFS as endpoints (Materials and Methods). Results for the best splitting thresholds are shown in Fig. 4. For the breast cancer cohort, NHOC showed robust results in terms of OS, with a best splitting threshold in both OS and DFS of NHOC = 0.499 (Fig. 4 A and C). Interestingly for OS, the most relevant metric, the C index, reached an outstanding value of 1 (for DFS it was 0.899). Thus, no patients with tumors having their SUV_max closer than half the radius ( $n$ = 30) died from the disease. In NSCLC, NHOC separated the patients well, and the best splitting threshold, NHOC = 0.64, led to a C index of 0.875 for OS. The separation in median OS between groups was 57.33 mo, while in DFS it was 36.62 mo. Therefore, patients for whom the ¹⁸F-FDG hotspot displays an increasing shift from the tumor centroid are associated to a worst outcome.

Fig. 4.

Kaplan–Meier curves obtained for best splitting thresholds corresponding to the NHOC metric. (A) Overall survival in breast cancer cohort. (B) Overall survival in NSCLC cohort. (C) Disease-free survival in breast cancer cohort. (D) Disease-free survival in NSCLC cohort.

To determine whether the prognostic value of NHOC might be related to possible necrotic regions, we evaluated the presence of hypometabolic voxels inside the delineated tumor, i.e., those situated inside the tumor that have a lower SUV value than the segmentation threshold. We found that these regions are present only in 7 of the 61 patients of our breast cancer cohort and in 13 of the 161 patients of our NSCLC cohort. Moreover, tumors where hypometabolic voxels are detected account on average for only 1.63% of the total voxels in the breast cancer cohort and 2.68% of the voxels in the NSCLC cohort. Consequently, in our cohorts, the NHOC significance is not associated to the occurrence of necrotic regions assessed by the manifestation of hypometabolic voxels in the images.

The presence of lymph node metastases is a strong predictor of outcome in breast and lung cancers, with nodal metastases having negative prognostic significance. We evaluated whether there was a significant difference in the NHOC value for the group of patients which showed nodal metastases on diagnosis and those who did not (Mann–Whitney test; SI Appendix, Fig. S9). The P values were 0.4683 for the breast cancer dataset and 0.1071 for the NSCLC, thus discarding any significant difference across groups and indicating that the information carried by NHOC is independent of the existence of lymph node metastases.

These results show the strength of NHOC as a prognostic biomarker in comparison with the classical metrics. For OS in the breast cancer cohort, only MTV approached the performance of NHOC; however, the C index for NHOC (C index = 1) outperformed the result for MTV (C index = 0.875). For DFS, none of the classical variables showed nonisolated thresholds leading to a statistically significant association between subgroups. In the lung cancer cohort, the NHOC metric outperformed, once again, the prognostic value of the classical variables. Regarding OS, there were ranges of thresholds of MTV, TLG, and SUV_max leading to statistically significant results with best C indexes of 0.736 (MTV), 0.682 (TLG), and 0.658 (SUV_max) still substantially lower than the value obtained for NHOC (0.875). Results for DFS were again similar, with only MTV and TLG achieving significance, with best values of 0.638 (MTV) and 0.607 (TLG), but still underperforming NHOC, with a C index of 0.651.

Discussion and Conclusion

Heterogeneity is one of the hallmarks of tumor malignancy (1, 2). Many mathematical models have been constructed accounting for different aspects of the development of heterogeneity through evolutionary dynamical processes in a number of cancer types (31 32–33). We did not intend here to develop a universal mathematical model to describe every aspect of tumor growth progression, but rather to focus on understanding the basic dynamics of the hotspot of metabolic activity, due to the potential applicability of the results. Different levels of complexity were considered in each of the two complementary models constructed, and both led to the same conclusions: 1) Tumors would evolve toward higher proliferation rate values, and 2) maximum metabolic activity would move toward the tumor edge as the tumor evolves with time.

It is interesting to note that conclusion 2 (but not conclusion 1) would have been arrived at through the use of a “classical” local Fisher–Kolmogorov model (34). In the context of that mathematical model, proliferation is inhibited in areas of high cell density and higher proliferation areas would switch from the tumor core to the periphery, as predicted also by the nonlocal model Eq. 1. This is, therefore, a robust finding of the study, and it is likely to be anticipated by other mathematical approaches. Conclusion 1 does not come as a surprise, since the maximum metabolic activity obtained from PET images (e.g., as measured by the SUV_max) has been known to contain prognostic information in different cancer types (22, 23). Thus, the fact that ${S U V}_{max}$ possesses a prognostic value in the clinical setting is compatible with the results of our models, where the activity grows with time over the tumor’s natural history until a maximum value is reached.

Many studies have correlated either classical PET-derived metrics such as MTV and TLG (35, 36) or complex spatial features of the distribution of SUV values (37) with the outcome of the disease. However, no study has analyzed the prognostic value of metrics derived from the location of the highest metabolic activity. The fact that this simple biomarker has a high prognostic value is remarkable and probably related to the robustness of the biological assumptions behind the mathematical models used to substantiate it. In fact, it is natural to expect that the presence of more aggressive glucose-avid cells, that might be unable to progress when located in saturated areas near the center of the tumor, may be a risk factor when placed in regions with much more capacity to settle and invade.

In our study we chose to take the voxel bearing the SUV_max as the location of maximum metabolic activity to compute NHOC. We could also have used SUV_peak (the maximum SUV appearing in the distribution when all of the voxels are averaged with their 26 neighbors), which is thought to be more stable and to better define an extended region of high uptake (38). However, SUV_max is often placed in the area defined by SUV_peak, thus leading to equivalent metrics (39). SUV_max is also easier to identify visually and is therefore easier to use in clinical practice, besides being the simplest option.

The fact that NHOC provides an accessible and powerful prognostic metric could be extended in different ways. First, it would be valuable to look at whether changes in this biomarker might provide a robust indication of an increase in malignancy for initially indolent tumors (e.g., benign lung nodules, low-grade gliomas, etc.) undergoing a malignant transformation. Second, an intriguing open question would be to determine whether the rate at which NHOC changes during patient follow-up correlates with the occurrence and fixation of specific mutations. Finally, it may be the case that changes with time of this metric, after different treatment modalities, could help in assessing the response through sequential PET studies as a measure of how much NHOC is reduced.

Hotspots in staging ¹⁸F-FDG PET–computed tomography (PET/CT) have previously been reported as preferential sites for relapse after chemoradiotherapy in several types of cancer, including NSCLC (40, 41). Some authors have found that these hotspots colocalize with hypoxic regions as inferred from HIF-1 $α$ tumor immunostaining (42) and are thus very relevant in the context of radioresistance (40, 43, 44). Furthermore, these hotspots have also been linked to the occurrence of somatic mutations in lung cancer (45, 46) similarly to what is described in this work. Although further research is needed, this rationale has led some authors to assert that patients may benefit from hotspot-based dose escalation (43, 47). The voxel of SUV_max, hence, has a preeminent significance both in prognosis and in therapy and, as it is now shown here, its relative position has a direct bearing on its relevance. In addition to NHOC tracking throughout the disease, NHOC estimation in hypoxia-specific PET imaging, employing other radiotracers such as fluoromisonidazole (¹⁸F-MISO) or fluoroazomycin arabinoside (¹⁸F-FAZA), could provide valuable information complementary to that of ¹⁸F-FDG PET/CT images.

Mathematical and computational models are progressively gaining their place among the tools that are used to study cancer (48). In silico models based on evolutionary dynamics may capture relevant aspects of tumor growth and have proved helpful in understanding tumor clonal heterogeneity, one of the main hallmarks of cancer (49). Mechanistic mathematical models of different levels of complexity have been shown to provide biomarkers of clinical significance (24, 50 51 52 53 54 55 56–57). This type of approach provides a rational alternative to radiomic and deep-learning studies, where a mechanistic explanation is often missing. The study described in this paper falls into the former category, demonstrating that an informed understanding of the system’s emergent properties can shed light on the deeper roots of its working.

It is worth mentioning that our mathematical approach, beyond its fundamental interest, has led to the proposal of a simple metric with high prognostic value that can be obtained from ¹⁸F-FDG PET studies. The NHOC biomarker was able to separate patients with breast cancer and non–small-cell lung cancer into two groups with significantly different survival (both overall and disease-free) and proved to be more powerful than traditional ¹⁸F-FDG PET/CT biomarkers (MTV, TLG, SUV_max) currently used in clinics. This demonstrated that the geometric location of the maximum metabolic activity, and not only its value, contains information of clinical significance: Tumors in which the hotspot is located farther from the center are found to have a worst prognosis.

This study opens many further avenues for research. The first one is the search for other biomarker definitions accounting for the location of highest metabolic activity. Second, it would be interesting to test our findings in other tumor histologies, such as lymphoma or melanoma (58, 59). PET is a mainstream technique, increasingly employed in clinics and in many imaging studies for which a broad spectrum of tumor histologies is available. Thus, the applicability of NHOC to other cancer types would be an interesting extension of our work.

In conclusion, by using two mathematical models incorporating evolutionary dynamics, we have shown that the maximum metabolic activity is expected to increase in magnitude and to move toward the tumor boundary as human solid tumors progress. On the basis of the theoretical predictions we have defined a metric, the NHOC, representing the normalized distance from the point of maximum activity to the tumor centroid and validated it as a prognostic biomarker in lung and breast cancer patients using PET imaging datasets. The biomarker NHOC outperformed classical PET-based biomarkers such as TLG, MTV, and SUV_max and provides a notable example of mathematically grounded research with applicability in oncology.

Materials and Methods

Patients.

Our study was based on data from two different studies, both of them approved by the Institutional Review Board (IRB) of Ciudad Real University General Hospital (HGUCR). Breast cancer patients were participants of a multicenter prospective study, and written informed consent was obtained from all patients. The inclusion criteria were 1) newly diagnosed locally advanced breast cancer with clinical indication of neoadjuvant chemotherapy, 2) lesion uptake higher than background (i.e., those having a SUV_max larger than twice the background activity readings), 3) absence of distant metastases confirmed by other methods prior to the request of PET/CT for staging, and 4) breast lesion size of at least 2 cm. Sixty-one patients (18% lobular carcinoma, 82% ductal carcinomas, 100% women, age rank 25 to 80 y, median 50 y) were included in this dataset. The TNM data were 54% T2, 18% T3, 28% T4; 28% N0, 55% N1, 6% N2, 11% N3; 100% M0.

One hundred seventy-five patients (153 men, 22 women, age rank 41 to 84 y, median 65 y) were included in the study from a dataset of lung cancer patients who received surgery in the period June 2007 to December 2016. Histologies were 63 squamous-cell carcinomas and 112 adenocarcinomas. Staging information was 69 stage I, 70 stage II, 33 stage III, 3 stage IV. The N staging was 107 patients N0, 46 N1, and 22 N2. All patients had M0. PET protocol and machine were as in subgroup 1. The inclusion criterion was established that minimal lesion size should be greater than 2.0 cm. From those initial patients, 14 were removed due to the unavailability of survival data.

The PET machine was a dedicated whole-body PET/CT scanner (Discovery SDTE-16s; GE Medical Systems) in three-dimensional (3D) mode. Image acquisition began 60 min after intravenous administration of approximately 370 MBq (10 mCi) of ¹⁸F-FDG; the images obtained had a voxel size of 5.47 $\times$ 5.47 $\times$ 3.27 ${m m}^{3}$ , with no gap between slices, and a matrix size of $128 \times 128$ . The inclusion criteria considered only newly diagnosed patients with availability of pretreatment PET/CT examination and a lesion uptake higher than background (SUV_max larger than twice the background), absence of distant metastases, and a lesion size of at least 2 cm.

The ¹⁸F-FDG PET Image Analysis and Computation of the Relevant Metrics

PET images in DICOM format were loaded into MATLAB for the image analysis. In each image, the tumor was manually selected and subsequently delineated in 3D by an automatic algorithm. The result of the delineation is a 3D matrix with the SUV values ( $S$ ) of the $N$ voxels of the tumor. Occasionally, there might be $N^{'}$ hypometabolic voxels enclosed by the active ones, which could be assigned to necrotic regions; these are used here only for geometric considerations. For all delineated tumors in PET images, we evaluated the following metrics:

MTV, i.e., volume of the delineated tumor, computed as the number $N$ of selected voxels multiplied by the volume of one voxel $V_{V}$ :

M T V = N \times V_{V} .

[2]

TLG, calculated as the sum of the SUV value multiplied by the volume of the voxel for all of the N voxels in the tumor:

T L G = \sum_{i = 1}^{N} S_{i} \times V_{V} .

[3]

SUV_max, maximum value of SUV in the tumor:

S U V_{m a x} = max {S_{i}} .

[4]

We refer to the SUV_max voxel position by its coordinates $x_{s m}$ , $y_{s m}$ , $z_{s m}$ and use it as a location for the ¹⁸F-FDG hotspot.

MSR, the radius of a hypothetical sphere having the same volume as the MTV and serving as a linear surrogate of volume:

M S R = {(\frac{3}{4 π} M T V)}^{1 / 3} .

[5]

Centroid: As a reference point for the center of the tumor that is independent of the observer and computable for any delineated tumor, we use the geometrical centroid of the 3D shape defined by the segmented tumor and all its interior points, including the $N^{'}$ hypometabolic ones. Its coordinates are computed by the mathematical definition

x_{c} = \frac{1}{N + N^{'}} \sum_{i = 1}^{N + N^{'}} x_{i},

[6a]

y_{c} = \frac{1}{N + N^{'}} \sum_{i = 1}^{N + N^{'}} y_{i},

[6b]

z_{c} = \frac{1}{N + N^{'}} \sum_{i = 1}^{N + N^{'}} z_{i} .

[6c]

HOC, distance from ¹⁸F-FDG hotspot (voxel with SUV_max) to tumor centroid measured as the Euclidean distance between both points:

H O C = \sqrt{{(x_{s m} - x_{c})}^{2} + {(y_{s m} - y_{c})}^{2} + {(z_{s m} - z_{c})}^{2}} .

[7]

NHOC, normalized distance from ¹⁸F-FDG hotspot to tumor centroid. To make the HOC size independent, we normalize it by the MSR as a linear measure of volume. Consequently, we get a metric for the shift of the ¹⁸F-FDG hotspot that is comparable along all tumors:

N H O C = \frac{H O C}{M S R} .

[8]

Even though all of the images were taken with PET/CT coregistration, only data from positron emission tomography were used in this work.

Kaplan–Meier Statistics.

We performed Kaplan–Meier analyses over these two cohorts of patients, using the log-rank and Breslow tests to assess the significance of the results. These methods compare two populations separated in terms of one parameter and study their statistical differences in survival. Specifically, OS and DFS Kaplan–Meier analyses were performed. A two-tailed significance level with P value lower than 0.05 was applied. The hazard ratio (HR) and its adjusted 95% confidence interval (CI) were also computed for each threshold using Cox proportional hazards regression analysis.

Splitting Thresholds.

For each variable, we searched for every value splitting the sample into two different subgroups, satisfying the condition that none of them be more than five times larger than the other. We then tested each of them as splitting thresholds through Kaplan–Meier analyses, obtaining the significance results shown in SI Appendix, Figs. S5–S8. The best splitting threshold was chosen as the nonisolated significant value giving the lowest P value in both log-rank and Breslow tests, as described in ref. 52.

Harrell’s C Index.

To assess the accuracy of prognostic models, Harrell’s concordance index score was also computed (60). This method compares the survival of two populations of patients (best prognosis versus worst prognosis) by studying all possible combinations of individuals belonging to different groups. Then, the percentage of right guesses is the reported result. Concordance indexes were computed using the noncensored sample and ranged from 0 to 1, with 1 indicating a perfect model (a purely random guess would give a concordance index of 0.5).

Variable Correlations.

Spearman correlation coefficients were used to assess the dependencies between pairs of variables. We considered significant correlation coefficients above 0.7 or below −0.7 as strong (direct or inverse, respectively) correlations between variables. In this way we were able to exclude possible confounding effects in our analysis.

Statistical Software.

SPSS (v. 22.0.00), MATLAB (R2018b; The MathWorks, Inc.), and R (3.6.3) software were used for all statistical analyses.

Nonlocal Fisher–Kolmogorov Model and Simulations.

The migration–proliferation integro-differential Eq. 1 in radial coordinates was solved numerically using the method of lines (61) combined with Newton–Cotes integration formulas to deal with the nonlocal term. In the simulations displayed in Fig. 1, the computational domain consisted of a radial variable $r \in [0, R_{max}]$ , where $R_{max} = 7$ cm was the maximum radius, and the proliferation rate $ρ \in [0, ρ_{m}]$ , where the maximum proliferation rate was $ρ_{m} = 0.06 d^{- 1}$ . The number of nodes in the discretized $r$ - $ρ$ mesh was $350 \times 180$ . Additional parameters were $D_{c} = 3.5 \cdot 1 0^{- 4} {c m}^{2}$ /d, $D_{ρ} = 1.3 \cdot 1 0^{- 8} d^{- 3}$ , $μ = 4 \cdot 1 0^{- 3} d^{- 1}$ , and $K = 8 \cdot 1 0^{7}$ cells/ ${c m}^{3}$ . The initial condition consisted of a highly localized lesion with a radius of 1 mm containing $1 0^{5}$ tumor cells and having a mean proliferation rate $ρ_{0} = 1.7 \cdot 1 0^{- 2} d^{- 1}$ and standard deviation $σ_{0} = 3 \cdot 1 0^{- 3} d^{- 1}$ .

The general expression for the proliferation density is

\begin{align} M (x, t) = \int_{0}^{ρ_{m}} (ρ - μ) (1 - \frac{1}{K} \int_{0}^{ρ_{m}} u (x, ρ^{'}, t) d ρ^{'}) u (x, ρ, t) d ρ \end{align}

[9]

and was used to compute the plots shown in Fig. 1B assuming spherical symmetry.

In Fig. 1C, the mean metabolic radius was defined as

R_{m e t} (t) = \frac{\int_{0}^{R_{max}} M (r, t) r^{3} d r}{\int_{0}^{R_{max}} M (r, t) r^{2} d r},

[10]

while the average proliferation rate was determined via

ρ_{a} (t) = \frac{\int_{0}^{ρ_{m}} \int_{0}^{R_{max}} ρ u (r, ρ, t) r^{2} d r d ρ}{\int_{0}^{ρ_{m}} \int_{0}^{R_{max}} u (r, ρ, t) r^{2} d r d ρ} .

[11]

In Fig. 1D, the distance from the tumor centroid to the point of maximum proliferation (HOC) was calculated at each time step via expression . The NHOC was computed by means of the ratio

N H O C (t) = H O C (t) / R_{m e t} (t)

Acknowledgements

This research has been supported by grants awarded to V.M.P.-G. by the James S. Mc. Donnell Foundation, United States of America, 21st Century Science Initiative in Mathematical and Complex Systems Approaches for Brain Cancer (collaborative award 220020560) and by the Junta de Comunidades de Castilla-La Mancha, Spain (grant number SBPLY/17/180501/000154). V.M.P.-G. and G.F.C. thank the funding from Ministerio de Ciencia e Innovación, Spain (grant number PID2019-110895RB-I00). This research has also been supported by a grant awarded to G.F.C. by the Junta de Comunidades de Castilla-La Mancha, Spain (grant number SBPLY/19/180501/000211). J.J.-S. received support from Universidad de Castilla-La Mancha (grant number 2020-PREDUCLM-15634). J.J.B received support from the University of Castilla-La Mancha (grant number 2018-CPUCLM-7798). A.M. received support from Asociación Pablo Ugarte (http://www.asociacionpablougarte.es). C.O.-S. received support from Asociación Española Contra el Cáncer (grant number 2019-PRED-28372).

The authors declare no competing interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at https://www.pnas.org/lookup/suppl/doi:10.1073/pnas.2018110118/-/DCSupplemental.

Data Availability.

All study data are included in this article and/or SI Appendix. The data extracted by PET imaging processing of all of the patients are accessible in our group’s webpage: http://matematicas.uclm.es/molab/EvolutionaryDynamicsAtTheTumorEdge_DATA.xlsx.

Change History

July 6, 2021: The Acknowledgments have been updated.

References

C. Swanton, Intratumor heterogeneity: Evolution through space and time. Canc. Res. 72, 4875–4882 (2012).

N. McGranahan, C. Swanton, Clonal heterogeneity and tumor evolution: Past, present, and the future. Cell 168, 613–628 (2017).

L. M. Merlo, J. W. Pepper, B. J. Reid, C. C. Maley, Cancer as an evolutionary and ecological process. Nat. Rev. Canc. 6, 924–935 (2006).

M. Greaves, C. C. Maley, Clonal evolution in cancer. Nature 481, 306–313 (2012).

M. C. Lloyd , Darwinian dynamics of intratumoral heterogeneity: Not solely random mutations but also variable environmental selection forces. Canc. Res. 76, 3136–3144 (2016).

Q. Li , Dynamics inside the cancer cell attractor reveal cell heterogeneity, limits of stability, and escape. Proc. Natl. Acad. Sci. U.S.A. 113, 2672–2677 (2016).

A. Álvarez Arenas, A. Podolski-Renic, J. Belmonte-Beitia, M. Pesic, G. F. Calvo, Interplay of Darwinian selection, Lamarckian induction and microvesicle transfer on drug resistance in cancer. Sci. Rep. 9, 9332 (2019).

D. Hanahan, R. A. Weinberg, Hallmarks of cancer: The next generation. Cell 144, 646–674 (2011).

B. Waclaw , A spatial model predicts that dispersal and cell turnover limit intratumour heterogeneity. Nature 525, 261–264 (2015).

D. Jia , Elucidating cancer metabolic plasticity by coupling gene regulation with metabolic pathways. Proc. Natl. Acad. Sci. U.S.A. 116, 3909–3918 (2019).

S. McGuirk, Y. Audet-Delage, J. St.-Pierre, Metabolic fitness and plasticity in cancer progression. Trends Canc. 6, 49–61 (2020).

F. Davnall , Assessment of tumor heterogeneity: An emerging imaging tool for clinical practice?Insights Imaging 3, 573–589 (2012).

W. A. Weber, Positron emission tomography as an imaging biomarker. J. Clin. Oncol. 24, 3282–3292 (2006).

A. Gallamini, C. Zwarthoed, A. Borra, Positron emission tomography (PET) in oncology. Cancers 6, 1821–1889 (2014).

E. Pauwels , FDG accumulation and tumor biology. Nucl. Med. Biol. 25, 317–322 (1998).

R. A. Gatenby, R. J. Gillies, Glycolysis in cancer: A potential target for therapy. Int. J. Biochem. Cell Biol. 39, 1358–1366 (2007).

M. G. V. Heiden, L. C. Cantley, C. B. Thompson, Understanding the Warburg effect: The metabolic requirements of cell proliferation. Science 324, 1029–1033 (2009).

H. Vesselle , Lung cancer proliferation correlates with [F-18]fluorodeoxyglucose uptake by positron emission tomography. Clin. Canc. Res. 6, 3837–3844 (2000).

H. Vesselle , Relationship between non-small cell lung cancer FDG uptake at PET, tumor histology, and Ki-67 proliferation index. J. Thorac. Oncol. 3, 971–978 (2008).

J. Hausser, U. Alon, Tumour heterogeneity and the evolutionary trade-offs of cancer. Nat. Rev. Canc. 20, 247–257 (2020).

F. Hamel, L. Ryzhik, On the nonlocal Fisher-kpp equation: Steady states, spreading speed and global bounds. Nonlinearity 27, 2735–2753 (2014).

R. Bos , Biologic correlates of 18 fluorodeoxyglucose uptake in human breast cancer measured by positron emission tomography. J. Clin. Oncol. 20, 379–387 (2002).

L. F. de Geus-Oei , Biological correlates of FDG uptake in non-small cell lung cancer. Lung Canc. 55, 79–87 (2007).

V. M. Pérez-García , Universal scaling laws rule explosive growth in human cancers. Nat. Phys. 16, 1232–1237 (2020).

J. Metzcar, Y. Wang, R. Heiland, P. Macklin, A review of cell-based computational modeling in cancer biology. JCO Clinical Cancer Informatics 2, 1–13 (2019).

C. Gridelli , Non-small-cell lung cancer. Nat. Rev. Dis. Primers 1, 15009 (2015).

N. Harbeck , Breast cancer. Nat. Rev. Dis. Primers 5, 66 (2019).

T. C. G. A. Network, Comprehensive molecular portraits of human breast tumours. Nature 490, 61–70 (2012).

L. Ding , Somatic mutations affect key pathways in lung adenocarcinoma. Nature 455, 1069–1075 (2008).

M. E. Phelps, Positron emission tomography provides molecular imaging of biological processes. Proc. Natl. Acad. Sci. U.S.A. 97, 9226–9233 (2000).

A. R. Anderson, A. M. Weaver, P. T. Cummings, V. Quaranta, Tumor morphology and phenotypic evolution driven by selective pressure from the microenvironment. Cell 127, 905–915 (2006).

I. Bozic, J. M. Gerold, M. A. Nowak, Quantifying clonal and subclonal passenger mutations in cancer evolution. PLoS Comput. Biol. 12, e1004731 (2016).

C. D. McFarland, L. A. Mirny, K. S. Korolev, Tug-of-war between driver and passenger mutations in cancer and other adaptive processes. Proc. Natl. Acad. Sci. U.S.A. 111, 15138–15143 (2014).

J. D. Murray, Mathematical Biology: I. An Introduction (Springer Science & Business Media, 2007), vol. 17.

S. Liao , Prognostic value of metabolic tumor burden on 18 F-FDG PET in nonsurgical patients with non-small cell lung cancer. Eur. J. Nucl. Med. Mol. Imag. 39, 27–38 (2012).

H. H. Chen, N. T. Chiu, W. C. Su, H. R. Guo, B. F. Lee, Prognostic value of whole-body total lesion glycolysis at pretreatment FDG PET/CT in non–small cell lung cancer. Radiology 264, 559–566 (2012).

T. Berghmans , Primary tumor standardized uptake value (SUVmax) measured on fluorodeoxyglucose positron emission tomography (FDG-PET) is of prognostic value for survival in non-small cell lung cancer (NSCLC): A systematic review and meta-analysis (MA) by the European lung cancer working party for the IASLC lung cancer staging project. J. Thorac. Oncol. 3, 6–12 (2008).

A. Sher , For avid glucose tumors, the SUV peak is the most reliable parameter for [18 F] FDG-PET/CT quantification, regardless of acquisition time. EJNMMI Res. 6, 21 (2016).

L. Ong , Comparison of SUVmax and SUVpeak for response assessment in stage 4 non-small cell lung cancer (NSCLC) by FDG PET/CT. J. Nucl. Med. 56 (suppl. 3), 233 (2015).

H. J. Aerts , Identification of residual metabolic-active areas within individual NSCLC tumours using a pre-radiotherapy 18Fluorodeoxyglucose-PET-CT scan. Radiother. Oncol. 91, 386–392 (2009).

S. F. Petit , Metabolic control probability in tumour subvolumes or how to guide tumour dose redistribution in non-small cell lung cancer (NSCLC): An exploratory clinical study. Radiother. Oncol. 91, 393–398 (2009).

A. van Baardwijk , The maximum uptake of 18f-deoxyglucose on positron emission tomography scan correlates with survival, hypoxia inducible factor-1

α

and glut-1 in non-small cell lung cancer. Eur. J. Canc. 43, 1392–1398 (2007).

J. Calais , Areas of high 18F-FDG uptake on preradiotherapy PET/CT identify preferential sites of local relapse after chemoradiotherapy for non–small cell lung cancer. J. Nuclear Med. 56, 196–203 (2015).

A. Chaput , Correlation between fluorodeoxyglucose hotspots on pretreatment positron emission tomography/CT and preferential sites of local relapse after chemoradiotherapy for head and neck squamous cell carcinoma. Head Neck 39, 1155–1165 (2017).

C. Caicedo , Role of [18F]FDG PET in prediction of KRAS and EGFR mutation status in patients with advanced non-small-cell lung cancer. Eur. J. Nucl. Med. Mol. Imag. 41, 2058–2065 (2014).

S. S. Yip , Associations between somatic mutations and metabolic imaging phenotypes in non-small cell lung cancer. J. Nucl. Med. 58, 569–576 (2017).

H. J. Aerts , Identification of residual metabolic-active areas within NSCLC tumours using a pre-radiotherapy FDG-PET-CT scan: A prospective validation. Lung Canc. 75, 73–76 (2012).

N. Beerenwinkel, R. F. Schwarz, M. Gerstung, F. Markowetz, Cancer evolution: Mathematical models and computational inference. Syst. Biol. 64, e1–e25 (2014).

S. Turajlic, A. Sottoriva, T. Graham, C. Swanton, Resolving genetic heterogeneity in cancer. Nat. Rev. Genet. 20, 404–416 (2019).

G. Lorenzo , Mechanistic modelling of prostate-specific antigen dynamics shows potential for personalized prediction of radiation therapy outcome. J. R. Soc. Interface 16, 20190195 (2019).

J. Pérez-Beteta , Tumor surface regularity at MR imaging predicts survival and response to surgery in patients with glioblastoma. Radiology 288, 218–225 (2018).

J. Pérez-Beteta , Glioblastoma: Does the pre-treatment geometry matter? A postcontrast T1 MRI-based study. Eur. Radiol. 27, 1096–1104 (2017).

K. Swanson, R. Rostomily, E. Alvord Jr, A mathematical modelling tool for predicting survival of individual patients following resection of glioblastoma: A proof of principle. Br. J. Canc. 98, 113–119 (2008).

P. Wang , microRNA-21 negatively regulates Cdc25A and cell cycle progression in colon cancer cells. Canc. Res. 69, 8157–8165 (2009).

V. M. Pérez-García, G. F. Calvo, J. Belmonte-Beitia, D. Diego, L. Pérez-Romasanta, Bright solitary waves in malignant gliomas. Phys. Rev. E Stat. Nonlinear Soft Matter Phys. 84, 021921 (2011).

S. C. Massey , Image-based metric of invasiveness predicts response to adjuvant temozolomide for primary glioblastoma. PloS One 15, e0230492 (2020).

A. L. Baldock , Patient-specific metrics of invasiveness reveal significant prognostic benefit of resection in a predictable subset of gliomas. PloS One 9, e99057 (2014).

J. Mettler , Metabolic tumor volume for response prediction in advanced-stage Hodgkin lymphoma. J. Nuclear Med. 60, 207–211 (2019).

C. Bisschop, E. de Heer, A. Brouwers, G. Hospers, M. Jalving, Rational use of 18F-FDG PET/CT in patients with advanced cutaneous melanoma: A systematic review. Crit. Rev. Oncol. Hematol. 153, 103044 (2020).

J. Harrell , Evaluating the yield of medical tests. J. Am. Med. Assoc. 247, 2543–2546 (1982).

S. Hamdi, W. E. Schiesser, G. W. Griffiths, Method of lines. Scholarpedia 2, 2859 (2007).

D. Molina , Lack of robustness of textural measures obtained from 3D brain tumor MRIs impose a need for standardization. PloS One 12, e0178843 (2017).

F. Staub, A. Baethmann, J. Peters, H. Weigt, O. Kempski, Effects of lactacidosis on glial cell volume and viability. J. Cerebr. Blood Flow Metabol. 10, 866–876 (1990).

J. Van de Steene , Definition of gross tumor volume in lung cancer: Inter-observer variability. Radiother. Oncol. 62, 37–49 (2002).

M. Lagendijk , Breast and tumour volume measurements in breast cancer patients using 3-D automated breast volume scanner images. World J. Surg. 42, 2087–2093 (2018).

T. Takenaka, K. Yamazaki, N. Miura, R. Mori, S. Takeo, The prognostic impact of tumor volume in patients with clinical stage IA non–small cell lung cancer. J. Thorac. Oncol. 11, 1074–1080 (2016).