# Hybrid computational modeling highlights reverse warburg effect in breast cancer-associated fibroblasts

Sahar Aghakhani, Sacha E Silva-Saffar, Sylvain Soliman, Anna Niarakis

► To cite this version:

Sahar Aghakhani, Sacha E Silva-Saffar, Sylvain Soliman, Anna Niarakis. Hybrid computational modeling highlights reverse warburg effect in breast cancer-associated fibroblasts. Computational and Structural Biotechnology Journal, 2023, 21, pp.4196-4206. 10.1016/j.csbj.2023.08.015. hal-04192259

HAL Id: hal-04192259

https://inria.hal.science/hal-04192259v1

Submitted on 31 Aug 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

Research article

# Hybrid computational modeling highlights reverse warburg effect in breast cancer-associated fibroblasts

Sahar Aghakhani<sup>a, b</sup>, Sacha E Silva-Saffar<sup>a</sup>, Sylvain Soliman<sup>b</sup>, Anna Niarakis<sup>a, b, *</sup>

<figure><img src="image_3.png" /></figure>

Check for
updates

## ARTICLE INFO

**Keywords:**

Hybrid modeling

Metabolic reprogramming

Breast cancer

Cancer-associated fibroblasts

## ABSTRACT

Cancer-associated fibroblasts (CAFs) are amongst the key players of the tumor microenvironment (TME) and are involved in cancer initiation, progression, and resistance to therapy. They exhibit aggressive phenotypes affecting extracellular matrix remodeling, angiogenesis, immune system modulation, tumor growth, and proliferation. CAFs phenotypic changes appear to be associated with metabolic alterations, notably a reverse Warburg effect that may drive fibroblasts transformation. However, its precise molecular mechanisms and regulatory drivers are still under investigation. Deciphering the reverse Warburg effect in breast CAFs may contribute to a better understanding of the interplay between TME and tumor cells, leading to new treatment strategies. In this regard, dynamic modeling approaches able to span multiple biological layers are essential to capture the emergent properties of various biological entities when complex and intertwined pathways are involved. This work presents the first hybrid large-scale computational model for breast CAFs covering major cellular signaling, gene regulation, and metabolic processes. It was generated by combining a cell- and disease-specific asynchronous Boolean model with a generic core metabolic network leveraging both data-driven and manual curation approaches. This model reproduces the experimentally observed reverse Warburg effect in breast CAFs and further identifies Hypoxia-Inducible Factor 1 (HIF-1) as its key molecular driver. Targeting HIF-1 as part of a TME-centered therapeutic strategy may prove beneficial in the treatment of breast cancer by addressing the reverse Warburg effect. Such findings in CAFs, in light of our previously published results in rheumatoid arthritis synovial fibroblasts, point to a common HIF-1-driven metabolic reprogramming of fibroblasts in breast cancer and rheumatoid arthritis.

## 1. Introduction

Cancer is the second leading cause of mortality worldwide, accounting for nearly 10 million deaths in 2020 [1]. Recognition of its burden on patients, families, communities, and health care systems as a major public health concern boosted research on the disease in an unprecedented way. Over the past two decades, numerous prevention strategies have been developed to avoid major risk factors, including screening and early detection protocols, successfully preventing an estimated 30–50% of cancer-related deaths [1]. Additionally, a variety

of treatment options are currently available including surgical interventions, chemotherapies, radiotherapies, immunotherapies, hormone therapies, targeted drug treatments, etc. However, a significant part of drug-based therapeutic approaches still fails, mainly due to drug resistance, representing the main challenge of cancer researchers to date. Resistance can occur when cancer cells present molecular alterations driving insensitivity to drugs before treatment, called intrinsic resistance, or when cancer cells adapt to a drug while it is being administered, known as extrinsic resistance. The latter is suspected to arise from changes in the tumor microenvironment (TME), as its

**Abbreviations:** ACSN, atlas of cancer signaling network; BC, breast cancer; C-flux, carbon flux; CAFs, cancer-associated fibroblasts; CALM, curation and annotation of logical models; DEA, differential expression analysis; DEG, differentially expressed genes; ECM, extracellular matrix; FAIR, findability accessibility interoperability and reproducibility; FBA, flux balance analysis; FC, fold-change; HGNC, HUGO gene nomenclature committee; HIF-1, hypoxia-inducible factor 1; MINERVA, molecular interaction networks visualization; MIRIAM, minimum information required in the annotation of models; OXPHOS, oxidative phosphorylation; PD, process description; PMID, PubMed ID; RA, rheumatoid arthritis; RASFs, rheumatoid arthritis synovial fibroblasts; SBGN, systems biology graphical notation; SBML-qual, systems biology marked up language-qualitative; TCA, citric acid cycle; TME, tumor microenvironment.

* Corresponding author at: GenHotel – European Research Laboratory for Rheumatoid Arthritis, Univ. Evry, Univ. Paris-Saclay, Evry, France.

E-mail address: anna.niarakis@inria.fr (A. Niarakis).

<footnote><sup>a</sup> GenHotel – European Research Laboratory for Rheumatoid Arthritis, Univ. Evry, Univ. Paris-Saclay, Evry, France</footnote>

<footnote><sup>b</sup> Lifeware Group, Inria Saclay Île-de-France, Palaiseau, France</footnote>
interplay with cancer cells is often partially or totally overlooked in the
development of cancer therapy.

Yet, the importance of the TME in cancer initiation and progression
has been widely recognized for years [2-4]. While cancer develops
through genetic and epigenetic alterations, tumor growth, survival, and
metastasis are regulated through complex interactions with stromal cells
of the TME [5]. Among them, multiple studies in various cancers have
demonstrated the key role of cancer-associated fibroblasts (CAFs) [6-8].
Quiescent fibroblasts are responsible for the structural integrity of the
extracellular matrix (ECM), its nutrient supply, and contribute to wound
healing. However, CAFs present a tumor-like phenotype by engaging in
cancer cell proliferation and invasion, angiogenesis, inflammation, im-
mune system modulation, and ECM remodeling [9]. Healthy fibroblasts
appear to modify their phenotypic profile to adapt to their new envi-
ronment and progressively amplify the disease's disastrous character-
istics by shifting from passive responders to key effectors.

A variety of assumptions support the origin of CAFs activation and transformation, the most recent suggesting a role for metabolic alterations [10]. Indeed, various hallmarks of CAFs could be driven by a response to a wide range of paracrine signals (e.g., cytokines, chemokines, growth factors) emitted by cancer cells into the TME. The latter would induce CAFs to reprogram their energy production to aerobic glycolysis and produce high levels of energy-rich fuels for cancer cells to convert into biomass and proliferate faster in return. This process is known as the reverse Warburg effect [11,12].

An accurate mapping of all the knowledge related to CAFs involve-
ment in the initiation and progression of cancer at the signaling, genetic,
and metabolic levels would be beneficial to decipher the complex
mechanisms underlying their activation. As knowledge assembly is an
active field in systems biology, efforts have already been initiated with
the publication of the CAF-map within the Atlas of Cancer Signaling
Network (ACSN) [13]. The ACSN is a web-based multi-scale resource of
biological maps depicting molecular processes in cancer cells and TME.
It was developed in response to an increasing need for visualization and
simplification of complex information resulting from massive produc-
tion of biological data through high-throughput sequencing. The
CAF-map is based on the manual curation of cell-specific molecular
mechanisms from scientific literature and external databases. However,
CAFs are a heterogeneous population within the TME [14] and their
distribution across distinct cancer types is still under investigation.
Despite its extensive coverage of the role of CAFs, the CAF-map, even
updated in a cancer-specific manner, is limited regarding in-silico sim-
ulations by its static knowledge base function. Nevertheless, it can be
used as a basis to generate a cancer-specific CAF dynamic model to
extend knowledge with executable information. Indeed, contributions of
dynamic approaches to decipher metabolic reprogramming of fibro-
blasts have been highlighted recently to elucidate pathogenic mecha-
nisms and accelerate identification of novel therapeutic targets [10].

Computational modeling is essential to gain insight on the emergent
behavior of biological entities when complex and interconnected path-
ways are involved. Qualitative models based on logical relationships
between components provide an appropriate description for systems
whose mechanistic processes are unknown or lack quantitative data [15,
16]. They allow a parameter-free study of the underlying dynamical
properties of large biological pathways. Boolean models are the simplest
form of qualitative models, yet extremely powerful [17]. They assume
that all biological components are described by binary values (ON/OFF
or 1/0) and their interactions by Boolean rules ("AND", "OR", "NOT").
Therefore, this formalism is very suitable for modeling signaling or gene
regulation mechanisms carrying signal flow. In contrast, quantitative
models provide great insight into biological systems by capturing their
kinetic parameters [18]. For instance, constraint-based modeling allows
fast calculations of large networks under steady-state assumptions [19].
This formalism is very appropriate to assess metabolic properties.
Evidently, a hybrid formalism able to cover these various layers at once
would be widely beneficial to study complex diseases such as cancer.

