Metastatic cells exploit their stoichiometric niche in the network of cancer ecosystems

The network structure of metastases

We study large-scale patterns in metastases across 28 organs based on 9303 published metastatic records (Fig. 1, A and B, and the Supplementary Materials), which allow us to quantify organotropism and reveal variation in the degree of generality and specificity in both primary tumors and metastatic sites. The topological measures of nestedness of the metastatic network quantify the variability in organotropism [temperature, 6.46; nestedness based on the overlap and decreasing fill (NODF2), 47.4; and weighted-interaction nestedness estimator (WINE), 0.84; Fig. 1C], which is significantly greater than the null models (fig. S1). Analogous to habitat patches (33), nestedness can arise as a macroscopic property (34) when there is a gradient in suitability leading to differences in colonization or extinction probabilities (35). In addition, a truncated power-law degree distribution, especially in the source organ degree distribution (Fig. 1D and table S1), indicates a large heterogeneity in organ connectivity, where some organs are highly connected to others through metastases, as is characteristic of ecological networks (36), although it should be interpreted carefully, considering the small size of the metastatic network. The network also shows significant modularity, evidencing densely connected groups of organs (z-score > 55, Fig. 1E and fig. S1). These modules were observed when primary tumors occur in organs with high cellular turnover (e.g., skin, stomach, small intestine, and lung), forming clusters with acceptor organs having lower cellular turnover but among the most metabolically active, such as the brain, thyroids, or adrenal glands (37). This pattern indicates potential drivers associated with cell proliferation, competition, and resource exploitation (29).

Moreover, to study the role of the vascular system in defining spatial relationships between organs in the sense of Ewing’s structural hypothesis, we assembled a binary vascular matrix representing spatial relationships between organs as a function of shared main blood vessels. We tested the statistical association between the vascular matrix and the metastatic incidence matrix (Kendall’s rank correlation). The correlation between the matrices of metastatic incidence and vascularity for 19 source/acceptor organs with defined local vascular topology is significant but small in magnitude (Kendall’s rank correlation, tau = 0.126, z = 2.7, P = 0.007). These network patterns provide evidence of a connection between physiological and metastatic attributes.

Macronutrient heterogeneity, primary tumors, and successful metastases

Successful colonization emerges from the compatibility of the metastatic cell’s niche and the acceptor’s microenvironmental conditions (38). To test the P-driven metastatic colonization hypothesis, we reason that a higher likelihood of colonization would be possible in organs with resources to support higher energy production based on the metastatic cell’s requirements needed to fulfil the demands of a fast proliferative strategy (17, 18, 39, 40). The P content is a measure of phosphate-rich molecules, such as high-energy molecules (ATP), membrane phospholipids, signaling pathways, and RNA/DNA phosphodiester backbones (21, 22), and it could give clues about an organ’s conditions to support a demanding metabolism. Assuming metabolic constraints, phosphorus content is predicted to be correlated with the incidence of primary tumors and with metastatic potential. Specifically, in organs with low habitat quality (low tissue P content), cancer cells forming a primary tumor have stronger advantages due to their metabolic shift; hence, a negative association between age-adjusted cancer incidence rates and organ P content is expected, in a general sense, regardless of the organ being source or acceptor of metastases. In addition, a negative relationship between tissue P content and its number of links as a source (source degree, k_S), a positive association with acceptor degree (k_A), and the P difference (source-acceptor, ΔP) negatively correlated with the metastatic relative occurrences (number of occurrences for a source-acceptor pair over the total number of records). To assess the habitat quality of organs, we use literature estimates of phosphorus (P) content per gram of organ/tissues (41) for the source and acceptor organs. This information was available for 19 organs, involving 4934 metastatic records in the present database. The age-adjusted cancer incidence rates in primary tumors [from (42)] have a significant negative correlation with the organ’s P content (N = 11, Pearson’s ρ = −0.78, 95% confidence interval: −0.941 to –0.346, P = 0.0044; fig. S2). Regarding the metastatic cases, we found that primary tumors in low-P organs are more invasive in terms of the number of acceptor organs they reach, and this is manifested in a significant negative effect of P on k_S [negative binomial regression, effect estimate (1,16) = −593.99, z = −2.72, P = 0.007], while no significant effect of P content on k_A was detected [negative binomial regression, effect estimate (1,16) = 95.03, z = 0.74, P = 0.46]. The larger the phosphorus content in acceptor tissues in comparison to the source (negative ΔP), the higher the relative occurrence of metastasis [zero-inflated regression, R² = 0.96, occurrences (Poisson): effect estimate of ΔP = −0.42, SE = 0.013, z = −32.13, P < 0.001; absences (binomial): effect estimate of ΔP = 0.61, SE = 0.103, z = 5.96, P = 2.47 × 10⁻⁹] (Fig. 2). Together, these results indicate that P limits both primary tumor incidence and metastatic spread. In particular, supporting the P hypothesis, a low P in the source relative to an acceptor organ may be necessary to generate metastases.