Previous attempts in this direction include the PROM method, inte-
grating high-throughput data into constraint-based modeling [20], the
FlexFlux tool, combining flux balance analysis (FBA) and qualitative
synchronous simulations [21], or the regulatory dynamic enzyme-cost
flux balance analysis method (r-deFBA) combining metabolic dynamic
modeling and transcriptional regulation [22]. Additionally, we recently
developed a framework successfully generating a hybrid computational
model of rheumatoid arthritis synovial fibroblasts, covering cellular
signaling, gene regulation, and metabolism by coupling a cell- and
disease-specific Boolean regulatory model with a constraint-based
generic metabolic network [23]. Our approach improves on previous
attempts by extracting information from asynchronous trap-spaces
rather than synchronous stable states, providing a more comprehen-
sive analysis of the regulatory model. Additionally, the previous need for
user-defined transcription factors to target gene relationships or initial
values and qualitative states to continuous intervals equivalences, a
daunting task in large-scale regulatory models, is not required in our
framework.

Computational modeling approaches have already been suggested to
improve understanding of cancer notably with models of local invasion
and migration [24], drug transport and delivery [25], or response to
treatment [26]. The relevance of hybrid modeling approaches has also
been recognized with the development of models integrating discrete
and continuous variables to study tumor development and treatment
[27]. However, regardless of the formalism (i.e., qualitative, quantita-
tive, or hybrid), modeling efforts have not yet assessed the role of CAFs
in the TME, let alone in a cancer-specific manner.

As breast cancer (BC) has overtaken lung cancer as the world's most
frequently diagnosed cancer since 2020 [28], the study of CAFs in this
specific cancer has been initiated and some biological data are available.

To unravel the pathogenic mechanisms of cancer and accelerate the
identification of innovative therapeutic targets, we aim to leverage our
previously published coupling framework [23] and apply it to the
integrative study of breast CAFs in the TME. Built upon previous efforts
to formalize CAF-related knowledge in a detailed representation of
disease mechanisms, we provide a state-of-the-art knowledge base
further translated into a cell- and disease-specific Boolean model. By
combining the latter regulatory model with a generic constraint-based
metabolic network, this work presents the first large-scale hybrid
model for breast CAFs. It addresses gene regulation, cellular signaling,
and metabolic machineries in breast CAFs along with their involvement
in the initiation and maintenance of the hallmarks of BC. All analyses are
compiled and thoroughly annotated in Jupyter notebooks and R scripts
available on a GitLab repository (https://gitlab.com/genhotel/
breast-cafs-reverse-warburg-effect) and a Zenodo permanent archive
[29].

## 2. Methods

### 2.1. *Update of the original CAF-map*

The ACSN is an online database of multi-scale biological maps depicting molecular processes in cancer cells and TME [13]. All molecular interaction maps can be downloaded in various formats for further exploration and modifications. For this work, the CAF-map was downloaded from the ACSN 2.0 website in the standard XML format. It represents the major molecular interactions depicting the role of CAFs in the TME. However, the CAF-map displays generic CAF-specific interactions and is not specific to a certain type of cancer. It includes 681 components connected through 581 reactions and is based on the manual curation of 358 peer-reviewed articles. It is consistent with Systems Biology Graphical Notation (SBGN) standards [30] and thoroughly annotated with pertinent information and references. The CAF-map was imported on CellDesigner (version 4.4.2), a structured diagram editor for drawing regulatory and biochemical networks [31].

The original CAF-map was manually updated on multiple levels,
following precise community-driven guidelines to ensure findability, accessibility, interoperability, and reproducibility, known as FAIR principles [32]. The CAF-map's layout was edited to shift from a round cell with functional modules (e.g., "integrin signaling pathways", "motility", "growth factors production", "cytokines and chemokines production", "core signaling") to a compartmentalized map with biological compartments (i.e., extracellular space, cytosol, nucleus, mitochondria, endoplasmic reticulum, secreted compartment, and phenotypes). With the help of specialized literature and external pathway databases (e.g., KEGG PATHWAY [33], PANTHER [34]), clear attribution of species to each compartment was obtained including modification or addition of interactions when needed. This representation allows for a better illustration of signal flow from top (displaying extracellular ligands complexation with plasmic receptors), through the cytoplasm (signaling and metabolic machineries), nucleus (gene regulation), to the bottom (secreted molecules and phenotype activation). Additionally, phenotypes were deleted either to allow for a detailed description of the pathways of interest (i.e., "TCA", "glycolysis", and "ketone bodies degradation" phenotypes), when they lacked added value in CAFs involvement in TME and cancer development (i.e., "microtubule polymerization", "actin polymerization", and "septine polymerization" phenotypes), or when their only interactions were to be activated by another phenotype (i.e., "stress fibril formation" and "matrix effects" phenotypes). Semantics were updated to comply with SBGN Process Description (PD) standards [35] for graphical visualization when not fully respected. It included adorning macromolecules with their residue (e.g., "phosphorylated", "acetylated") or state of modification (i.e., "active", "inactive"). Components were named in accordance with HUGO Gene Nomenclature Committee (HGNC) IDs [36] for regulatory components and BiGG IDs [37] for metabolic species when not already. Relevant annotations regarding components and reactions (e.g., PubMed IDs (PMID), HUGO identifiers, pathways of interest) were previously stored in the note section of CellDesigner. Consistent with MIRIAM (Minimum Information Required In The Annotation of Models) standards [38], such notes were retrieved and re-inserted in the dedicated MIRIAM section of CellDesigner with the qualifier "bqbiol: isDescribedBy". This qualifier is commonly used to link a component or a reaction to the literature or data describing it (e.g., PMIDs, DOI, GEO, KEGG identifier). Finally, an update of entities and reactions was carried based on the latest published mechanistic information and corrections or additions were applied if necessary. Only specific generic CAFs relevant studies were selected to enrich the molecular interaction map according to Curation and Annotation of Logical Models (CALM) standards [39].

The state-of-the-art and standardized CAF-map V2 was made available as an online interactive map on the standalone web server MINERVA (Molecular Interaction NETwoRks VisuAlization) [40]. This platform allows for visual exploration, analysis, and management of large-scale molecular networks encoded in systems biology formats (e. g., CellDesigner, SBML, and SBGN). MINERVA also provides automated content annotation and verification along with overlaying experimental data.

## 2.2. Boolean model inference from the CAF-map V2

We used the CaSQ tool (version 1.1.0) [41] and its default parameters to generate a Boolean model from the CAF-map V2 in the standard Systems Biology Marked up Language-qualitative (SBML-qual) format [42]. Briefly, this map-to-model framework infers a Boolean model based on the topology and semantics of the molecular interaction map. Note that CaSQ allows to retain references, annotations, and layout of the molecular map in the associated model, facilitating interoperability.

The CAF-model was further made available on the Cell Collective repository of biological models [43] along with BioModels.

## 2.3. Breast cancer contextualization of the CAF-model

As outlined above, the original CAF-map from the ACSN [13] was designed as a knowledge base comprising generic CAF-specific molecular mechanisms. Despite its updates and corrections, the CAF-map V2 remains a global map, gathering information from multiple types of cancer. This is consequently true for the Boolean model inferred from the latter.

Thus, a data-driven cancer-specific contextualization of the CAF-model was conducted. To this purpose, the private and pre-processed dataset EGAD00001003808 from the European Genome-phenome Archive [44] was analyzed. It consists of RNA-Seq data from 47 CAF samples sorted from fresh BC including 28 CAFs-S1 and 19 CAFs-S4 subtypes. Both CAFs-S1 and CAFs-S4 subpopulations are myofibroblastic activated CAFs whose expressions are strictly restricted to cancer and characterized by high levels of fibroblast activating protein [45]. However, CAFs-S1 are mainly involved in tumor-like phenotypes (e.g., immunosuppression, tumor growth and proliferation, inflammation, ECM remodeling), whereas CAFs-S4 are responsible for generic core signaling, motility, and perivascular signatures. Thus, CAFs-S4 represented the control group against the aggressive CAFs-S1 group of interest. The tissue harvesting protocol, as well as the data mapping, alignment, quality control, and normalization processes are detailed in [9]. Differential expression analysis (DEA) was performed afterwards using the Limma package [46] in R (version 4.2.2) to identify differentially expressed genes (DEG) between breast CAFs-S1 and CAFs-S4. Standard significance threshold of adjusted p-value > 0.05 and absolute fold change (FC) > 1.5 were applied. Most significant DEG were depicted with an absolute FC > 1. Biological processes related to the most significant DEG were assessed by gene ontology enrichment analysis on the Gene Ontology Resource powered by PANTHER [47]. Identified DEGs were discretized to fix the CAF-model's initial conditions in accordance with BC-specific biology. Up-regulated DEG present in the regulatory model were set to an initial value of 1. Similarly, down-regulated DEG present in the regulatory model were set to an initial value of 0. Finally, the inputs of a logical model being the nodes which do not present any upstream regulators, they are suspected to exert a significant control on the model's dynamics due to the linearity of signal transduction. Thus, fixing their initial values is crucial to reproduce breast CAF-specific conditions. For inputs that would not have been fixed by the DEA, initial conditions may be extracted from BC-specific peer-reviewed literature. Overall, this effort to combine data-driven and manually curated breast CAF-specific information to initialize the regulatory CAF-model allows to confidently contextualize it to reproduce, as closely as possible, breast CAF-specific conditions.

## 2.4. Breast CAF-model regulatory behavior validation

Once contextualized, the breast CAF-model's dynamic behavior was thoroughly assessed to confirm its biological relevance through two different approaches. First, a generic validation approach was carried to confirm the involvement of individual components into breast CAF-specific signaling or gene regulatory pathways. In greater detail, a literature review was undertaken and experimental evidence for breast CAF-specific activity were retrieved from in-vitro and in-vivo studies in both humans and murine models of BC (e.g., single or multiple knock-out and knock-in experiments, genetic recombination experiments). The latter studies experimental conditions were used to initialize the breast CAF-model (e.g., a gene knock-out experiment was translated as an initial value of 0 for said gene in the CAF-model and a gene knock-in was translated as an initial value of 1). This mechanistic verification was performed in a synthetic state of the model where remaining nodes were set to 0 (i.e., absent/inactive). Specific components to test values were then switched from 0 to 1 or 1-0 depending on the experimental scenario to assess their particular contributions to the pathway. Simulations were performed on the Cell Collective interactive platform [43] in the
asynchronous updating mode, with a simulation speed of one, and a sliding window of 30. Secondly, the global behavior of the model was assessed. This comprehensive analysis, conducted at the level of the CAF-model's cellular phenotypes (e.g., "growth factor production", "immune system modulation", "hypoxia", "angiogenesis", "cytokine production") allowed to compare the model's overall behavior to biologically known breast CAFs cellular behavior. This time, the model's long-term behavior was evaluated through analysis of its trap-spaces under breast CAF-specific initial conditions identified earlier. Trap-spaces are regions of the state-space from which the system cannot escape. By definition, each trap-space contains at least one attractor of the model. However, trap-spaces may overlap with each other. Thus, minimal trap-spaces (later referred to as "trap-spaces"), or trap-spaces which do not include smaller trap-spaces, offer a good approximation of a model's attractors and faithfully capture the asymptotic behavior of Boolean models [48,49]. Trap-spaces were computed using BioLQM [50] directly using binary decision diagrams or an ASP-based solver. Their computation relies on the identification of positive and negative prime implicants for each component's function without performing simulation but rather through a symbolic approach implementing a constraint-solving method. In this context, each trap-space reflects a different subspace of breast CAFs cellular phenotypes. The latter were compared with known breast CAF overall cellular behavior extracted from the literature to confirm the model's asymptotic behavior.

## 2.5. Breast CAF-model coupling with a generic metabolic model

The breast CAF-model was further integrated with a generic core metabolic network, MitoCore [51], to cover an additional biological layer of metabolism. MitoCore is a manually curated constraint-based model of human central metabolism, including 74 metabolites, 324 metabolic reactions, and 83 transport reactions. Although initially parameterized for human cardiomyocytes, MitoCore's default parameters can be extended to a variety of biological contexts without involving cell-specific features. Indeed, cardiomyocytes can metabolize a wide range of substrates, have reactions common to many other cell types and represent the metabolism of the human heart, an organ of the utmost

importance for human health, disease, and toxicology. As a proof of concept for the generalizability of their model, MitoCore's modelers were able to successfully simulate fumarase deficiency, a disease of the nervous system, using default cardiomyocyte parameters [51]. Its realistic ATP production rates for different substrates have been proven [51], making it very valuable in the present study as opposed to larger genome-scale reconstruction of human metabolism exhibiting unrealistic results.

The coupling of the regulatory breast CAF-model with the generic MitoCore metabolic network was conducted by exploiting the asymptotic behavior of metabolic components present in both models to extract additional regulatory-based metabolic constraints. In particular, metabolic components exhibiting an asymptotic inactive behavior under regulatory influences (Fig. 1).

Briefly, cell- and disease-specific initial conditions were used to initialize the regulatory CAF-model before applying the value propagation algorithm [52] implemented in CoLoMoTo [53]. The latter allowed to evaluate the dynamical influence of initial conditions on the global behavior of the model to decrease its complexity. In this restricted space, minimal trap-spaces of the regulatory network were identified, reflecting its asymptotic behavior, and particularly the asymptotic behavior of metabolic components (i.e., metabolic enzymes and metabolites) influenced by cell- and disease-specific signaling and gene regulatory processes. The maximal value of trap-spaces relative to metabolic components were used to constrain associated MitoCore's metabolic fluxes. For every metabolic component with a projected maximal regulatory trap-space equal to 0, the flux of its associated MitoCore metabolic reactions (i.e., catalyzed reactions for enzymes or producing reactions for metabolites) were constrained to 0. This approach assumes that a metabolic enzyme with a projected maximal trap-space value equal to 0 (i.e., trap-spaces always equal to 0) expresses the inactivation of the associated enzyme by signaling or gene regulation pathways. The reactions catalyzed by the enzyme will not happen. Similarly, a metabolite-associated maximal trap-space value of 0 denotes the absence of production of the associated metabolite by signaling or gene regulation pathways. Its producing reactions did not occur. Alternatively, metabolic-associated trap-space values different

<figure><img src="image_4.png" /><figcaption>Fig. 1. General architecture of the hybrid modeling framework.</figcaption></figure>
than 0 would reflect a potential activation or production of the associated species. This qualitative information does not reflect the feasibility nor the kinetics of the associated reactions. It is not sufficient to extract additional metabolic constraints. By extracting additional constraints from metabolic components with a proven asymptotic inactive behavior under cell- and disease-specific regulatory conditions, our framework enables contextualization of MitoCore's generic metabolic model in a breast CAF-specific manner. For further details, the general architecture of the framework for hybrid modeling is provided and thoroughly detailed in [22].

## 2.6. Metabolic network analysis

A widely used method for analyzing large-scale reconstructions of metabolic networks is FBA [54]. Its main advantage lies in the need for little information regarding the kinetic parameters of biological components as it calculates metabolic fluxes by assuming steady state conditions. Two FBAs were performed using the python package CobraPy [55] to highlight a potential change in fibroblast metabolism under disease-specific conditions. First FBA was conducted without any additional constraints from the regulatory network and considered as the "control", reflecting the healthy fibroblast state. Second FBA was conducted with additional metabolic constraints extracted from the breast CAF-specific regulatory model's trap-spaces, reflecting breast CAFs condition. The FBA objective function was set to maximal cellular energy production in the form of ATP to reflect the central role of metabolism. The latter was manually defined as the sum of the three ATP producing reactions, i.e., step seven and ten of glycolysis, and complex V of oxidative phosphorylation (OXPHOS). As MitoCore lacks cellular and tissular specificity, a numerical interpretation of metabolic flux values is not possible. FBA results were rather interpreted in terms of flux distribution through the ratios of ATP production from glycolytic (i.e., fluxes of step seven and ten of glycolysis) or oxidative reactions (i.e., flux of OXPHOS complex V reaction) relative to total ATP production (i.e., objective function flux) to highlight the main metabolic pathway of ATP production. Carbon fluxes (C-flux) of uptake and secretion were also analyzed to investigate the exchanges between the cell (i.e., breast CAF) and its environment (i.e., breast TME). They respectively represent the total cellular carbon influx from a specific uptake reaction and the proportion of total cellular carbon efflux from a specific secretion reaction. Finally, a comparison of internal metabolic fluxes is performed, leading to identify metabolic alterations from a greater than 2-fold variation.

## 2.7. Identification of regulatory molecular drivers

Beyond the overall impact of breast CAFs cellular signaling and gene regulatory pathways upon their metabolic machinery, we aim to identify key regulators. As stated above, the regulatory network's inputs are suspected to significantly influence the model's dynamics due to the linearity of signal flow. Hence, successive individual knock-outs and knock-ins of regulatory inputs were conducted to identify potential key regulators. All regulatory inputs initially fixed at 1 under breast CAF-specific initial conditions were successively fixed at 0 while the other components remained at breast CAF-specific values, mimicking experimental knock-outs. Similarly, all regulatory inputs fixed at 0 under breast CAF-specific initial conditions were successively fixed at 1 while the others remained unchanged, mimicking experimental knock-ins. Each new set of initial conditions was used to initialize the CAF-model and the framework for hybrid modeling was applied as outlined above. Each FBA result illustrates the distribution of metabolic fluxes under new regulatory influences. Again, the proportion of total energy production in the form of cellular ATP production from glycolytic and oxidative pathways was reported to identify a potential alteration under new regulatory conditions.

# 3. Results

## 3.1. CAF-map V2