Phosphorus diversity across organs and structural properties of the metastatic network

Last, we sought to test whether the observed phosphorus variability across organs (Fig. 2) could contribute to explaining the macroscopic patterns observed in the metastasis network (Fig. 1). We simulated six scenarios of metastases in which the probability of metastasis depends on a subset of the observed organ’s attributes: stem cell divisions of the source organ, blood flow through the organ, topological-vascular constraints producing propagules filtering in the liver in the case of gastrointestinal (GI) tumors, filter flow analysis sensu (13, 14), and the source-acceptor difference in P (ΔP) (Fig. 3A). We include stem cell divisions (cumulative stem cell divisions per lifetime) as a driver of the probability of having a primary tumor in a particular organ (43) because it is a necessary condition for a metastatic occurrence. All the other variables are linked to the acceptor organ. Consistently, including ΔP in the simulation produces a significantly higher correlation between observed and simulated metastatic relative occurrences (Fig. 3B). In addition, the simulation scenarios with ΔP as a codeterminant of metastasis generate nested and modular metastatic networks similar to the observed values (Fig. 3C, figs. S3 and S4, and table S4). This suggests that the observed phosphorus difference as a proxy of habitat quality across organs contributes to generating observed macroscopic patterns of the network, along with the combined effect of stem cell divisions or propagule pressure mediated by blood flow.

Ecological constraints and the stoichiometric niche of metastatic cancer cells

Analogies from ecology (44, 45) can bring new perspectives to understanding metastasis. Hence, the seed and soil hypothesis has guided our conceptual understanding of metastasis (5, 46) but at the cost of focusing on single source-acceptor links and overemphasizing the attributes of the acceptor “soil” organ relative to the source organ. By integrating the structural, seed-and-soil, and growth rate hypotheses powered by a network perspective, we provide a more general statistical assessment emphasizing the vital role played by tissue stoichiometry (20) in priming metastatic propagules and defining suitable conditions for metastasis. These microscopic molecular and cellular phenomena produce observed macroscopic patterns of metastasis.

Analogous to findings on the influence of nutrients in species invasions (47), our results indicate that phosphorus content is decisive for metastatic progression. The difference in phosphorus content between source and acceptor organs may act as an environmental filter, with metastatic propagules being more successful in habitats/organs with higher P content (16). Tumor cells from low-P organs will have better chances of survival and proliferation in acceptor tissues due to their low basal P requirements, having substantial competitive advantages in P-rich organs due to their faster phosphorus metabolism (29, 48, 49). Cancer onset may emerge by chance (9), but the stoichiometric landscape could configure tumor progression from the primary site and potential metastatic cascades and the habitat suitability present in the acceptor tissue such as nutrients or amino acids (30), where metabolic flexibility and plasticity are associated with the premetastatic niche (50). However, the link between the metabolic shift of cancer cells and P requirements may be more complex than we think. Recent evidence from murine models suggests that although cancer cells in primary orthotopic breast tumors showed the expected reduction in oxidative metabolism, this was increased in lung metastases (51). Further experimental studies exploring the energetic landscape across organs, under different health and neoplastic conditions, are required to further advance our understanding of the mechanisms underlying ATP production and its relation to P requirements in cancer cells.

A legacy effect, however, seems to be in action such that the priming of the stoichiometric niche of the pioneer migrant diaspora acquired in the source tissue affects its future success of invading a secondary organ (16), such as has been suggested for lung cancer (17) and shown for metastatic breast cancer (39, 40). This is analogous to a maternal effect, as known in evolutionary ecology, whereby there is a causal influence of the maternal genotype or phenotype (i.e., the primary organ) on the offspring phenotype (i.e., the primary tumor and its metastatic propagules). In the case of cancer, since the original healthy cell already grows under a given phosphorus availability, that basal nutrient requirement is passed on, as a constraint, to the tumor cells and their metastases.

Phosphorus dynamics have been associated with overall cancer risk in the Swedish population (52), in lung cancer (53, 54), and with an increased risk of lethal and high-grade prostate cancer (55, 56). Consequently, phosphorus could be a marker for tumor progression (22, 57, 58), and the expression of phosphorus-associated genes, e.g., ABCC11, could serve as biomarkers of metastasis. Phosphates are also involved in protein phosphorylation (59) and act as a mitogenic factor inducing proliferation and activating cellular growth (60). Inorganic phosphate may promote the establishment of cancer cells via cell-mediated angiogenesis, dependent on the forkhead box protein C2, osteopontin, and vascular endothelial growth factor α (61). Moreover, the vitamin D–dependent regulation of cell morphology could be associated with phosphorus metabolism (62, 63) and coupled with the effect of chemokines, such as CXCL8, regulating the epithelial-mesenchymal transition (EMT) and immune infiltration (64). Although the mechanistic pathway between phosphorus and EMT is still unknown, the ecological stoichiometry of P may contribute to an integrative framework for unveiling epigenetic mechanisms underpinning EMT (65, 66) and the discovery of novel biomarkers and therapeutic targets.