The CAF-map V2 is available on the MINERVA platform.<fnref n="1" /> This extensive knowledge base illustrates the major cellular signaling, gene regulation, and metabolic pathways as well as molecular mechanisms and phenotypes involving CAFs in cancer initiation and progression. It is fully compliant with SBGN PD standards for visualization [35], MIRIAM for annotation [38], and CALM for biocuration [39]. It includes 649 species (308 proteins, 95 genes, 114 RNAs, 43 simple molecules, and 89 molecular complexes), 19 phenotypes, and 544 reactions. A dozen additional peer-reviewed articles enabled corrections and extended the mechanical processes depicted in the map for a total of 368 PMIDs. This high number of references enhances confidence in presented mechanisms. Biological compartments depicting extracellular space, nucleus, mitochondrion, endoplasmic reticulum, secreted compartments, and phenotypes were further added in the map to account for biological compartmentalization and cell transport. Beyond the considerable effort of the initial CAF-map to include major cancer functional modules (e.g., "growth factors production", "cytokines", "chemokines production", "matrix regulation") and molecular pathways (e.g., "growth factors signaling pathways", "inflammatory signaling pathways", "integrin signaling pathways"), coverage in disease-specific mechanisms was expanded in the CAF-map V2. For instance, metabolic pathways initially only considered as phenotypes were fully detailed, along with addition of the calcium ion signaling pathway in the new endoplasmic reticulum compartment. Visualizations of the original CAF-map and CAF-map V2 are provided in Fig. 2.

## 3.2. Breast CAF-model

Translation of the CAF-map V2 using CaSQ generated a dynamic Boolean model of 463 nodes (including 62 inputs) and 793 interactions. The latter is publicly available on the Cell Collective platform<fnref n="2" /> and BioModels repository under the identifier MODEL2307090001. However, components and interactions depicted in the model are not cancer-specific but generic to CAFs.

The DEA conducted to contextualize the model identified 3678 DEG in CAFs-S1 vs. CAFs-S4 among the 18252 mapped genes. 1866 DEG were significantly up-regulated in CAFs-S1 vs. CAFs-S4, leading to fix an initial condition of 1 for 71 nodes of the CAF-model. Most significantly up-regulated DEG in CAFs-S1 are mainly involved in chemotaxis, ECM organization, locomotion, response to growth factors, cell adhesion, migration, motility, differentiation, and proliferation. Such observations support the classification of CAFs-S1 as the subpopulation of CAFs carrying key aggressive functions. 1812 DEG were significantly down-regulated in CAFs-S1 vs. CAFs-S4, leading to fix an initial condition of 0 for 34 nodes of the CAF-model. Most significantly down-regulated DEG in CAFs-S1 are involved in ion transport and signaling, specifically calcium ions. Finally, 41 CAF-model's inputs remained unfixed, leading to assign their values through manual curation of peer-reviewed breast CAF-specific literature. The complete list of breast CAF-specific initial conditions along with their source of attribution is outlined in Table S1 of supplementary material.

To validate the behavior of the model, generic model simulations were first compared to breast CAF-specific experimental scenarios extracted from the literature. Details of each experimental scenario, associated CAF-model initialization, and dynamic results are presented in Table S2 of the supplementary materials. Overall, regarding this

<footnote n="1">¹ https://pathwaylab.elixir-luxembourg.org/minerva/index.xhtml?id=CAF-map_V2</footnote>

<footnote n="2">² https://research.cellcollective.org/#30b2359b-38be-42cc-8885-c3ab9b8f8063</footnote>
A

<figure><img src="image_5.png" /></figure>

B

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 2.** CellDesigner visualization of the (A) original CAF-map from ACSN [12] and (B) CAF-map V2 with associated zoom-ins on glucose-related pathways. The CAF-map V2 depicts a standard-compliant formal representation and layout revision with a clear signal flow. Components are precisely attributed to biologically relevant compartments. Glucose pathways are compartmentalized and are no longer represented only as phenotypes but properly detailed. Up-to-date mechanistic information was further added.</figcaption></figure>

generic evaluation, 29 experimental scenarios were confirmed by the CAF-model out of 41 biological scenarios. Ten scenarios were not reproducible mostly due to a lack of mechanistic detail regarding specific interactions in the literature, leading to a missing or incomplete representation in the CAF-map V2 and the associated breast CAF-model. Finally, two generic scenarios were not validated as other pathways were needed concomitantly. Subsequent global model’s behavior evaluation under breast CAF-specific initial conditions, i.e., identification of its complete asymptotic behavior, depicted 128 trap-spaces (see Table S3 of the supplementary materials). Each trap-space reflects a different subspace of breast CAFs cellular phenotypes. Most components values are stable within all trap-spaces (always fixed at 0 or 1) but others vary within trap-spaces. The projection of trap-spaces restricted to the model’s ontological phenotypes (i.e., distinct cellular outcome) identifies a single trap-space (Table 1). The latter illustrates activity of aggressive phenotypes in breast CAFs, e.g., angiogenesis, fibroblast proliferation, hypoxia, matrix degradation, tumor growth and invasion, while negative regulators are inactive. These findings are consistent with the behavior of breast CAFs as described in scientific literature [56].

The asymptotic state of ontological phenotypes is the result of their regulators combined effect exerted through their Boolean rules. Thus,

<caption>**Table 1**</caption>

Projection of the breast CAF-model’s trap-spaces restricted to its ontological phenotypes.

<table><thead><tr><td>Phenotype</td><td>Trap-space</td></tr></thead><tbody><tr><td>Actomyosin contractility</td><td>1</td></tr><tr><td>Angiogenesis</td><td>1</td></tr><tr><td>Cytokine production</td><td>1</td></tr><tr><td>ECM regulation</td><td>1</td></tr><tr><td>Fibroblast markers</td><td>1</td></tr><tr><td>Fibroblast proliferation</td><td>1</td></tr><tr><td>Fibrosis</td><td>1</td></tr><tr><td>Formation of basement membrane</td><td>1</td></tr><tr><td>Growth factors production</td><td>1</td></tr><tr><td>Hypoxia</td><td>1</td></tr><tr><td>Immune system modulation</td><td>1</td></tr><tr><td>Markers of fibroblast activation</td><td>1</td></tr><tr><td>Matrix degradation</td><td>1</td></tr><tr><td>Migration into the tumor</td><td>1</td></tr><tr><td>Negative regulators of CAFs</td><td>0</td></tr><tr><td>Reactive stroma</td><td>1</td></tr><tr><td>ROS production</td><td>1</td></tr><tr><td>Tumor growth</td><td>1</td></tr><tr><td>Tumor invasion</td><td>1</td></tr></tbody></table>
the behavior of groups of biomarkers associated with each ontological
phenotype can be identified to further validate the behavior of the breast
CAF-model. For instance, the constant activation of matrix metal-
loproteinases (e.g., MMP2, MMP9), known to lead to matrix degradation
[57], within all trap-spaces, accounts for the asymptotic activation of the
associated phenotype. Their role in the coordination of ECM is also well
established [56], which supports the consistency of their trap-space with
the "ECM regulation" phenotype. Likewise, most interleukins (e.g.,
IL1A, IL1B, IL11 IL6) being active under breast CAF-specific conditions
confirms their experimentally observed function of immune system
modulation and cytokine production drivers [58]. The asymptotic
activated state of the hypoxic phenotype is explained by the sustained
activation of its main molecular driver, Hypoxia -Inducible Factor 1
(HIF-1) [59]. The latter is further involved with ROS production [59],
confirmed in the trap-spaces analysis. Likewise with the angiogenic and
reactive stroma phenotypes and their key player, VEGF [60]. Migration,
tumor growth, and invasion appear to be driven by growth factors and
Wnt signaling pathway, consistent with literature [61]. Furthermore,

activation of key fibroblastic markers (e.g., FAP, PDGF) in parallel with
inactivation of known negative regulators of breast CAFs (e.g., miR101,
miR141, CAV3, SIRT3) under initial disease-specific conditions explain
the state of both related phenotypes. Finally, fibroblast proliferation
requires EGF and TGFB signaling pathways [62]. Fibrosis and basement
membrane formation seems governed by collagens [63]. Actin, myosin,
and TNC drive actomyosin contractility [64], all confirming experi-
mental observations in breast CAFs.

3.3. *Hybrid breast CAF-model*

When projected on the breast CAF-model’s metabolic components, i.
e., metabolic enzymes and metabolites, a single trap-space is identified.
It is a stable state, which indicates that, under influences from breast
CAF-specific cellular signaling and gene regulation networks, the
asymptotic behavior of all metabolic components is fully stable. As
described in the hybrid modeling framework, additional metabolic
constraints were extracted from metabolic compounds with a proven

<figure><img src="image_9.png" /><figcaption>**Fig. 3.** Visualization of major active metabolic pathways for maximal cellular ATP production in (A) control and (B) breast CAF-specific conditions according to Flux Balance Analysis.</figcaption></figure>
inactive asymptotic behavior. This concerns 7 metabolic enzymes and 15 metabolites leading to constrain the flux of 66 unique metabolic reactions to 0 (see Table S4 and Table S5 of supplementary materials). The high number of constrained reactions is explained by enzymes commonly catalyzing a single metabolic reaction whereas metabolites were produced by numerous reactions. Additionally, certain metabolic reactions (e.g., R_PDHm, R_ICDHxm, R_GLUNm) were constrained both through their catalyzing enzymes and their produced metabolites. It reflects the consistency of breast CAF-specific cellular signaling and gene regulation machinery upon its metabolic processes. This approach allows to contextualize MitoCore, a generic core metabolic network, in the cell- and disease-specific context of breast CAFs.

### 3.4. Metabolic network analysis

Visualization of both control and breast CAF-specific FBA results are provided in Fig. 3A and B, details of metabolic fluxes in both conditions are provided in Tables S6 and S7 of supplementary materials.

Optimal fluxes for ATP production in a healthy fibroblast are OXPHOS fluxes, accounting for 96% of cellular ATP production. Main uptaken carbonated molecules are hexadecanoate (40.92% of total C-flux), glucose (29.48%), lactate (9.42%), and HCO₃ (9.34%), primary energy sources for most cells. Main secreted carbonated molecule is CO₂ (99.96%). Under breast CAF-specific regulatory conditions, optimal fluxes for ATP production are glycolytic fluxes, as they now explain 85.05% of cellular energy production in the form of ATP. Main uptaken carbonated molecules are glucose (87.72% of total C-flux) and aspartate (10.01%). Main secreted carbonated molecules are lactate (88.69%) and alanine (7.99%). Comparison of internal metabolic fluxes revealed reprogramming of major metabolic pathways for ATP production. Globally increased glycolytic fluxes along with increased lactate secretion in breast CAFs reflect a highly glycolytic metabolism. Low oxidative metabolism is demonstrated by decreased OXPHOS and TCA-associated fluxes along with decreased secretion of oxidative by-products (e.g., CO₂ and H₂O). A hypoxic environment along with anaerobic metabolism is revealed through decreased O₂ uptake and increased H⁺ secretion, associated with environmental acidity. Besides metabolic pathways for energy production, macromolecular building blocks pathways are altered in breast CAFs. Particularly, amino-acids uptake by breast CAFs is decreased (e.g., proline, glycine), along with their decreased degradation (e.g., proline), and increased secretion (e.g., serine, lysine). Fatty acids uptake is further decreased (e.g., butanoic acid) while their secretion is increased (e.g., palmitic acid). Cardiolipin synthesis is decreased in breast CAFs, coherently with former findings as they are known to regulate OXPHOS. Further pathways including folate cytosolic, reductive carboxylation, and butanoate metabolism, appear to be impacted, most likely resulting indirectly from metabolites redirection through other altered pathways or through the application of our framework's additional metabolic constraints. Finally, mitochondrial transporters are affected due to the reprogramming of mitochondrial pathways discussed above.

### 3.5. HIF-1, regulatory driver of metabolic reprogramming in breast CAFs

After computationally reproducing the experimentally observed metabolic reprogramming in breast CAFs, the challenge consists in identifying its main regulatory drivers. Successive individual knock-outs and knock-ins of breast-CAF specific regulatory inputs resulted in 147 new sets of initial conditions, and therefore 147 new FBAs (complete list of FBA results per set of initial conditions by ratio of total cellular ATP produced through glycolytic or oxidative pathways is provided in Table S8 of supplementary materials). Among the performed FBA, only one set of regulatory initial conditions reproduces a healthy metabolic profile for energy production in breast CAFs, it is condition 57 (C57). Indeed, by knocking-out HIF-1, metabolic pathways seem to recover a healthy-like distribution with cellular ATP being generated from

oxidative rather than glycolytic pathways. Other regulatory inputs do not seem to directly affect metabolic fluxes distribution in breast CAFs.

## 4. Discussion

To investigate the complex mechanisms underlying the key role of CAFs in TME for breast cancer development and maintenance, dynamic modeling is essential. In this regard, we developed the first large-scale Boolean model for breast CAFs. It covers signaling, gene regulation, and their impact upon metabolism in a cell- and disease-specific context. As a basis, the CAF-map from the ASCN resource [13] was updated to fully meet community-driven guidelines for representation [30,35], curation [39], and annotation [38] of biological data. This extensive manually constructed knowledge base provided the basis to infer a Boolean model through CaSQ tool's map-to-model framework [41]. The CAF-model was further contextualized to reproduce breast CAF-specific conditions. It was achieved by combining both data-driven and manual curation approaches to ensure high cellular specificity and confidence in the depicted interactions and initial conditions. Biological coherence of the breast CAF-model's behavior was assessed against experimental evidence extracted from the scientific literature at two distinct levels (i. e., generic pathways and global cellular behavior), each time reproducing experimental observations. Overall, generating and parameterizing a logical model from manual network construction along with data analysis and bio-curation allows to avoid issues associated with automatic and non-curated methodologies (e.g., improper reconstruction of reactions, incorrect compartmentalization or representation, lack of omics data in poorly addressed biological fields).

As displayed in our previously published workflow [23], integration of the breast CAF-specific regulatory model with MitoCore's generic reconstruction of human central metabolism [51] in a context-specific manner was achieved through extraction of additional constraints from metabolic components with a proven inactive state. We further demonstrated the adaptability of our hybrid modeling framework as it successfully coupled a new regulatory network with MitoCore. Built upon previous attempts to couple regulatory with metabolic networks [20-22], we proposed integrating a core metabolic network through interpretation of regulatory trap-spaces. It allows the assessment of the impact of its long-term regulatory behavior on the updating of metabolic sub-systems under cell- and disease-specific conditions. This approach, while capturing high-level of regulatory complexity, overcomes the need for translating discrete synchronous qualitative states into user-defined continuous intervals as well as relying solely on automated integration of high-throughput data. Overall, combining various layers of biological machinery, namely signaling, gene regulation, and metabolism, enabled a rich characterization of the role of CAFs in breast TME.

Simulations of our hybrid network constitute the first modeling attempt to address breast CAFs regulatory impact upon their metabolism. Findings depicted a hypoxic and acidic TME along with a highly glycolytic metabolism and almost null oxidative fluxes for energy production in breast CAFs. Increased production and secretion of energy-rich fuels (e.g., pyruvate, lactate, amino-acids, fatty acids) along with decreased secretion of oxidative by-products (e.g., CO₂, H₂O) were reported, fully reproducing the experimentally observed reverse Warburg effect in CAFs [65]. Fibroblasts would undergo metabolic reprogramming and turn into "metabolic slaves", generating high levels of energy-rich fuels through pentose phosphate pathway or amino acids synthesis. Indeed, beyond metabolic pathways directly associated with energy production in the form of ATP (i.e., glycolysis, OXPHOS), FBA results demonstrated various additional altered metabolic pathways in breast CAFs, including fatty acids, amino acids, and cardiolipin synthesis. In the context of the reverse Warburg effect, these nutrients would feed cancer cells in macromolecular building blocks and sustain their aggressive phenotype. Following a reverse Warburg effect, breast CAFs seem to enhance the production and secretion of amino acids to
provide nutrients for neighboring cancer cells protein synthesis pro-
cesses. The exact mechanisms by which CAFs increase amino acid syn-
thesis are not fully understood but may involve alterations in signaling
pathways, gene expression, and nutrient availability. Additionally, fatty
acids play a critical role in cancer cell growth and energy metabolism.
Breast CAFs appear to enhance fatty acid synthesis and release,
providing a source of lipids for cancer cells. This increased fatty acid
synthesis in CAFs may be driven by metabolic alterations, such as acti-
vation of the Akt/mTOR signaling pathway or upregulation of lipogenic
enzymes like fatty acid synthase (FASN) but further studies are needed.
Finally, cardiolipin is a unique phospholipid found predominantly in the
inner mitochondrial membrane where it plays a crucial role in mito-
chondrial function and energy metabolism. According to our in-silico
simulations, breast CAFs seem to exhibit alterations in cardiolipin syn-
thesis, which may impact mitochondrial function and contribute to
cancer cells metabolic reprogramming. Indeed, disrupted cardiolipin
synthesis in CAFs could lead to mitochondrial dysfunction, affecting
energy production and promoting a more favorable microenvironment
for cancer cell survival and growth. Such metabolic alterations are
suspected to be driven directly by signals from cancer cells, illustrating
an extensive crosstalk within the TME. Subsequent simulations identi-
fied HIF-1 as the regulatory activating stimuli of breast CAFs metabolic
reprogramming. HIF-1 is known as a master transcriptional factor
involved in homeostasis and cellular response to hypoxia. Indeed,
chronic hypoxia is a major hallmark of solid tumors as they quickly
outgrow their blood supply to support their continuous growth and
proliferation, leaving parts of the tumor with almost null concentration
of oxygen. Angiogenesis, metastasis, and drug resistance benefit from
this hypoxic state [66,67]. To sustain such aggressive behaviors in this
challenging environment, cancer cells need additional fueling. As a
result, this chronic hypoxic and acidic environment, generated by cancer
cells to maintain their tumorigenic phenotype, was recently suspected to
induce breast CAFs reverse Warburg effect and promote BC progression
[66]. Our findings fully support this hypothesis and further identify
HIF-1 as a main molecular regulator. In response to cancer cells hypoxic
paracrine signals, breast CAFs would activate transcription of glycolytic
genes along with glucose transporters, suppress oxygen consumption,
by-pass oxidative pathways, and induce a reverse Warburg effect. Thus,
we hypothesize that targeting the metabolic reprogramming of fibro-
blasts through HIF-1's pro-glycolytic and anti-oxidative transcriptional
activity within a therapeutic strategy including the TME could benefit
the treatment of BC. However, although we covered a variety of regu-
latory input combinations and while inputs are suspected to exert the
most significant role upon our networks, we did not systematically
assess the impact of internal components. An extensive, if not exhaustive
(for computational reasons), analysis of the combinations involving in-
ternal components could shed light on the precise mechanism of action
of HIF1 in breast CAFs.