Nutrient availability may affect not only growth dynamics at metastatic sites but also the migration process (31). Cell migration and invasion are the energy-intensive processes, requiring changes in cytoskeleton organization, synthesis of metalloproteinases, and microtubular remodeling during migration across the extracellular matrix from the primary tumor microenvironment and in the acceptor organ (67–70). These increased energetic costs are reflected in a positive correlation between migration potential and intracellular ATP:adenosine 5′-diphosphate ratio (68) and the selection of migration paths that require less energy (71). Thus, it is expected that migrating cells will follow resource gradients (70). In this context, our data suggest that cell migration could, in particular, be affected by gradients in key resources for ATP synthesis and growth, such as phosphorus. Nevertheless, our data are limited, and further tests of our predictions together with a comprehensive understanding and prediction of metastatic routes will come once we are able to produce systematic data on the human body ecosystem, from cells, as well as organs, to epidemiological and macroecological patterns.

Ecological insights have proven helpful in the study of cancer dynamics (72, 73). A network approach coupled with an ecological framework provides new insights into metastatic phenomena (74) by studying macroscopic attributes (nestedness, degree distribution, and modularity), which are common in other complex systems (75–77). Harnessing niche and stoichiometric ecological theories (19, 29, 32) show that the emerging macroscopic patterns of metastatic incidence result from local ecologies orchestrating the habitat suitability for metastatic spread. Thus, ecological analysis of individuals as ecosystems could be essential to understanding tumor biology, evolution, and clinical progression, providing a framework for potential therapeutic targets.

Metastatic records

We studied the metastatic network by constructing a bipartite network between organs as nodes where a primary neoplasm in a source organ generated metastasis in an acceptor organ. In operational terms, the metastatic link is quantified as the number of metastases from a source to an acceptor organ. We included 9303 metastatic occurrences from five different sources (78–82). These medical records are from the United States, Switzerland, Germany, and Slovenia, taken between 1885 and 2009 (view references for years and occurrences). On the basis of Disibio and French (78), 33 anatomical zones (hereafter, “organs”) were initially recorded, maintaining most of the organs present in (78). Given the constraints of the data, we were unable to filter by sex or other demographic categories, but we decided to exclude organs commonly associated with sex (breast, prostate, penis, testicles, uterus, vulva, vagina, and ovaries). In addition, we grouped organs as follows to reconcile contrasting terminology between studies: large intestine (colon, rectum, and anus), small intestine (appendix, duodenum, and small intestines), and head and neck (neck, lip, pharynx, salivary glands, and tonsil). After this data filtering, a total of 28 organs and 9303 occurrences were included in the analysis on the basis of autopsies and tomographies in the case of muscular cancers.

Network analysis

We define the metastatic process as a graph G = (S, A, E), where S and A denote the set of source and acceptor organs, respectively, and E identifies the links connecting them, which, in this case, represents occurrences. We studied the weighted network W = G = (S, A, E) of source-acceptor organ interactions S × A with S = A = {s_i, a_j, …, N}, with i, j ∈ N and N corresponds to the 28 organs. The network is studied as a matrix where each entry represents the metastatic process quantified by n_i,j corresponding to the number of occurrences of the metastatic pair source-acceptor, i.e., the number of reported secondary growth metastases in an acceptor organ i that originated from a primary tumor in a source organ j. For network calculations, each cell value n_i,j was transformed to a metastatic relative occurrence value, calculated as w_i,j = n_i,j/T, with T being the total number of cases (9303 cases).

We quantify three main macroscopic attributes of the metastatic incidence network: nestedness (76, 83), the distribution of metastatic links (degree distribution) (84), and modularity (85). All indices to estimate these attributes are available in the R package bipartite (86). We quantify nestedness using the nestedness temperature, NODF2, and the WINE (87). For the degree distribution, we used a modified version of the bipartite::degreedistr() function to find the best fit among exponential, power-law, and truncated power-law distributions based on R², Akaike’s information criterion, and the Bayesian information criterion. For modularity, we compute modularity and module detection based on the Beckett algorithm as suggested for the case of weighted bipartite networks (85) method available in bipartite.