Similar in-silico experiments were previously conducted in [23] with
a hybrid model of rheumatoid arthritis synovial fibroblasts (RASFs). We
suggested a cancer-like reverse Warburg effect in the rheumatic joint
where RASFs would reprogram their metabolism to produce more ATP,
maintain their aggressive phenotype, and feed neighboring energetically
demanding cells (e.g., chondrocytes, macrophages, dendritic cells) with
fuels and nutrients. Metabolic reprogramming of fibroblasts appears to
be a crucial element in the pathogenesis of two complex diseases as
different as BC and rheumatoid arthritis (RA). Due to the similarity of
the framework applied to contextualize the generic MitoCore network in
both diseases, a comparison of altered metabolic pathways can be per-
formed (see Table S9 of supplementary materials for RASF-specific FBA
results). Overall, metabolic pathways directly related to energy pro-
duction are altered in the same manner (i.e., increased glycolysis,
decreased OXPHOS and TCA). Oxidative by-products are accordingly
altered in both breast CAFs and RASFs (e.g., increased lactate and
secretion, decreased CO₂ and H₂O secretion). The fibroblasts environ-
ment, respectively BC TME and RA joint, appears to be similarly

modified (e.g., decreased O₂ uptake and increased H⁺ secretion). Metabolic pathways not directly involved in ATP production are additionally affected. For instance, disease-specific regulatory conditions similarly reprogram amino-acids and fatty acids pathways in RA and BC-associated fibroblasts probably due to their importance in the biosynthesis of macromolecules for aggressive cells. Indeed, secretion of building blocks is additionally increased in both situations. Certain pathways are similarly altered but raise the question of their interest under different environmental conditions. For instance, reductive carboxylation is similarly increased in breast CAFs and RASFs. Acting as a novel glutamine pathway, it supports the growth of cells depicting mitochondrial deficiencies. Warburg originally hypothesized that cancer-like cells presented a mitochondrial defect [68] but later work refuted it [69]. Such studies have not yet been conducted in both types of fibroblasts. Further experimental investigations are needed to decipher their mitochondrial status and identify a potential benefit from reductive carboxylation reactions. Same goes for cytosolic misc, needing additional experimental studies to decipher their precise beneficiaries in BC and RA. Altogether, mitochondrial and cytosolic transporters or shuttles pathways are affected, not necessarily in the same way, but all due to the reprogramming of upstream metabolic pathways producing their metabolite of interest. Finally, butanoate metabolism is not affected similarly in breast CAFs and RASFs. However, as this pathway is typically involved in processes associated with intestinal fermentation, its alterations do not seem to be significant in a cancerous or autoimmune context.

Beyond the many shared altered metabolic pathways, the key mo-
lecular regulator was also identified as HIF-1 in RASFs reverse Warburg
effect where its mechanism of action was suspected at the transcrip-
tional level. Already recognized in RA as a driver of inflammation,
angiogenesis, and cartilage destruction [70], targeting HIF-1 has not yet
been proposed within a therapeutic strategy aiming at the resolution of
metabolic reprogramming in fibroblasts. In BC, therapeutic opportu-
nities targeting HIF-1 appeared until very recently to be limited to its
metastatic or driver of tumor-proliferation activity [71]. Growing in-
terest in metabolic targeting to address pro-tumor characteristics
resulted in a few insights such as Honokiol as an inhibitor of
HIF-1-mediated glycolysis to halt BC cells growth [72]. However, its
interest in breast CAFs has not yet been investigated. Finally, studies in
other types of human fibroblasts have recognized a key role for HIF-1 (e.
g., in anti-aging and regeneration in dermal fibroblasts [73], attenuating
fibrosis and delaying vascular remodeling in systemic
sclerosis-associated fibroblasts [74]) but its targeting in the resolution of
fibroblasts metabolic reprogramming has not yet been studied.

According to our findings and given the shared metabolic alterations and regulatory molecular driver, we assume the existence of a common mechanism directing the phenotypic transformation of fibroblasts through a HIF-1-driven metabolic reprogramming in BC and RA. In both situations, regulation of CAFs and RASFs through a reverse Warburg effect enables them to adapt and survive in a new hypoxic and acidic environment, maintain their surroundings, and actively participate in the amplification of the associated disease's debilitating symptoms.

5. Conclusion

Hybrid modeling formalisms are highly relevant for the study of
complex diseases such as cancer. By covering multiple levels of inter-
twined processes and biological entities, they contribute to deciphering
the complex mechanisms at the origin of carcinogenesis and further
develop therapeutic strategies considering tumor cells surroundings. In
this regard, we present the first large-scale hybrid breast CAF-model
covering major cellular signaling, gene regulation, and metabolic pro-
cesses. It combines a cell- and disease-specific asynchronous Boolean
model with a generic core metabolic network. Generation of this hybrid
model leveraged both data-driven and manual curation approaches to
ensure high cellular specificity and confidence in the depicted
interactions and contextualization. Hybrid simulations reproduce the experimentally observed reverse Warburg effect in breast CAFs and further identify HIF-1 as its primary molecular driver. Following paradigm hypoxic signals emitted by cancer cells in the TME, breast CAFs would reprogram their energy production pathways and produce high-levels of energy-rich fuels and nutrients for neighboring tumor cells. Developing innovative treatment strategies targeting the reverse Warburg effect in CAFs through HIF-1 may represent a promising path for the treatment of breast cancer. This process appears to be further shared with rheumatoid arthritis synovial fibroblasts.

## 6. Perspectives

The presented approach, while capturing the complexity of signaling and gene regulation processes upon breast CAFs metabolism, provides a linear view of biological events. Considering the regulatory network's update from metabolic outputs may provide a more comprehensive understanding of events associated with the reverse Warburg effect in breast CAFs. Additionally, identification of regulatory molecular drivers for breast CAFs metabolic reprogramming focuses solely on regulatory inputs. The effect of internal core components is not considered although it could enable identification of HIF-1's precise mechanism of action. Finally, the regulatory inputs are tested individually but biological joint influences cannot be dismissed. Combined knock-outs and knock-ins could be considered, leading to therapeutic strategy suggestions targeting several components at once.

## Funding

This work was supported by the doctorate program of the University of Paris-Saclay, France (Sahar Aghakhani), Genopole (Sahar Aghakhani, Sacha E Silva-Saffar, and Anna Niarakis), and Inria (Sylvain Soliman and Anna Niarakis). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

## CRediT authorship contribution statement

**Sahar Aghakhani:** Methodology/Study design, Validation, Formal analysis, Investigation, Data curation, Writing – original draft, Writing – review & editing, Visualization. **Sacha E Silva-Saffar:** Validation, Formal analysis, Investigation, Data curation. **Sylvain Soliman:** Conceptualization, Methodology/Study design, Validation, Resources, Writing – review & editing, Supervision, Project administration, Funding acquisition. **Anna Niarakis:** Conceptualization, Methodology/Study design, Resources, Writing – review & editing, Supervision, Project administration, Funding acquisition.

## Declaration of Competing Interest

None Declared.

## Acknowledgements

The authors would like to thank Marek Ostaszewski (Luxembourg Centre for Systems Biomedicine) for his technical help handling Cell-Designer annotations and hosting the CAF-map V2 on MINERVA, as well as Fatima Mechta-Grigoriou and Yann Kieffer (Institut Curie) for granting access to the breast CAFs RNA-Seq data (European Genome-phenome Archive dataset EGAD00001003808).

## Appendix A. Supporting information

Supplementary data associated with this article can be found in the online version at doi:10.1016/j.csbj.2023.08.015.

## References

[1] WHO. WHO. Cancer. World Health Organization: WHO, 2022.

[2] De Wever O, Mareel M. Role of tissue stroma in cancer cell invasion. J Pathol 2003; 200:429-47. https://doi.org/10.1002/path.1398.

[3] Giraldo NA, Sanchez-Salas R, Peske JD, Vano Y, Becht E, Petitprez F, et al. The clinical role of the TME in solid cancer. Br J Cancer 2018;120:45-53. https://doi.org/10.1038/s41416-018-0327-z.

[4] Neophytou CM, Panagi M, Stylianopoulos T, Papageorgis P. The Role of Tumor Microenvironment in Cancer Metastasis: Molecular Mechanisms and Therapeutic Opportunities. Cancers 2021;13:2053. https://doi.org/10.3390/cancers13092053.

[5] Guo S, Deng C-X. Effect of stromal cells in tumor microenvironment on metastasis initiation. Int J Biol Sci 2018;14:2083-93. https://doi.org/10.7150/ijbs.25720.

[6] Asif PJ, Longobardi C, Hahne M, Medema JP. The role of cancer-associated fibroblasts in cancer invasion and metastasis. Cancers 2021;13:4720. https://doi.org/10.3390/cancers13184720.

[7] Joshi RS, Kanugula SS, Sudhir S, Pereira MP, Jain S, Aghi MK. The role of cancer-associated fibroblasts in tumor progression. Cancers 2021;13:1399. https://doi.org/10.3390/cancers13061399.

[8] Ping Q, Yan R, Cheng X, Wang W, Zhong Y, Hou Z, et al. Cancer-associated fibroblasts: overview, progress, challenges, and directions. Cancer Gene Ther 2021; 28:984-99. https://doi.org/10.1038/s41417-021-00318-4.

[9] Costa A, Kieffer Y, Scholer-Dahirel A, Pelon F, Bourachot B, Cardon M, et al. Fibroblast heterogeneity and immunosuppressive environment in human breast cancer. Cancer Cell 2018;33(463-479):e10. https://doi.org/10.1016/j.ccell.2018.01.011.

[10] Aghakhani S, Zerrouk N, Niarakis A. Metabolic reprogramming of fibroblasts as therapeutic target in rheumatoid arthritis and cancer: deciphering key mechanisms using computational systems biology approaches. Cancers 2020;13:35. https://doi.org/10.3390/cancers13010035.

[11] Warburg O, Wind F, Negelein E. The metabolism of tumors in the body. J Gen Physiol 1927;8:519-30. https://doi.org/10.1085/jgp.8.6.519.

[12] Pavlides S, Whitaker-Menezes D, Castello-Cros R, Flomenberg N, Witkiewicz AK, Frank PG, et al. The reverse Warburg effect: Aerobic glycolysis in cancer associated fibroblasts and the tumor stroma. Cell Cycle 2009;8:3984-4001. https://doi.org/10.4161/cc.8.23.10238.

[13] Kuperstein I, Bonnet E, Nguyen H-A, Cohen D, Viara E, Grieco L, et al. Atlas of Cancer Signalling Network: a systems biology resource for integrative analysis of cancer data with Google Maps. Oncogenesis 2015;4:e160. https://doi.org/10.1038/oncsis.2015.19.

[14] Wang Z. CAF heterogeneity and dynamics. Nat Cell Biol 2022;24:1686. https://doi.org/10.1038/s41556-022-01054-z.

[15] Hall B, Niarakis A. Data Integration in Logic-based Models of Biological Mechanisms. MDPI AG, 2021.

[16] Hemedan AA, Niarakis A, Schneider R, Ostaszewski M. Boolean modelling as a logic-based dynamic approach in systems medicine. Comput Struct Biotechnol J 2022;20:3161-72. https://doi.org/10.1016/j.csbj.2022.06.035.

[17] Kauffman SA. Metabolic stability and epigenesis in randomly constructed genetic nets. J Theor Biol 1969;22:437-67. https://doi.org/10.1016/0022-5193(69)90015-0.

[18] Saadatpour A, Albert R. A comparative study of qualitative and quantitative dynamic models of biological regulatory networks. EPJ Nonlinear Biomed Phys 2016:4. https://doi.org/10.1140/epjnbp/40366-016-0031-y.

[19] Yasemi M, Jolicoeur M. Modelling cell metabolism: a review on constraint-based steady-state and kinetic approaches. Processes 2021;9:322. https://doi.org/10.3390/pr9020322.

[20] Chandrasekaran S, Price ND. Probabilistic integrative modeling of genome-scale metabolic and regulatory networks in Escherichia coli and Mycobacterium tuberculosis. Proc Natl Acad Sci 2010;107:17845-50. https://doi.org/10.1073/pnas.1005139107.

[21] Marmiesse L, Peyraud R, Cottret L. FlexFlux: combining metabolic flux and regulatory network analyses. BMC Syst Biol 2015:9. https://doi.org/10.1186/s12918-015-0238-z.

[22] Liu L, Bockmayr A. Regulatory dynamic enzyme-cost flux balance analysis: A unifying framework for constraint-based modeling. J Theor Biol 2020;501:110317. https://doi.org/10.1016/j.jtbi.2020.110317.

[23] Aghakhani S, Soliman S, Niarakis A. Metabolic reprogramming in Rheumatoid Arthritis Synovial Fibroblasts: A hybrid modeling approach. PLOS Comput Biol 2022;18:e1010408. https://doi.org/10.1371/journal.pcbi.1010408.

[24] Cohen DPA, Martignetti L, Robine S, Barillot E, Zinovyev A, Calzone L. Mathematical Modelling of Molecular Pathways Enabling Tumour Cell Invasion and Migration. PLOS Comput Biol 2015;11:e1004571. https://doi.org/10.1371/journal.pcbi.1004571.

[25] Zhan W, Alamer M, Xu XY. Computational modelling of drug delivery to solid tumour: Understanding the interplay between chemotherapeutics and biological system for optimised delivery systems. Adv Drug Deliv Rev 2018;132:81-103. https://doi.org/10.1016/j.addr.2018.07.013.

[26] Béal J, Pantolini L, Noël V, Barillot E, Calzone L. Personalized Logical Models to Investigate Cancer Response to Braf Treatments in Melanomas and Colorectal Cancers. Cold Spring Harbor Laboratory, 2020.

[27] Chamseddine IM, Rejnika AK. Hybrid modeling frameworks of tumor development and treatment. WIREs Syst Biol Med 2019:12. https://doi.org/10.1002/wbm.1461.

[28] Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataran I, Jemal A, et al. Global Cancer Statistics 2020: GLOBOCAN estimates of incidence and mortality
worldwide for 36 cancers in 185 countries. CA: A Cancer J Clin 2021;71:209-49.
https://doi.org/10.3322/caac.21660.

[29] Aghakhani, S., Soliman, S., Niarakis, A. CODE - Breast C.A.F. Hybrid Model. Zenodo, 2023. (https://zenodo.org/record/7874613).

[30] Novère NL, Hucka M, Mi H, Moodie S, Schreiber F, Sorokin A, et al. The systems biology graphical notation. Nat Biotechnol 2009;27:735-41. https://doi.org/10.1038/nbt1558.

[31] Funahashi A, Morohashi M, Kitano H, Tanimura N. CellDesigner: a process diagram editor for gene-regulatory and biochemical networks. BIOSILICO, 1; 2003. p. 159-62. https://doi.org/10.1016/s1478-5382(03)02370-9.

[32] Wilkinson MD, Dumontier M, Aalbersberg IJJ, Appleton G, Axton M, Baak A, et al. The FAIR Guiding Principles for scientific data management and stewardship. Scientific. Data 2016:3. https://doi.org/10.1038/sdata.2016.18.

[33] Kanehisa M. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res 2000;28:27-30. https://doi.org/10.1093/nar/28.1.27.

[34] Mi H, Thomas P. PANTHER Pathway: An Ontology-Based Pathway Database Coupled with Data Analysis Tools. Methods in Molecular Biology. Totowa, NJ: Humana Press.; 2009. p. 123-40.

[35] Moodie S, Le Novere N, Demir E, Mi H, Schreiber F. Systems biology graphical notation: process description language level 1. Nat Preced 2011. https://doi.org/10.1038/npre.2011.3721.3.

[36] Seal RL, Braschi B, Gray K, Jones TEM, Tweedie S, Haim-Vilmovsky L, et al. Genenames.org: the HGNC resources in 2023. Nucleic Acids Res 2022;51:D1003-9. https://doi.org/10.1093/nar/gkac888.

[37] King ZA, Lu J, Dräger A, Miller P, Federowicz S, Lerman JA, et al. BiGG Models: A platform for integrating, standardizing and sharing genome-scale models. Nucleic Acids Res 2015;44:D515-22. https://doi.org/10.1093/nar/gkv1049.

[38] Novère NL, Finney A, Hucka M, Bhalla US, Campagne F, Collado-Vides J, et al. Minimum information requested in the annotation of biochemical models (MIRIAM). Nat Biotechnol 2005;23:1509-15. https://doi.org/10.1038/nbt1156.

[39] Niarakis A, Kuiper M, Oostaszewski M, Malik Sheriff RS, Casals-Casas C, Thieffry D, et al. Setting the basis of best practices and standards for curation and annotation of logical models in biology—highlights of the [BC]2 2019 CoLoMoTo/SysMod Workshop. Brief Bioinforma 2020;22:1848-59. https://doi.org/10.1093/bib/bbaa046.

[40] Gawron P, Oostaszewski M, Satagopam V, Gebel S, Mazein A, Kuzma M, et al. MINERVA—a platform for visualization and curation of molecular interaction networks. Npj Syst Biol Appl 2016:2. https://doi.org/10.1038/npjsb.2016.20.

[41] Aghamiri SS, Singh V, Naldi A, Helikar T, Soliman S, Niarakis A. Automated inference of Boolean models from molecular interaction maps using CaSQ. Bioinformatics 2020;36:4473-82. https://doi.org/10.1093/bioinformatics/bta484.

[42] Chaouiya C, Bérenguier D, Keating SM, Naldi A, van Iersel MP, Rodriguez N, et al. SBML qualitative models: a model representation format and infrastructure to foster interactions between qualitative modelling formalisms and tools. BMC Syst Biol 2013;7:135. https://doi.org/10.1186/1752-0509-7-135.

[43] Helikar T, Kowal B, McClenathan S, Bruckner M, Rowley T, Madrahimov A, et al. The Cell Collective: Toward an open and collaborative approach to systems biology. BMC Syst Biol 2012;6:96. https://doi.org/10.1186/1752-0509-6-96.

[44] Freeberg MA, Fromont LA, D'Altri T, Romero AF, Ciges JI, Jene A, et al. The European Genome-phenome Archive in 2021. Nucleic Acids Res 2021;50:D980-7. https://doi.org/10.1093/nar/gkab1059.

[45] Bonneau C, Eliès A, Kieffer Y, Bourachot B, Ladoire S, Pelon F, et al. A subset of activated fibroblasts is associated with distant relapse in early luminal breast cancer. Breast Cancer Res 2020:22. https://doi.org/10.1186/s13058-020-01311-9.

[46] Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47-e47. https://doi.org/10.1093/nar/gkv007.

[47] Mi H, Muruganujan A, Ebert D, Huang X, Thomas PD. PANTHER version 14: more genomes, a new PANTHER GO-slim and improvements in enrichment analysis tools. Nucleic Acids Res 2018;47:D419-26. https://doi.org/10.1093/nar/gky1038.

[48] Zañudo JGT, Albert R. An effective network reduction approach to find the dynamical repertoire of discrete dynamic networks. Chaos: An Interdisciplinary. J Nonlinear Sci 2013;23:025111. https://doi.org/10.1063/1.4809777.

[49] Klarner R, Bockmayr A, Siebert H. Computing Symbolic Steady States of Boolean Networks. Lecture Notes in Computer Science. Cham: Springer International Publishing.; 2014. p. 561-70.

[50] Naldi A. BioLQM: A Java Toolkit for the Manipulation and Conversion of Logical Qualitative Models of Biological Networks. Front Physiol 2018:9. https://doi.org/10.3389/fphys.2018.01605.

[51] Smith AC, Eyassu F, Mazat J-P, Robinson AJ. MitoCore: a curated constraint-based model for simulating human central metabolism. BMC Syst Biol 2017:11. https://doi.org/10.1186/s12918-017-0500-7.

[52] Hernandez C, Thomas-Chollier M, Naldi A, Thieffry D. Computational Verification of Large Logical Models—Application to the Prediction of T Cell Response to Checkpoint Inhibitors. Front Physiol 2020:11. https://doi.org/10.3389/fphys.2020.558606.

[53] Naldi A, Hernandez C, Levy N, Stoll G, Monteiro PT, Chaouiy C, et al. The CoLoMoTo interactive notebook: accessible and reproducible computational analyses for qualitative biological networks. Front Physiol 2018:9. https://doi.org/10.3389/fphys.2018.00680.

[54] Orth JD, Thiele I, Palsson BØ. What is flux balance analysis? Nat Biotechnol 2010; 24:245-8. https://doi.org/10.1038/nbt.1614.

[55] Ebrahim A, Lerman JA, Palsson BO, Hyduke DR. COBRAPy: Constraints-based reconstruction and analysis for python. BMC Syst Biol 2013:7. https://doi.org/10.1186/1752-0509-7-74.

[56] Fernández-Nogueira P, Fuster G, Gutierrez-Uzquiza A, Gascón P, Carbó N, Bragado P. Cancer-associated fibroblasts in breast cancer treatment response and metastasis. Cancers 2021;13:3146. https://doi.org/10.3390/cancers13133146.

[57] Singer CF, Kronsteiner N, Marton E, Kubista M, Cullen KJ, Hirtenlehner K, et al. MMP-2 and MMP-9 Expression in Breast Cancer-Derived Human Fibroblasts is Differentially Regulated by Stromal-Epithelial Interactions. Breast Cancer Res Treat 2002;72:69-77. https://doi.org/10.1023/a:1014918512569.

[58] Fasoulakis Z, Kolios G, Papamanolis V, Kontomanolis EN. Interleukins associated with breast cancer. Cureus 2018. https://doi.org/10.7759/cureus.3549.

[59] Chiavarina B, Whitaker-Menezes D, Migneco G, Martinez-Outschoorn UE, Pavlides S, Howell A, et al. HIF1-alpha functions as a tumor promoter in cancer-associated fibroblasts, and as a tumor suppressor in breast cancer cells. Cell Cycle 2010;9:3534-51. https://doi.org/10.4161/cc.9.17.12908.

[60] Fukumura D, Xavier R, Sugiura T, Chen Y, Park E-C, Lu N, et al. Tumor Induction of VEGF Promoter Activity in Stromal Cells. Cell 1998;94:715-25. https://doi.org/10.1016/s0092-8674(00)81731-6.

[61] Ruan Y, Ogana H, Gang E, Kim HN, Kim Y-M. Wnt signaling in the tumor microenvironment. Advances in Experimental Medicine and Biology. Cham: Springer International Publishing.; 2020. p. 107-21.

[62] Hofland LJ, van der Burg B, van Eijck CHJ, Sprij DM, van Koetsveld PM, Lamberts SWJ. Role of tumor-derived fibroblasts in the growth of primary cultures of human breast-cancer cells: Effects of epidermal growth factor and the somatostatin analogue octreotide. Int J Cancer 1995;60:93-9. https://doi.org/10.1002/jc.2910600114.

[63] Willumsen N, Jensen C, Green G, Nissen NI, Neely J, Nelson DM, et al. Fibrotic activity quantified in serum by measurements of type III collagen pro-peptides can be used for prognosis across different solid tumor types. Cell Mol Life Sci 2022:79. https://doi.org/10.1007/s00018-022-04226-0.

[64] Katoh D, Kozuka Y, Noro A, Ogawa T, Imanaka-Yoshida K, Yoshida T. Tenascin-C Induces Phenotypic Changes in Fibroblasts to Myofibroblasts with High Contractility through the Integrin αvβ1/Transforming Growth Factor β/SMAD Signaling Axis in Human Breast Cancer. Am J Pathol 2020;190:2123-35. https://doi.org/10.1016/j.aajpath.2020.06.008.

[65] Liang L, Li W, Li X, Jin X, Liao Q, Li Y, et al. 'Reverse Warburg effect' of cancer-associated fibroblasts (Review). Int J Oncol 2022:60. https://doi.org/10.3892/ijo.2022.5357.

[66] Muz B, de la Puente P, Azab F, Azab AK. The role of hypoxia in cancer progression, angiogenesis, metastasis, and resistance to therapy. Hypoxia 2015:83. https://doi.org/10.2147/hp.s93413.

[67] Kim I, Choi S, Yoo S, Lee M, Kim I-S. Cancer-associated fibroblasts in the hypoxic tumor microenvironment. Cancers 2022;14:3321. https://doi.org/10.3390/cancers14143321.

[68] Warburg O. On the origin of cancer cells. Science 1956;123:309-14. https://doi.org/10.1126/science.123.3191.309.

[69] Weinhouse S. The Warburg hypothesis fifty years later. Zeitschrift Für Krebsforschung Und Klinische. Onkologie 1976:87. https://doi.org/10.1007/bf00284370.

[70] Hua S, Dias TH. Hypoxia-Inducible Factor (HIF) as a Target for Novel Therapies in Rheumatoid Arthritis. Front Pharmacol 2016:7. https://doi.org/10.3389/fphar.2016.00184.

[71] Luo S, Jiang Y, Anfu Zheng, Zhao Y, Wu X, Li M, et al. Targeting hypoxia-inducible factors for breast cancer therapy: A narrative review. Front Pharmacol 2022:13. https://doi.org/10.3389/fphar.2022.1064661.

[72] Yi X, Qi M, Huang M, Zhou S, Xiong J. Honokiol Inhibits HIF-1α-Mediated Glycolysis to Halt Breast Cancer Growth. Front Pharmacol 2022:13. https://doi.org/10.3389/fphar.2022.796763.

[73] Pagani A, Kirsch BM, Hopfner U, Aitzetmueller MM, Brett EA, Thor D, et al. Deferiprone Stimulates Aged Dermal Fibroblasts via HIF-1α Modulation. Aesthetic Surg J 2020;41:514-24. https://doi.org/10.1093/asj/sja142.

[74] Xiong A, Liu Y. Targeting Hypoxia Inducible Factors-1α As a Novel Therapy in Fibrosis. Front Pharmacol 2017:8. https://doi.org/10.3389/fphar.2017.00326.