To discard spurious macroscopic attributes resulting from sampling effects, we compared the observed estimates against two types of null models with 1000 replicates each [package vegan (88)]. The first and second null models, in addition to maintaining the total number of occurrences, keep the marginal sum of source (method “r0_ind”) and acceptor organ (method “c0_ind”) occurrences, respectively. In fig. S1, we present the raw values of nestedness and modularity and the corresponding z-score for each model. To compare the observed network metrics (nestedness, degree distribution, and modularity) with null models, we performed a z test between the observed raw estimators against the estimators calculated from the null networks.

Vascular, physiological, and stoichiometric data

As an approach to evaluating Ewing’s structural hypothesis, we assembled a binary vascular matrix representing spatial relationship between organs where an occurrence denotes the existence of a shared main blood vessel, hence both organs coexisting in a vascular territory. For example, the adrenal glands share blood supply with the esophagus and the diaphragm through the inferior phrenic arteries. We tested the statistical association between the resulting vascular matrix and the metastatic incidence matrix using Kendall’s rank correlation.

Physiological data were only available for 21 of the 28 organs (table S2). We evaluate the hypothesis that cancer metastatic incidence of primary tumors in source organs and the source organ degree are associated with the cumulative number of stem cell divisions in a lifetime. In the case of acceptor organs, we evaluated whether the blood flux (in milliliters per minute per gram) could have a significant effect on the success of neoplastic propagules invading acceptor tissue, i.e., explaining acceptor metastatic incidence and the acceptor organ degree.

To assess the seed and soil hypothesis and unveil the effect of stoichiometry (organ nutrient composition) on metastatic incidence, we obtained literature-based estimates of macronutrient content (in grams) for a subset of organs adjusted by organ size, hereafter organ’s mass-adjusted P content, abbreviated as P.g (grams of P over grams organ’s mass; n = 19; table S3). Hence, we built an n × n subnetwork with the source-acceptor difference in the nutrients in each entry (ΔP = P.g^source – P.g^acceptor). To quantitatively assess the pattern, we performed binomial regressions between the source or acceptor organ’s degree as a function of the organ’s mass-adjusted P content. In addition, we evaluated the effect of ΔP on occurrences and absences of metastases using a zero-inflated regression model with a Poisson distribution (log link) for occurrences and a binomial distribution (logit link) for absences, package pscl (89). For visualization, we computed nonlinear quantile regression for the first, fifth, and ninth deciles.

Incidence of primary tumors data

To evaluate the association between phosphorus content and the incidence of primary tumors, we matched the organ’s P.g records with age-adjusted cancer incidence for N = 11 organs obtained from the Surveillance, Epidemiology, and End Results Program for the years 1999 to 2003, including males and females (42). After logarithmic transformation on both variables, we used Pearson’s correlation to test the strength of the association.

Simulated metastases and emerging properties

To evaluate whether the observed physiological and stoichiometric relationships could reproduce the observed network structure in terms of nestedness and modularity, we simulate four different metastatic scenarios (see below), with 1000 replicates for each scenario. The number of organs in each scenario corresponds to the data availability (table S3). For each simulated metastatic network, we run 10,000 random metastatic events. A metastatic link is defined by the probability of metastasis P_j,i from a source organ i ∈ S to an acceptor organ j ∈ A depending on the following:

Scenario I: The lifetime cumulative number of stem cell divisions in the source organ i (cumlscd_i)

Scenario II: Scenario I multiplied by the normalized difference in phosphorus content between the corresponding source and acceptor pair

Scenario III: The blood flow (flow) to an acceptor organ j

Scenario IV: Scenario III multiplied by the difference in normalized phosphorus content between the corresponding source and acceptor pair

To account for blood filtering in the liver coming from GI organs, two other scenarios consider that when by chance a GI primary tumor and nonliver organs are chosen during the simulation, the probability of metastasis Pr_j,i is halved [filter flow sensu (13)]. Then, scenarios V and VI are the same as scenarios III and IV, respectively, except for primary gastric tumors which P_j,i is halved because of a filtering of blood in the liver.

Since P can take negative values, we transformed it by adding the absolute value of min(ΔP). The ratio $\frac{\text{Δ}{P}_{j,i}}{{\displaystyle \sum _{i=1}^A}\text{Δ}{P}_{j,i}}$ was then calculated across acceptor organs for each source. We did not simulate a network scenario with all three factors (lscdcum, flow, and ΔP) included because of only five organs with records for all these parameters.

For each simulated network, metastatic relative occurrences were calculated as the number of successful metastatic links for a source-acceptor pair over the total number of links. In addition, topological descriptors of nestedness (temperature, NODF, and WINE) and modularity were calculated for each simulated network. We compared these metrics from the simulated metastases with the observed topological descriptors recalculated on the basis of the organs considered in the simulation. We computed Pearson’s correlation between the two matrices to compare the observed and simulated metastatic incidence matrix for each case. Last, we compared the observed nestedness and modularity indices with values calculated for the simulated networks between scenarios for which we obtained mean/median values and the range.

This PDF file includes: