Re: 폐포의 대식세포 M1, M2분화 과정... 2024 nature!!

[문형철]

1. nature  
2. scientific reports  
3. articles  
4. article

Dissecting metabolic landscape of alveolar macrophage

Download PDF

• Article
• Open access
• Published: 05 December 2024

Dissecting metabolic landscape of alveolar macrophage

• Sunayana Malla, 
• Karuna Anna Sajeevan, 
• Bibek Acharya, 
• Ratul Chowdhury & 
• Rajib Saha 

Scientific Reports volume 14, Article number: 30383 (2024) Cite this article

• 2295 Accesses

• 1 Citations

• Metricsdetails

Abstract

The highly plastic nature of Alveolar Macrophage (AM) plays a crucial role in the defense against inhaled particulates and pathogens in the lungs. Depending on the signal, AM acquires either the classically activated M1 phenotype or the alternatively activated M2 phenotype. In this study, we investigate the metabolic shift in the activated phases of AM (M1 and M2 phases) by reconstructing context specific Genome-Scale Metabolic (GSM) models. Metabolic pathways such as pyruvate metabolism, arachidonic acid metabolism, chondroitin/heparan sulfate biosynthesis, and heparan sulfate degradation are found to be important driving forces in the development of the M1/M2 phenotypes. Additionally, we formulated a bilevel optimization framework named MetaShiftOptimizer to identify minimal modifications that shift one activated state (M1/M2) to the other. The identified reactions involve metabolites such as glycogenin, L-carnitine, 5-hydroperoxy eicosatetraenoic acid, and leukotriene B4, which show potential to be further investigated as significant factors for developing efficient therapy targets for severe respiratory disorders in the future. Overall, our study contributes to the understanding of the metabolic capabilities of the M1 and M2 phenotype of AM and identifies pathways and reactions that can be potential targets for polarization shift and also be used as therapeutic strategies against respiratory diseases.

초록

폐 알베올라 대식세포(AM)의 높은 가소성은 

흡입된 입자 및 병원체에 대한 방어 메커니즘에서 핵심적인 역할을 합니다. 

신호에 따라 AM은 전통적으로 활성화된 

M1 형질 또는 대안적으로 활성화된 M2 형질을 획득합니다. 

본 연구에서는 

AM의 활성화 단계(M1 및 M2 단계)에서 발생하는 대사 변화를 재구성된 

맥락 특이적 게놈 규모 대사(GSM) 모델을 통해 조사했습니다. 

피루vate 대사, 아라키돈산 대사, 콘드로이틴/헤파란 설페이트 생합성, 헤파란 설페이트 분해와 같은 대사 경로는 

M1/M2 형질 발달의 중요한 촉진 요인으로 확인되었습니다. 

또한, 한 활성화 상태(M1/M2)에서 

다른 상태로 전환시키는 최소 수정 사항을 식별하기 위해 

이중 최적화 프레임워크인 MetaShiftOptimizer를 제안했습니다. 

식별된 반응에는 

글리코겐인, L-카르니틴, 5-하이드로페록시 이코사테트라엔산, 류코트리엔 B4와 같은 대사물이 포함되며, 

이는 미래에 심각한 호흡기 질환의 효율적인 치료 표적 개발을 위한 중요한 요인으로 

추가 연구가 필요할 수 있습니다. 

전체적으로, 본 연구는 

AM의 M1 및 M2 형질의 대사적 능력을 이해하는 데 기여하며, 

극성 전환의 잠재적 표적이 될 수 있는 경로와 반응을 식별했으며, 

이는 호흡기 질환에 대한 치료 전략으로 활용될 수 있습니다.

Similar content being viewed by others

Macrophage metabolic reprogramming during chronic lung disease

Article 12 November 2020

Lipid metabolism-related genes are involved in the formation of macrophage extracellular traps in allergic airway inflammation

Article 09 January 2025

Multi-omic meta-analysis identifies functional signatures of airway microbiome in chronic obstructive pulmonary disease

Article 27 July 2020

Introduction

Alveolar macrophages (AMs) serve as the first line of defense against respiratory pathogens and exhibit a high degree of plasticity1. Depending on their interactions with pathogens, AMs can differentiate into distinct subsets, with the two main phenotypes being classically activated (M1) and alternatively activated (M2) macrophages2,3. M1 macrophages, typically induced by microbial components such as lipopolysaccharide (LPS) and Th1 pro-inflammatory cytokines, are essential for the bacterial clearance and the recruitment of additional immune cells4,5,6,7. In contrast, M2 macrophages are promoted by interleukin-4 (IL-4), interleukin −13(IL-13), and immune complexes (ICs) and are associated with anti-inflammatory functions, including the resolution of inflammation and tissue repair8. However, depending on the different stimuli, M2 macrophages can be further divided into four different subtypes: M2a, M2b, M2c and M2d9. These subtypes are highly diverse, and their functions range from recruiting of other immune cells to promoting tissue repair9. Given these insights, it is crucial to further explore the role of these metabolic pathways within the specific context of the lung microenvironment to better understand their functional significance in AM polarization10.

During infection, alveolar macrophage polarization can shift dynamically based on the type of pathogen and stage of immune response4. Alveolar macrophages polarize towards the M1 phenotype in response to microbial stimuli, such as bacterial lipopolysaccharides (LPS) or viral components. On the other hand, the macrophage polarization may shift towards the M2 phenotype5when the infection is controlled, and the inflammatory response needs to be regulated. The shift in phenotypes of AM plays an important role in regulating the body’s immune response and metabolism11,12. Yet, the regulatory mechanisms governing the polarization are not completely understood. Due to the challenges associated with experimental studies, systems biology approaches can be utilized for reconstructing context specific (i.e., M1 or M2 phase) Genome-Scale Metabolic (GSM) models using available omics data. Systems biology has been proven to be very useful for pragmatic modelling and theoretical exploration of complex biological systems13. By the integration of high-throughput omics data (metabolomics, proteomics, or transcriptomics), the reconstruction of human GSM models of pancreatic cancer, tuberculosis infected AMs, obesity and diabetes, and neurodegenerative diseases have led to discovery of novel therapeutic targets and better understanding of the metabolic shifts14,15,16,17.

Overall, GSM models can be considered as computational and context specific (species, cells, tissue, etc.) knowledgebases capable of dissecting systemic metabolic phenomena18. A GSM model contains all annotated metabolic reactions and pathways within a specific biological system19. Using Flux Balance Analysis (FBA) and Flux Variability Analysis (FVA) the specific values and the ranges of reaction fluxes can be predicted for a given condition/timepoint. A recent computational modeling study by Gelbach et al9on M1 and M2 subtypes of human colorectal cancer cells via generation of GSM models marks a significant step towards understanding and manipulating the polarization mechanism from M2 to M1 of macrophages, related to this specific type of cancer. Other computational modeling of macrophages includes, the GSM of AM curated by Bordbar et al., from the global human model, Recon1, that was further used to model tuberculosis infected AM14,16. However, the activated states (i.e. M1 and M2 phenotypes) were not further explored. In addition, Wang et al. developed a method of creating draft metabolic models (mCADRE) of 126 human tissues including macrophages20. Overall, there is a lack of understanding of the metabolic reprogramming of AM (i.e., in healthy and activated states). While comparative analysis with other types of macrophages is highly beneficial, understanding the behavior of AMs exclusively to the lung microenvironment is critical to develop effective therapeutic targets for the lung diseases.

In this study, three context specific GSM models were developed by integrating transcriptomics data of healthy AM, and its activated phases: M1 and M2. Figure 1shows the overview of steps involved in generation, curation, and analysis of the three GSM models. Transcriptomics data integration techniques such as mCADRE20, INIT21, iMAT21and E-flux22were used for GSM model reconstruction from the global human metabolic model, Human123. The metabolic models were further curated by using a tool called OptRecon (inhouse tool, currently unpublished) that can be used to identify and resolve Thermodynamically Infeasible Cycles (TICs). These models were utilized to comprehensively analyze the metabolic landscape of activated alveolar macrophages, both in comparison to their healthy counterparts and to each other. In doing so, we highlight interesting activity in pathways such as pyruvate metabolism, glycolysis, carnitine shuttle (mitochondria) pathway, bile acid synthesis (BAS), chondroitin/heparan biosynthesis and heparan sulphate degradation as possible significant drivers of M1/M2 polarization in lungs. Additionally, we identified two sets of 30 reactions for M1 to M2 switch and vice versa, which include glycogenin formation and transport, the formation of leukotriene B4, potassium ion exchange, and movement of L-carnitine from the cytoplasm to mitochondria by developing a bilevel optimization framework (MetaShiftOptimizer). These reactions when altered (upregulated/downregulated or deleted) could more successfully shift M2 phenotype to M1; while for the reverse, it becomes more challenging. Going forward, the context specific activated AM GSM models will be used to study interaction with respiratory pathogens and identify effective therapeutic measures. In addition, a metabolic modeling framework of the relevant immune cells such as AM, Neutrophils, and Mast cells, could be developed to analyze the innate metabolic immune responses during pathogenesis in the lungs.

소개

폐포 대식세포(AMs)는

호흡기 병원체에 대한 첫 번째 방어선으로 작용하며

높은 가소성을 나타냅니다1.

병원체와의 상호작용에 따라 AMs는 서로 다른 하위 집합으로 분화될 수 있으며,

주요 두 가지 표현형은 전통적으로 활성화된(M1) 대식세포와

대안적으로 활성화된(M2) 대식세포입니다2,3.

M1 대식세포는

리포폴리사카라이드(LPS)와 Th1 염증성 사이토카인 등 미생물 성분에 의해 주로 유도되며,

세균 제거와 추가 면역 세포의 모집에 필수적입니다4,5,6,7.

반면

M2 대식세포는

인터루킨-4(IL-4), 인터루킨-13(IL-13), 면역 복합체(ICs)에 의해 촉진되며,

염증 해소 및 조직 복구와 같은 항염증 기능과 연관되어 있습니다8.

그러나

자극에 따라 M2 대식세포는 M2a, M2b, M2c, M2d의 네 가지 하위 유형으로 추가 구분될 수 있습니다9.

이 하위 유형은 매우 다양하며,

기능은 다른 면역 세포의 모집에서 조직 복구 촉진까지 다양합니다9.

이러한 통찰력을 고려할 때,

폐 미세환경의 특정 맥락에서 이러한 대사 경로의 역할을 더욱 탐구하는 것은

AM 극성화의 기능적 의미를 이해하는 데 필수적입니다10.

감염 시 폐포 대식세포의 분극은

병원체의 유형과 면역 반응의 단계에 따라 동적으로 변화할 수 있습니다4.

폐포 대식세포는

세균성 리포폴리사카라이드(LPS)나 바이러스 성분과 같은 미생물 자극에 반응하여

M1 형질로 분극합니다.

반면,

감염이 통제되고 염증 반응이 조절되어야 할 때

대식세포 분극은 M2 형질로 전환될 수 있습니다5.

AM의 형질 전환은 신체 면역 반응과 대사 조절에 중요한 역할을 합니다11,12.

그러나

분극화를 조절하는 메커니즘은 완전히 이해되지 않았습니다.

실험 연구의 어려움으로 인해,

시스템 생물학 접근법은 기존 오믹스 데이터를 활용해

맥락 특정(예: M1 또는 M2 단계) 유전체 규모 대사(GSM) 모델을 재구성하는 데 활용될 수 있습니다.

시스템 생물학은

복잡한 생물학적 시스템의 실용적 모델링과 이론적 탐구에 매우 유용함이 입증되었습니다13.

고해상도 오믹스 데이터(대사체학, 단백체학, 전사체학)의 통합을 통해

췌장암, 결핵 감염 AM, 비만 및 당뇨병, 신경퇴행성 질환의 인간 GSM 모델 재구성은

새로운 치료 표적 발견과 대사 변화에 대한 이해를 높이는 데 기여했습니다14,15,16,17.

전체적으로 GSM 모델은 계산적이며 맥락 특정적(종, 세포, 조직 등)인 지식 기반으로, 체계적 대사 현상을 분석할 수 있는 도구로 간주될 수 있습니다18. GSM 모델은 특정 생물학적 시스템 내의 모든 주석이 달린 대사 반응과 경로를 포함합니다19. 유동 균형 분석(FBA)과 유동 변동성 분석(FVA)을 통해 특정 조건/시점에서 반응 유동량의 특정 값과 범위를 예측할 수 있습니다. Gelbach 등9의 최근 계산 모델링 연구는 인간 대장암 세포의 M1 및 M2 하위 유형을 대상으로 GSM 모델을 생성하여, 이 특정 유형의 암과 관련된 대식세포의 M2에서 M1로의 분극화 메커니즘을 이해하고 조작하는 데 중요한 단계를 마련했습니다. 대식세포에 대한 다른 계산 모델링에는 Bordbar 등9이 글로벌 인간 모델 Recon1에서 추출한 AM의 GSM이 포함되며, 이는 결핵 감염된 AM을 모델링하는 데 추가로 사용되었습니다14,16. 그러나 활성화 상태(즉, M1 및 M2 표현형)는 추가로 탐구되지 않았습니다. 또한 Wang 등20은 126개 인간 조직(대식세포 포함)의 초안 대사 모델(mCADRE)을 생성하는 방법을 개발했습니다. 전체적으로 AM의 대사 재프로그래밍(즉, 건강한 상태와 활성화된 상태)에 대한 이해가 부족합니다. 다른 유형의 대식세포와의 비교 분석은 매우 유익하지만, 폐 미세환경에 특화된 AM의 행동을 이해하는 것은 폐 질환에 대한 효과적인 치료 표적을 개발하는 데 필수적입니다.

이 연구에서는 건강한 AM 및 그 활성화 단계(M1과 M2)의 전사체 데이터를 통합하여 세 가지 맥락 특이적 GSM 모델을 개발했습니다. 그림 1은 세 가지 GSM 모델의 생성, 정리, 분석 단계의 개요를 보여줍니다. mCADRE20, INIT21, iMAT21 및 E-flux22와 같은 전사체 데이터 통합 기술을 사용하여 글로벌 인간 대사 모델인 Human123에서 GSM 모델을 재구성했습니다. 대사 모델은 Thermodynamically Infeasible Cycles (TICs)를 식별하고 해결할 수 있는 도구인 OptRecon (내부 도구, 현재 미발표)을 사용하여 추가로 정리되었습니다. 이 모델들은 활성화된 폐 알베올라 대식세포의 대사 풍경을 건강 상태의 대조군 및 서로 비교하여 종합적으로 분석하는 데 활용되었습니다. 이 과정에서 피루vate 대사, 글리코lysis, 카르니틴 셔틀(미토콘드리아) 경로, 담즙산 합성(BAS), 콘드로이틴/헤파란 생합성 및 헤파란 황산염 분해와 같은 경로에서 M1/M2 분극화의 잠재적 주요 요인으로 작용할 수 있는 흥미로운 활성을 강조했습니다. 또한, M1에서 M2로의 전환 및 그 반대의 과정을 위한 두 세트의 30개 반응을 식별했습니다. 이 반응에는 글리코겐 형성 및 운반, 류코트리엔 B4 형성, 칼륨 이온 교환, 그리고 MetaShiftOptimizer라는 이중 최적화 프레임워크를 개발하여 세포질에서 미토콘드리아로 L-카르니틴의 이동이 포함됩니다. 이 반응들이 변화(상향 조절/하향 조절 또는 삭제)될 경우 M2 형질을 M1로 전환하는 데 더 효과적일 수 있으며, 반대의 경우 더 어려워집니다. 향후, 호흡기 병원체와의 상호작용을 연구하고 효과적인 치료 조치를 식별하기 위해 맥락 특이적 활성화된 AM GSM 모델이 사용될 것입니다. 또한, 폐 병리 과정에서 선천성 대사 면역 반응을 분석하기 위해 AM, 중성구, 비만세포와 같은 관련 면역 세포의 대사 모델링 프레임워크가 개발될 수 있습니다.

Fig. 1

Schematic of the workflow for the generation of healthy and activated Alveolar macrophages and the steps involved in analysis of the metabolic shift during polarization. (a) The schematic of the steps involved in the generation of context specific GSM models by integrating omics data into global human reconstruction, Human1. (b) Optimization techniques such as Flux Balance Analysis and Flux Variability Analysis were used to obtain flux ranges to analyze the metabolic behavior of each of the phases. (c) A t-SNE plot that shows the flux distribution of M1 and M2 phase and Modified M2 (M2 after incorporation of changes as predicted by MetaShiftOptimizer) to obtain the shift from M2 to M1 phase. The same process is applied to obtain the shift from M1 to M2 phase as describe in Results and Methods. .

건강한 활성화된 폐포 대식세포의 생성 과정 및 분극화 과정에서 대사 변화 분석에 포함된 단계의 흐름도.

(a) 오믹스 데이터를 글로벌 인간 재구성 모델 Human1에 통합하여 맥락 특이적 GSM 모델을 생성하는 단계의 흐름도.

(b) 각 단계의 대사 행동을 분석하기 위해 유동 균형 분석(Flux Balance Analysis) 및 유동 변동성 분석(Flux Variability Analysis)과 같은 최적화 기법을 사용하여 유동 범위를 도출했습니다.

(c) M1 및 M2 단계와 MetaShiftOptimizer로 예측된 변화를 반영한 수정된 M2(M2)의 유동 분포를 보여주는 t-SNE 플롯으로, M2에서 M1 단계로의 전환을 확인합니다. 동일한 과정은 결과 및 방법 섹션에 설명된 대로 M1에서 M2 단계로의 전환을 얻기 위해 적용됩니다.

Results

Metabolic model reconstruction of alveolar macrophage (AM) metabolism

To develop the context specific metabolic models of healthy AM, M1 phase, and M2 phase, we needed to integrate gene expression‎ values of metabolic genes24,25,26onto the Human123 model. These transcriptomic profiles were acquired from GEO databases (GSE8823, GSE40885, and GSE41649 for AM, M1 phase and M2 phase, respectively). These datasets effectively capture the biology of various states of alveolar macrophages, as demonstrated by Becker et al.27, where the GSE40885 and GSE41649 datasets were utilized to identify macrophage polarization signatures in the lungs. Consequently, we utilized these datasets to perform gene set enrichment analysis using GSEA desktop software28 to gain an initial understanding of the differences in pathway activities based solely on transcriptomic profile. The analysis captured the activity of glycolysis, bile acid synthesis, and peroxisomal pathways to be enriched for M1, while cholesterol metabolism and protein secretion were enriched for M2. In addition, OXPHOS and Fatty acid metabolism were found to be enriched in both the phenotypes (the details on the GSEA results are available in the GitHub repository). While gene set enrichment analysis provided an initial insight into pathways that could be important in M1 and M2 phase, we used GSM models to further explore how the pathway level transcriptomic difference can be manifested into metabolic differences and capabilities of each AM phenotype.

Among various methods available to develop context specific models via the integration of omics data, INIT (Integrative Network Inference for Tissues)29, mCADRE (Context-specificity Assessed by Deterministic Reaction Eval‎uation)20, iMAT (integrative Metabolic Analysis Tool)21, and E-flux22were used in this study. While iMAT, INIT, and mCADRE are switch approaches in which reactions are turned on/off depending on gene expression‎ levels, E-flux is a valve approach in which each reaction is constrained based on the gene expression‎ values. iMAT madidates the presence/absence of a specific reaction depending on the relevant gene(s) having higher expression‎ levels at a specific condition21. INIT considers evidence from various sources, such as protein levels and gene expression‎, assigning weights to reactions in a way that maximizes network consistency with available data while allowing for flexibility in metabolite secretion or accumulation29. mCADRE algorithm creates a model by removing the non-core metabolic reactions while ensuring the consistency of high confidence set of core reactions with the ability to produce essential metabolite from simple precursors such as glucose20. On the other hand, a valve approach such as E-flux, uses gene expression‎ levels to regulate the flux of the corresponding reactions22.

The healthy AM model obtained upon the implementation of iMAT consists of 4,554 reactions (governed by 2,173 metabolic genes) and 3,967 metabolites (2,003 unique) distributed across eight intracellular compartments (Extracellular, Peroxisome, Mitochondria, Cytosol, Lysosome, Endoplasmic reticulum, Golgi apparatus, Inner mitochondria, and Nucleus); while the model generated by E-flux consists of 8,073 reactions and 5,380 metabolites (2,823 unique) across these eight compartments with the same number of metabolic genes. Similar comparisons were performed between the models obtained from INIT and mCADRE and are shown in Fig. 2(a, b, and c). INIT and mCADRE algorithms have the risk of either underestimating or overestimating the metabolic capabilities of any cellular systems which could potentially be problematic for human cells and can lead to divergence from their true physiological behavior20,29. In addition, the sensitivity of iMAT approach to user defined threshold typically leads to higher number of reactions to be omitted or sometimes leads to exclusion of important reactions from the pruned model. In our case, these switch approaches resulted a pruned model capable of producing biomass but failed to include important pathways such as NO production, glycerolipid metabolism, heme synthesis, and porphyrin metabolism or consisted of several incomplete pathways such as fatty acid oxidation (FAO), fatty acid synthesis (FAS), and amino sugar and nucleotide metabolism. The E-flux-generated model, on the other hand, consists of comparatively more comprehensive list of reactions, metabolites, and ensured the expected activity of all the important pathways mentioned earlier. Additional information on the models built by each of these methods can be found in Supplementary Table 1, 2 and 3 along with the pathway distribution figure in Supplementary Data 1.

Fig. 2

Four different integration methods, namely iMAT, E-flux, INIT, and mCADRE are implemented for the context specific GSM model reconstruction from Human1. (a) The number of reactions and metabolites incorporated in each situation (healthy, M1 phase and M2 phase) by different methods is displayed through spider plots. In all three plots orange represents the number of metabolites and maroon represents the number of reactions. Each number of reactions and metabolites is scaled by 1000 (for example 3,967 is plotted as 3.697 and so on) to have better visualization. Methods like iMAT and mCADRE generate significantly smaller GSM models in comparison to INIT and E-flux. (b) A bar graph displaying the number of reactions carrying fluxes for each state of alveolar macrophage shows that E-flux GSM models has the highest number of reactions carrying fluxes in all of the cases. Hence, for the purpose of this study, we used E-flux generated GSM models.

Full size image

Similar observations were obtained while implementing the above-mentioned algorithms with the expression‎ values of 2,951 and 2,390 metabolic genes to reconstruct GSM models for M1 and M2 phase respectively (additional information on the gene expression‎ values, distribution of pathways, reactions, and metabolites for all the models are available in Supplementary Data 1). Hence, the models developed with the implementation of E-flux allowed us to exhaustively investigate the metabolic shift occurring during polarization. The details on the GSM model comparison for all the generated models is available in Supplementary Data 1.

Validation of the metabolic capabilities of phase specific AM GSM models

The E-flux models were next validated to ensure their ability to simulate biologically significant processes by reproducing important metabolite production rates and characteristic behavior as reported in literature. The biomass function/reaction of these models was developed based on the AM cellular maintenance requirements such as proteins, lipids, DNA repair, ATP maintenance, and RNA turnover (adapted from Bordbar et al. 2010)16. AMs are tissue resident macrophages that populate the lung environment during birth and typically last for the lifespan of the individuals16. Flux Balance Analysis (FBA) is a mathematical approach to predict the flow of metabolites through a GSM model with a specific goal/objective, typically maximizing of growth/biomass. Since AM population remains stable under steady state condition, the choice of growth for AM models does not make sense. Hence, for all the simulations performed by our AM models, we optimized the production rates of key metabolites, while maintaining the maximal level of biomass (i.e., 0.03 h−1). The healthy AM model was optimized for ATP production and NO production that yielded the flux of 0.6 mol/g cell DW/h and 0.03 mol/g cell DW/h respectively. These in silico values were very close to the values of 0.71 mol/g cell DW/h and 0.037 mol/g cell DW/h, respectively, as obtained from in vitro experiment16,30. Hence, the healthy AM model obtained via E-flux algorithm is capable of accurately recapitulating experimentally reported production rates of important metabolites.

Similar to the healthy AM model, the activated phase models should also be able to recapitulate the relevant metabolic reprogramming. Based on our in silico observations, the maximal flux range for growth (i.e., biomass flux) was found to be different across all states (healthy, M1 and M2 phase) with the minimum value found (0.03 h−1) for the healthy state. Hence, the biomass level was kept at this minimum value for the activated phases and was used as a maintenance function only. Despite lack of experimental evidence of production rates of any specific metabolites, the comparative enhanced/inhibited pathway and transcriptional activity were widely reported in literature for M1 and M2 phase31,32. To this end, with the help of the flux ranges obtained through FVA with the healthy AM GSM model as the base/control, the M1/M2 phase reaction fluxes were compared. FVA is an extension of FBA that identifies the range of possible fluxes for each reaction in a metabolic model33,34. These fluxes were categorized into five different categories, namely, complete overlap: widened flux space, complete overlap: shrunk flux space, partial overlap: increase, no overlap: definite increase in forward direction, and no overlap: definite increase in reverse direction. Complete overlap with increase (widened flux space) indicates a collective upregulation in the reaction fluxes, and similarly complete overlap with decrease (shrunk flux space) indicates a collective downregulation in the reaction fluxes. Partial overlap indicates that there are some flux values within the ranges that are shared, while other values are unique to each range. No overlap scenarios (in forward or reverse direction) represent completely unique flux ranges. Figure 3 shows the flux distribution in seven different pathways that play an important role during AM polarization including glycolysis, OXPHOS, PPP, and TCA cycle.

상위 단계별 AM GSM 모델의 대사 능력 검증

E-flux 모델은 문헌에 보고된 중요한 대사체 생산 속도와 특이적 행동을 재현함으로써 생물학적으로 의미 있는 과정을 시뮬레이션할 수 있는지 검증되었습니다. 이 모델의 생물량 기능/반응은 AM 세포 유지 요구사항인 단백질, 지질, DNA 복구, ATP 유지, RNA 회전율(Bordbar et al. 2010에서 수정됨)을 기반으로 개발되었습니다.16 AM은 출생 시 폐 환경에 정주하는 조직 거주성 대식세포로, 일반적으로 개체의 수명 동안 존재합니다16. 유동 균형 분석(FBA)은 특정 목표/목적(일반적으로 성장/생물량 최대화)을 달성하기 위해 GSM 모델을 통해 대사체 흐름을 예측하는 수학적 접근법입니다. AM 인구 밀도는 안정 상태에서 안정적으로 유지되므로, AM 모델에 성장 목표를 설정하는 것은 의미가 없습니다. 따라서 우리 AM 모델에서 수행된 모든 시뮬레이션에서 우리는 생물량 최대 수준(즉, 0.03 h−1)을 유지하면서 핵심 대사물질의 생산 속도를 최적화했습니다. 건강한 AM 모델은 ATP 생산과 NO 생산을 최적화하여 각각 0.6 mol/g 세포 DW/h와 0.03 mol/g 세포 DW/h의 유량을 나타냈습니다. 이 in silico 값은 in vitro 실험16,30에서 얻은 0.71 mol/g cell DW/h와 0.037 mol/g cell DW/h와 매우 유사했습니다. 따라서 E-flux 알고리즘을 통해 얻은 건강한 AM 모델은 실험적으로 보고된 중요한 대사산물 생산 속도를 정확히 재현할 수 있습니다.

건강한 AM 모델과 유사하게, 활성화 단계 모델도 관련 대사 재프로그래밍을 재현할 수 있어야 합니다. 우리 in silico 관찰에 따르면, 성장(즉, 생물량 유동)의 최대 유동 범위는 모든 상태(건강, M1 및 M2 단계)에서 서로 다르게 나타났으며, 건강한 상태에서 최소 값(0.03 h−1)이 관찰되었습니다. 따라서 활성화 단계에서는 생물량 수준을 이 최소 값으로 유지했으며 유지 기능으로만 사용되었습니다. 특정 대사물의 생산 속도에 대한 실험적 증거가 부족함에도 불구하고, M1 및 M2 단계에서 비교적 강화되거나 억제된 경로 및 전사 활성은 문헌에서 널리 보고되었습니다31,32. 이를 위해 건강한 AM GSM 모델을 기준/통제 모델로 사용하여 FVA를 통해 얻은 유동 범위(flux ranges)를 바탕으로 M1/M2 단계 반응 유동량을 비교했습니다. FVA는 대사 모델 내 각 반응의 가능한 유동량 범위를 식별하는 FBA의 확장 방법입니다33,34. 이 유동은 다섯 가지 다른 범주로 분류되었습니다. 즉, 완전 중첩: 유동 공간 확장, 완전 중첩: 유동 공간 축소, 부분 중첩: 증가, 중첩 없음: 전방 방향에서 확실한 증가, 중첩 없음: 역방향에서 확실한 증가입니다. 완전 중첩과 증가(유동 공간 확장)는 반응 유동의 집단적 상향 조절을 나타내며, 마찬가지로 완전 중첩과 감소(유동 공간 축소)는 반응 유동의 집단적 하향 조절을 의미합니다. 부분 중첩은 일부 유동 값이 공유된 범위 내에 존재하며, 다른 값은 각 범위에서 고유하다는 것을 의미합니다. 중첩 없음 시나리오(전방 또는 역방향)는 완전히 고유한 유동 범위를 나타냅니다. 그림 3은 AM 극화 과정에서 중요한 역할을 하는 7개의 서로 다른 경로(글리코lysis, OXPHOS, PPP, TCA 사이클 등)의 유동 분포를 보여줍니다.

Fig. 3

Alveolar Macrophage acquires unique metabolic characteristics depending upon the phenotype. In the M1 phase, the reactions of glycolysis are enhanced which are highlighted by the green arrows and the PPP reactions which is a major contributor for NAPH production is also enhanced. Similarly, the pathways highlighted by yellow arrows in M2 phase are found to be enhanced. Each pie chart represents metabolic reprogramming of AM in the specific pathway in either M1 phase or M2 phase. Each component of the pie chart represents one of the four categories as color coded in the figure. The associated percentage in the pie chart represents the percentage of overall reactions of a specific pathway falling into each of the categories.

폐포 대식세포는 표현형에 따라 독특한 대사 특성을 획득합니다. 

M1 단계에서는 글리코lysis 반응이 강화되며, 이는 녹색 화살표로 표시되어 있습니다. 또한 NAPH 생성에 주요 기여 요인인 PPP 반응도 강화됩니다. 마찬가지로 M2 단계에서 노란색 화살표로 표시된 경로들도 강화된 것으로 확인되었습니다. 각 파이 차트는 M1 단계 또는 M2 단계에서 AM의 특정 경로별 대사 재프로그래밍을 나타냅니다. 파이 차트의 각 구성 요소는 그림에 색상으로 표시된 네 가지 카테고리 중 하나를 나타냅니다. 파이 차트에 표시된 관련 백분율은 특정 경로의 전체 반응 중 각 카테고리에 속하는 반응의 비율을 나타냅니다.

Full size image

As reported in literature for M1 phase, we observe increased activity (complete overlap and partial overlap) for glycolysis and gluconeogenesis pathway30,31. Pentose Phosphate Pathway (PPP) shows 48% of the reactions with increased fluxes which include important reactions such as formation of ribulose 5-phosphate and its conversion to ribose 5-phosphate with production of NADPH. The PPP, through the production of NAPDH, is reported to contribute to cellular functions that greatly impact cellular energy balance and overall health30,31. We also found the increased reaction fluxes for succinate and itaconate in the TCA cycle which act as the signaling molecules for inflammatory responses during M1 phenotype30,31. Similarly, multiple studies by Wang et al35, Zhang et al36, Wculek et al37, and Hooftman et al38provide evidence regarding the relationship between itaconate and macrophage polarization32. In line with some of the recent studies11,12 suggesting that M1 and M2 macrophages primarily derive their energy from oxidative phosphorylation (OXPHOS) and the TCA cycle, we also noted increased flux activities in reactions that transport reducing equivalents (in the form of NADH) from the cytoplasm into the mitochondria. This includes the reversible conversions of fumarate to malate and citrate to isocitrate, as well as the oxidation of acetyl-CoA. Similarly, for the M2 phenotype, enhanced activity in OXPHOS, TCA cycle along with fatty acid oxidation (FAO) was observed. Additionally, succinate and itaconate production through TCA cycle, was found to be inhibited, contrary to what we found for M1 phenotype.

Next, we used Flux Sum Analysis (FSA)39,40, as a proxy of metabolite pool size, to observe the differences between M1 and M2 phenotypes. FSA was conducted for metabolites such as ATP, NADPH, succinate, itaconate, and NO. Among these metabolites, the overall ATP and NADPH pool sizes were found to be significantly higher in the M2 phenotype with the major contribution from OXPHOS. However, glycolysis seems to be a significant but not the only contributor in M1 the phenotype. Similarly, in line with literature reports31,35,38, the pool sizes for succinate and itaconate showed higher levels during the M1 phase compared to the M2 phase. The results of FSA are available in Supplementary Data 7 in detail.

Hence, the metabolic traits of each phase as predicted by our models are well supported by literature and thus establish the credibility/predictability of our context-specific GSM models of the activated phases. Supplementary Table 4 of Supplementary Data 1 shows the full list of metabolic capabilities of each GSM models. Additional information including the lists of all the metabolites and reactions of the final models are available in Supplementary Data 2.

De novo metabolic reprogramming in AM polarization mechanismActivated phases in comparison to healthy AM

AM in the lungs act as the first line of defense against respiratory pathogens since these phagocytize pollutants and pathogens that act as a trigger to activate an innate immune response11. M1 and M2 macrophages acquire distinct phenotypes which are usually driven by different stimuli. M1 macrophages are characterized by their high production of pro-inflammatory cytokines (e.g., TNF-α, IL-6, and IL-1β), reactive oxygen species (ROS), and nitric oxide (NO)5. In contrast, M2 macrophages, known as "alternatively activated," are associated with anti-inflammatory and tissue-repair functions7. They are typically induced by stimuli such as IL-4 and IL-13 and are characterized by the secretion of anti-inflammatory cytokines like IL-10 and TGF-β. M2 macrophages support wound healing, tissue remodeling, and resolution of inflammation, promoting Th2-type responses8. While earlier literature reported distinct upregulated glycolysis and inhibited OXPHOS/TCA cycle in M1 macrophages and upregulated OXPHOS/TCA cycle along with enhanced glutaminolysis and polyamine synthesis in M2 macrophages9,10,11. However, recent studies, such as those by Woods et al. 2020 and Khaing et al. 2020, have suggested that M1 macrophages may not exclusively depend on glycolysis for energy metabolism in comparison to bone marrow-derived macrophages (BMDMs)41,42. In this study, particular attention and effort were directed toward identifying the activity of pathways in activated AMs with respect to healthy AM. Additionally, the pathways beyond the central carbon metabolism such as pyruvate metabolism, arachidonic acid metabolism, chondroitin/heparan sulphate biosynthesis, and heparan sulphate (HS) degradation were exhaustively investigated. Despite the increasing interest in AM polarization and their unique contribution to the progression and suppression of diseases, there is a lot to be learned about the role of the above-mentioned pathways in lung pathogenesis. Hence, upon validating the models, we compared the metabolic shift in activated AMs to the healthy state by using the Flux Variability Analysis (FVA) as explained earlier in the model validation section.

We observed that glycolysis did not show distinct inhibition in M2 phase in addition to the upregulated OXPHOS/TCA, while the OXPHOS/TCA activity was also not completely down in M1 phase. Based on the FSA predictions (as mentioned in the previous section), the energy metabolism seems to be dependent on both the glycolysis and OXPHOS/TCA, while the degree of contribution of each pathway varies based on the phenotype. One possible reason for the absence of distinct up/downregulation of these key pathways could also be the diversity of M2 subtypes. Recent findings indicate that M2 phenotypes is further divided into several subtypes, out of which M2d subtype is capable of proinflammatory responses43. Hence, we postulate all M2 subtype population may not exhibit inhibited glycolysis activity based on the in silico predictions. On the other hand, upregulated activity of carnitine shuttle pathway is unique to M2 phase as we did not observe similar activity in M1 phase. In fact, for M1 phase, over 80% of the reactions from carnitine shuttle (mitochondria) metabolism showed decrease in flux ranges. Figure 4 portrays the complex regulatory pathways in AM polarization and shows the pathways that are up- or down-regulated in either state.

문헌에서 보고된 M1 단계에서와 마찬가지로, 우리는 글리코lysis 및 글루코네오제네시스 경로에서 활동 증가(완전 중첩 및 부분 중첩)를 관찰했습니다30,31. 펜토스 인산 경로(PPP)는 반응의 48%에서 증가된 유동량을 보이며, 이는 리불로즈 5-인산 형성 및 NADPH 생성과 함께 리불로즈 5-인산으로의 전환과 같은 중요한 반응을 포함합니다. PPP는 NADPH 생성을 통해 세포 에너지 균형과 전반적인 건강에 큰 영향을 미치는 세포 기능에 기여한다고 보고되었습니다30,31. 우리는 TCA 회로에서 수크신산과 이타콘산의 반응 유동량이 증가했으며, 이는 M1 형질에서 염증 반응의 신호 분자로 작용합니다30,31. 마찬가지로 Wang 등35, Zhang 등36, Wculek 등37, Hooftman 등38의 다수 연구는 이타콘산과 대식세포 분극화 간의 관계를 입증했습니다32. 최근 연구11,12에서 M1과 M2 대식세포가 주로 산화적 인산화(OXPHOS)와 TCA 회로에서 에너지를 얻는다는 제안과 일치하여, 우리는 세포질에서 미토콘드리아로 환원 등가물(NADH 형태)을 운반하는 반응의 유동 활동이 증가했음을 관찰했습니다. 이에는 푸마레이트에서 말레이트로의 가역적 전환, 시트레이트에서 이소시트레이트로의 전환, 아세틸-CoA의 산화 등이 포함됩니다. 마찬가지로 M2 형질에서는 OXPHOS, TCA 회로 및 지방산 산화(FAO)의 활동이 증가했습니다. 또한 TCA 회로를 통해 생성되는 수크신산과 이타콘산 생산이 억제되었으며, 이는 M1 형질에서 관찰된 결과와 반대였습니다.

다음으로, 대사체 풀 크기의 대리 지표로 Flux Sum Analysis (FSA)39,40를 사용하여 M1과 M2 형질 간의 차이를 관찰했습니다. FSA는 ATP, NADPH, 수산화물, 이타콘산, NO와 같은 대사체에 대해 수행되었습니다. 이 대사체 중 전체 ATP 및 NADPH 풀 크기는 M2 형질에서 유의미하게 높았으며, 이는 주로 OXPHOS에 의해 기여되었습니다. 그러나 M1 형질에서는 글리코lysis가 중요한 기여 요인이지만 유일한 요인은 아닙니다. 마찬가지로 문헌 보고31,35,38과 일치하게, succinate와 itaconate의 풀 크기는 M1 단계에서 M2 단계보다 높은 수준을 보였습니다. FSA 결과는 보충 자료 7에 상세히 수록되어 있습니다.

따라서, 우리 모델이 예측한 각 단계의 대사 특성은 문헌과 잘 일치하며, 이는 활성화 단계에 대한 맥락 특정 GSM 모델의 신뢰성과 예측 가능성을 확립합니다. 보완 자료 1의 보완 표 4에는 각 GSM 모델의 대사 능력 전체 목록이 표시되어 있습니다. 최종 모델의 모든 대사체 및 반응 목록을 포함한 추가 정보는 보완 자료 2에 수록되어 있습니다.

AM 분극화 메커니즘에서의 de novo 대사 재프로그래밍활성화 단계와 건강한 AM의 비교

폐 내 AM은 호흡기 병원체에 대한 첫 번째 방어선으로 작용하며, 오염물질과 병원체를 식균작용으로 제거하여 선천성 면역 반응을 활성화하는 트리거 역할을 합니다11. M1과 M2 대식세포는 서로 다른 자극에 의해 유도되는 독특한 표현형을 획득합니다. M1 대식세포는 프로염증성 사이토킨(예: TNF-α, IL-6, IL-1β), 활성산소종(ROS), 및 질소산화물(NO)의 높은 생산으로 특징지어집니다5. 반면, '대체 활성화'로 알려진 M2 대식세포는 항염증 및 조직 재생 기능과 연관되어 있습니다7. 이들은 IL-4 및 IL-13과 같은 자극에 의해 유도되며, IL-10 및 TGF-β와 같은 항염증성 사이토카인의 분비로 특징지어집니다. M2 대식세포는 상처 치유, 조직 재구성, 염증 해소 등을 지원하며 Th2형 반응을 촉진합니다8. 이전 연구에서는 M1 대식세포에서 글리코lysis가 증가하고 OXPHOS/TCA 사이클이 억제되며, M2 대식세포에서는 OXPHOS/TCA 사이클이 증가하고 글루타민olysis 및 폴리아민 합성이 강화된 것으로 보고되었습니다9,10,11. 그러나 Woods 등(2020)과 Khaing 등(2020)의 최근 연구는 M1 대식세포가 골수 유래 대식세포(BMDMs)에 비해 에너지 대사에서 글리코lysis에 독점적으로 의존하지 않을 수 있음을 제안했습니다41,42. 본 연구에서는 활성화된 AM과 건강한 AM 간의 경로 활성도를 식별하는 데 특별히 주목하고 노력했습니다. 또한 중앙 탄소 대사 beyond the central carbon metabolism을 넘어 피루vate 대사, 아라키돈산 대사, 콘드로이틴/헤파란 황산염 생합성, 헤파란 황산염(HS) 분해와 같은 경로를 포괄적으로 조사했습니다. AM 분극화와 질병의 진행 및 억제에 대한 독특한 기여에 대한 관심이 증가하고 있음에도 불구하고, 위에서 언급된 경로가 폐 병리학에 미치는 역할에 대해 아직 많은 것이 알려지지 않았습니다. 따라서 모델을 검증한 후, 모델 검증 섹션에서 설명된 대로 유동성 변동 분석(FVA)을 사용하여 활성화된 AM의 대사 변화와 건강한 상태를 비교했습니다.

우리는 M2 단계에서 OXPHOS/TCA의 증가와 함께 글리코lysis가 명확한 억제를 보이지 않았으며, M1 단계에서도 OXPHOS/TCA 활성이 완전히 감소하지 않았음을 관찰했습니다. FSA 예측(이전 섹션에서 언급된 바와 같이)에 따르면, 에너지 대사 과정은 글리코lysis와 OXPHOS/TCA에 모두 의존하지만, 각 경로의 기여도는 표현형에 따라 달라집니다. 이러한 핵심 경로의 명확한 상향/하향 조절이 관찰되지 않은 한 가지 가능한 이유는 M2 하위 유형의 다양성일 수 있습니다. 최근 연구 결과에 따르면 M2 표현형은 여러 하위 유형으로 추가로 구분되며, 이 중 M2d 하위 유형은 염증 반응을 유발할 수 있습니다43. 따라서 우리는 모든 M2 하위 유형 인구군이 in silico 예측에 따라 글리코lysis 활성 억제를 나타내지 않을 수 있다고 가정합니다. 반면, 카르니틴 셔틀 경로의 활성화는 M2 단계에 독특하며, M1 단계에서는 유사한 활성을 관찰하지 못했습니다. 실제로 M1 단계에서는 카르니틴 셔틀(미토콘드리아) 대사 경로의 반응 중 80% 이상이 유동량 감소 범위를 보였습니다. 그림 4는 AM 분극화의 복잡한 조절 경로를 보여주고 각 상태에서 상향 또는 하향 조절된 경로를 표시합니다.

Fig. 4

Important AM pathways. Glycolysis, TCA cycle, and OXPHOS play major roles in energy production with the help of pathways such as the carnitine shuttle (mitochondria), which shows enhanced activity during the anti-inflammatory phase. On the contrary, Bile Acid Synthesis, and Arachidonic Acid Metabolism are heightened to induce acidic conditions to minimize pathogen survival. Pyruvate Metabolism play key roles in the immune response of the cell.

Full size image

The activity of the bile acid synthesis and arachidonic acid metabolism was found to be enhanced exclusively in M1 phenotype. Specifically, the reactions involved in the formation of oxylipins such as 5,15-DiHETE and 5,6-Ep-15S-HETE from arachidonic acid metabolism were found to have increased fluxes44,45,46. Since the M1 model is LPS induced, the presence of oxylipins and bile acids are in accordance with the immune response of pro-inflammatory AM. Similarly, the upregulated activity of chondroitin/heparan sulfate biosynthesis in the M1 model can also be linked to the proinflammatory responses. The formation of heparan sulfate is a crucial step for the recruitment, adhesion, crawling, and transmigration of leukocytes from the circulation to the site of inflammation47,48,49. On the other hand, increased heparan sulphate (HS) degradation was observed in the M2 model. Although the role of formation and degradation of HS has been a topic of interest during lung injury and inflammation, it has not been systematically explored47,48,49. Hence, the in silico activity of Chondroitin/heparan biosynthesis postulates that this pathway is enhanced due to inflammatory response in the M1 phase, while the HS degradation is enhanced in order to aid to the versatile function in M247,48,49.

In addition to the metabolite mentioned above, pyruvate is an important metabolite involved in maintenance of cellular and immune function. It is a key mediator connecting several important pathways such as glycolysis, TCA cycle, amino acid metabolism (such as arginine)50,51. Despite being such a key modulator, the pyruvate metabolism was found to be inhibited both in the M1 and M2 phase. We observed that the reactions contributing to the direct formation of pyruvate were mostly inhibited in both of these activated macrophages. While the key reaction such as PEP to pyruvate (at the end of glycolysis) and pyruvate to oxaloacetic acid (OAA) (at the beginning of TCA) maintained high fluxes, other reactions that contribute to pyruvate production via different mechanisms (e.g., from lactate or methylglyoxal) was inhibited. Hence, the activity of pyruvate metabolism emerged as a possible polarization driver in AMs which is further explored in the next section. Additional information on the flux comparison is available in Supplementary Data 3.

Activated phases in comparison to each other

We next compared the FVA solution of activated phenotypes to each other in similar manner as explained above. Earlier (while comparing healthy AM fluxes with activated phases), we reported that despite being a key intermediatory metabolite, the metabolic shift in AM resulted in limited pyruvate production in both M1 and M2 phase with respect to healthy AM. However, while comparing the M1 fluxes with those from M2, we observed all the reactions contributing to pyruvate production in the cytoplasm were inhibited (complete overlap, shrunk flux space), whereas the mitochondrial reactions are enhanced (complete overlap, widened flux space). In Fig. 5, it is shown that in the cytoplasm, only one reaction (malate to pyruvate) has enhanced fluxes with higher fluxes toward the production of D-Lactate. On the contrary, enhanced fluxes were observed in multiple mitochondrial reactions that indicate increased pyruvate and L-lactate production. The pyruvate produced in mitochondria is directly used up for OAA production which promotes OXPHOS and TCA cycle, while lactate plays an important role in the maintenance of acid–base balance in the cell and in the maintenance and resolution of inflammation50,51. A study by Abusalamah et al. (2020) suggests that incorporating pyruvate as sodium pyruvate in growth media for macrophages inhibited immune response of the cell and also had positive impact on the bacterial growth51. Hence, the reprogramming of the pathway as suggested by our in silico analysis shows that the cells actively inhibit multiple reactions that contribute to pyruvate production while upregulating specific reactions that channel pyruvate into mitochondria for energy production. Further experimental studies in human AM could establish if pyruvate metabolite is not only an important factor but also recognizes the regulation of pyruvate metabolism as a key step in pathogenesis.

Fig. 5

Pyruvate Metabolism activity in activated phase M2 when compared to M1 phase. The reactions indicated by green arrow are enhanced in the M2 phase and show clear redirection of metabolic activity through pyruvate production.

Full size image

In addition to pyruvate production reprogramming, glycogen emerged as a significant factor in acquiring specific phenotype. We observed the category with "definite increase in forward and reverse direction" consisted of mainly reactions related to glycogen production/consumption. While this observation could also be due to glycogen being a part of the biomass function, the potential role of glycogen in modulating immune responses is also highly probable, especially due to biomass reaction merely acting as a maintenance function for all of our analysis. The upregulated reactions mainly include the production of glycogenin G8 from glycogenin in cytoplasm that results in glycogenin G4G4, which is a critical precursor in the synthesis of glycogen. To further explore the M1 and M2 metabolic models for highly impactful reactions, we constrained the models by turning off the reactions completely (also known as reaction knockout) or by severely limiting the flux of each of these reactions. We found reaction involved in conversion of UTP and alpha-D-galactose-1-phosphate into UDP-galactose within the cytoplasm as critical for sustaining vital connections in the metabolic models during activated phenotypes. When these specific reactions were completely knocked out, we observed severe negative impact in the GSM models, while constraining or relaxing the reaction fluxes resulted in more desirable output. Hence, in the next section we will discuss our bilevel optimization framework to identify reactions that could be potential drivers of AM polarization52,53.

Finally, to analyze whether the pathway activities observed so far could be corroborated from a metabolite perspective, we also conducted a metabolite-centric analysis that identifies highly significant metabolites based on the transcriptomics date of each phenotypic state39. To this end, metabolites with positive normalized score were reported as significant metabolites in this study. The algorithms identified numerous metabolites from leukotriene metabolism, arachidonic acid metabolism, and chondroitin/heparan sulphate biosynthesis and degradation as significant for the M1 phase, while the major metabolites for the M2 phase were from glycolysis, OXPHOS and TCA cycle. This observation helped further substantiate the in silico predictions as mentioned earlier. The complete list of significant metabolites is available in the Supplementary Data 4.

Minimal modifications for AM polarization shift

It is extremely crucial for AM to acquire right phenotype for proper immune response. The imbalance in the M1/M2 cells can be deleterious to the lungs that can cause prolonged and unwanted inflammation in the absence of the process that shuts it down54. In addition, without the necessary inflammation, AMs cannot effectively activate other immune cells to fight invading organisms55. For example, the interaction between tuberculosis and AM is sometimes reported to promote M2 cells as opposed to M1, and reports on progression of cancer cells also mention the positive role of the M2 phenotype56. Hence it is very important for AM to shift toward the phenotype which is best suited to fight the invading pathogens. To this end, we formulated a bilevel optimization framework named as MetaShiftOptimizer (as described in the Methods) that predicts a set of reactions when altered (i.e., up-/down-regulation and minimal deletion) can force the shift from pro-inflammatory (M1) state to anti-inflammatory (M2) state and vice versa. The identified reactions in both the scenarios can be linked to important immune responses and have been linked to signaling pathways involved in macrophage polarization55. For example, one of the reactions identified by MetaShiftOptimzer is involved in the production and storage of glycogenin. Glycogen metabolism is known to play a significant role in regulation of macrophage mediated acute inflammatory responses and can be linked to several signaling networks, particularly those involved in metabolic regulation and cellular energy balance such as AMP-activated Protein Kinase (AMPK) Pathway, and mTOR pathway55,57. Similarly, Leukotriene B4 is also known as a potent inflammatory mediator that is closely linked to NF-KB, MAPK, PI3K/AKT pathways55,57. The possible links of different reactions and pathways to various signaling mechanism is illustrated in Fig. 6.

AM 극성 전환을 위한 최소한의 수정

AM이 적절한 면역 반응을 위해 올바른 표현형을 획득하는 것은 극히 중요합니다.

M1/M2 세포의 불균형은

이 과정을 차단하는 메커니즘이 결여될 경우

폐에 해로운 영향을 미쳐 장기적이고 원치 않는 염증을 유발할 수 있습니다54.

또한,

필요한 염증이 결여되면

AM은 침입한 병원체와 싸우기 위해 다른 면역 세포를 효과적으로 활성화할 수 없습니다55.

예를 들어,

결핵과 AM의 상호작용은 때로는 M1 대신 M2 세포를 촉진한다는 보고가 있으며,

암 세포의 진행에 대한 연구에서도 M2 형질의 긍정적인 역할이 언급되었습니다56.

따라서

AM이 침입한 병원체와 싸우는 데 가장 적합한 형질로 전환되는 것이 매우 중요합니다.

이를 위해 우리는 방법론에서 설명된 대로 MetaShiftOptimizer라는 이중 최적화 프레임워크를 수립했습니다.

이 프레임워크는

반응을 변경(즉, 상향/하향 조절 및 최소 삭제)할 때

염증성(M1) 상태에서 항염증성(M2) 상태로 전환하거나

그 반대로 강제할 수 있는 반응 집합을 예측합니다.

두시나리오에서 식별된 반응은 중요한 면역 반응과 연관되어 있으며,

대식세포 분극화에 관여하는 신호 전달 경로와 연결되어 있습니다55.

예를 들어,

MetaShiftOptimizer가 식별한 반응 중 하나는

글리코겐의 생산 및 저장 과정에 관여합니다.

글리코겐 대사는

대식세포 매개 급성 염증 반응의 조절에 중요한 역할을 하며,

대사 조절 및 세포 에너지 균형과 관련된 신호 전달 네트워크,

특히 AMP-활성화 단백질 키나제(AMPK) 경로 및 mTOR 경로55,57와 연결될 수 있습니다.

同様に,

류코트리엔 B4는

NF-KB, MAPK, PI3K/AKT 경로와 밀접하게 연관된 강력한 염증 매개체로 알려져 있습니다55,57.

다양한 반응과 경로가 다양한 신호 전달 메커니즘과 연결될 수 있는 가능성은

그림 6에 설명되어 있습니다.

Fig. 6

The five signaling pathways associated with macrophage polarization which are highlighted by pink color. The connection of each signaling pathway to the polarization to either M1 phase or M2 phase is denoted by brown and blue arrows respectively. The black arrows connect the signaling pathway to its key signaling protein(s) which are then connected to each activated phase. The complex interrelationship between signaling pathway and regulatory network is shown by connecting the pathway highlighted in yellow to the signaling pathway in pink by green arrows.

분홍색으로 강조 표시된 대식세포 분극화와 관련된 5개의 신호전달 경로입니다. 각 신호전달 경로가 M1 단계 또는 M2 단계로의 분극화와 연결되는 관계는 각각 갈색과 파란색 화살표로 표시되었습니다. 검은색 화살표는 신호전달 경로를 해당 경로의 핵심 신호전달 단백질(들)에 연결하며, 이 단백질들은 다시 활성화된 단계와 연결됩니다. 신호전달 경로와 조절 네트워크 간의 복잡한 상호작용은 노란색으로 강조 표시된 경로를 분홍색 신호전달 경로와 녹색 화살표로 연결하여 보여집니다.

Full size image

The MetaShiftOptimizer result indicates that the transition from the anti-inflammatory (M2) to the pro-inflammatory (M1) phase could be achieved by altering 30 reactions. For example, the transport of glycogenin G4G4, the transport of 5-hydroperoxyeicosatetraenoic acid (5(S)-HPETE) to mitochondria followed by the conversion of 5(S)-HETE to leukotriene B4 in arachidonic metabolism, the movement of L-carnitine from cytoplasm to mitochondria in the carnitine shuttle, and reactions in cholesterol metabolism58,59,60. Among these 30 reactions mentioned above, six reactions were to be upregulated, including the transport of glycogenin G4G4 and the transport of cofactors such as dTTP and dCTP to mitochondria. Two reactions (namely the transport of pitavastatin and exchange of potassium) were recommended for deletion, while the rest of the reactions were to be suppressed. By incorporating the suggested alterations, we achieved a modified version of M2 which shows flux modulation shifted more towards M1 state. On the other hand, the transition from the M1 to M2 state requires adjustments in a different set of 30 reactions across various pathways. These include the conversion of fructose-1,6-bisphosphate to fructose-6-phosphate and 2-phospho-D-glycerate to 3-phospho-D-glycerate—key initial reactions in glycolysis. Other necessary reactions involve the conversion of arachidonate to 5(S)-HPETE, which is subsequently metabolized to Leukotriene A4 (a precursor for LTB4 and LTC4) in arachidonic metabolism, the transformation of UDP-glucose to glycogen, and the production of adenosine60. In this scenario, eight reactions such as transport of 3-hydroxyisobutyrate from mitochondria to cytoplasm, the movement of succinate to cytoplasm, and the conversion of cyclic-3-hydroxymelatonin to hydroxide were suggested for upregulation, while one reaction that governs the movement of gliclazide was to be deleted. The rest of the reactions were needed to be downregulated. Further details regarding the reactions identified by the algorithm are available in Supplementary Data 5.

Following the adjustments made for the shift of M1 to M2 and vice versa, quantifying the difference in the metabolic flux is crucial to assess the impact of modifications on the overall metabolic behavior of the system. We used methods such as Hausdroff distance, Jaccard Distance, and Jaccard Index to measure the similarity/dissimilarity between M1 and M2 before and after adjustments (as described in Methods)61,62. Among the three methods mentioned, Hausdroff distance is a more reliable measure of dissimilarity in GSM models due to its interval-based representation (hence, can be used for FVA-predicted flux ranges), ability to incorporate uncertainty, and its sensitivity to the extreme values. The Hausdroff distance between the flux intervals of M1 and M2 is found to be 3.6789 e05 which decreases to 3.5018 e04 between normal M1 and modified M2. The distance between the normal M2 and the modified M1 did not show any significant change (both having a value of 3.6775 e05). This prompted us to further explore if the change is more prominent for select pathways essential for each activated phase, as reported in literature and this study, such as glycolysis/gluconeogenesis, OXPHOS, PPP, TCA cycle, OXPHOS, urea cycle, bile acid synthesis, pyruvate metabolism, and cholesterol metabolism were examined61,62. The Hausdroff distance for these select pathways before modification between M1 and M2 was found to be 4.30 e04 which decreases to 4.507 e3 after modification in M2. This clearly signifies that the trend of reduced distance between normal M1 and modified M2 holds true for these select pathways. However, the distance measure between normal M2 and modified M1 did not show any prominent change even for these select pathways. This resistance of M1 phase to transition to M2 phase could be due to the absence of crucial M2 reactions in M1. Study by Gelbach et al23 also mentioned that the activity of pathways such as chondroitin/heparan sulphate degradation and nucleotide sugar metabolism is unique to M2. Hence, we added 130 unique reactions present only in M2 in addition to the modification of 30 reactions as mentioned earlier and generated a new modified M1 model. We next calculated Hausdroff distance for the select pathways which was found to decrease from to 4.3037 e04 to 4.3009 e03 between modified M1 and normal M2. The details on the calculation of Hausdroff distance and other distance measures considered are available in Supplementary Data 5. The visual representation of flux distribution if the select pathways mentioned above, before and after modifications is shown using t-SNE plots in Fig. 7.

MetaShiftOptimizer 결과는

항염증성(M2) 단계에서 염증성(M1) 단계로의 전환이

30개의 반응을 변경함으로써 달성될 수 있음을 나타냅니다.

예를 들어,

글리코겐의 G4G4 운반,

5-하이드로페록시에이코사테트라엔산(5(S)-HPETE) 미토콘드리아로의 수송 후

아라키돈산 대사에서 5(S)-HETE의 류코트리엔 B4로의 전환,

카르니틴 셔틀에서 세포질에서 미토콘드리아로의 L-카르니틴 이동,

콜레스테롤 대사 반응58,59,60 등이 포함됩니다.

위에서 언급된 30개의 반응 중 6개 반응이 상향 조절되어야 하며, 이는 글리코겐인 G4G4의 운반 및 dTTP와 dCTP와 같은 보조인자의 미토콘드리아로의 운반을 포함합니다. 두 개의 반응(즉, 피타바스타틴의 운반 및 칼륨 교환)은 삭제 권장되었으며, 나머지 반응은 억제되어야 합니다. 제안된 변경 사항을 반영하여 M2의 수정된 버전을 얻었으며, 이는 유동 조절이 M1 상태로 더 이동된 특성을 보입니다. 반면, M1 상태에서 M2 상태로의 전환에는 다양한 경로에 걸쳐 다른 30개의 반응에 대한 조정이 필요합니다. 이에는 글리코lysis의 초기 반응인 프럭토스-1,6-비스포스페이트에서 프럭토스-6-포스페이트로의 전환과 2-포스포-D-글리세레이트에서 3-포스포-D-글리세레이트로의 전환이 포함됩니다. 기타 필수 반응에는 아라키돈산이 5(S)-HPETE로 전환되는 과정이 포함되며, 이는 아라키돈산 대사 과정에서 LTB4와 LTC4의 전구체인 류코트리엔 A4로 대사됩니다. 또한 UDP-글루코스가 글리코겐으로 전환되는 과정과 아데노신60의 생성도 포함됩니다. 이 시나리오에서, 미토콘드리아에서 세포질로 3-하이드록시이소부티레이트의 운반, 수크신산의 세포질로의 이동, 순환형 3-하이드록시멜라토닌의 수산화물로 전환과 같은 8개의 반응이 상향 조절될 것으로 제안되었으며, 글리클라지드의 이동을 조절하는 1개의 반응은 삭제될 것으로 제안되었습니다. 나머지 반응들은 하향 조절될 필요가 있었습니다. 알고리즘으로 식별된 반응에 대한 추가 세부 사항은 보충 자료 5에 수록되어 있습니다.

M1에서 M2로, 또는 그 반대로의 전환에 따른 조정을 반영한 후, 대사 유동량의 차이를 정량화하는 것은 시스템의 전체 대사 행동에 대한 수정 사항의 영향을 평가하는 데 중요합니다. 우리는 Hausdroff 거리, Jaccard 거리, Jaccard 지수 등의 방법을 사용하여 조정 전후의 M1과 M2 간의 유사성/불일치도를 측정했습니다(방법 섹션에 설명됨)61,62. 세 가지 방법 중 Hausdroff 거리는 GSM 모델에서 불일치도를 측정하는 더 신뢰할 수 있는 방법입니다. 이는 간격 기반 표현(따라서 FVA 예측 유동 범위에도 적용 가능), 불확실성 반영 능력, 극값에 대한 민감도 때문입니다. M1과 M2의 유동 간격 사이의 Hausdroff 거리는 3.6789 e05로 측정되었으며, 정상 M1과 수정된 M2 사이에서는 3.5018 e04로 감소했습니다. 정상 M2와 수정된 M1 사이의 거리는 유의미한 변화를 보이지 않았습니다(두 값 모두 3.6775 e05). 이 결과는 문헌 및 본 연구에서 보고된 각 활성화 단계에 필수적인 선택된 경로(글리코lysis/글루코네오게네시스, OXPHOS, PPP, TCA 사이클, OXPHOS, 요소 순환, 담즙산 합성, 피루vate 대사, 콜레스테롤 대사 등)에서 변화가 더 두드러지는지 추가로 탐구하도록 유도했습니다.61,62. M1과 M2 사이의 수정 전 선택적 경로의 Hausdorff 거리는 4.30 e04였으며, M2에서의 수정 후 4.507 e3로 감소했습니다. 이는 정상 M1과 수정된 M2 사이의 거리 감소 추세가 이 선택적 경로에서도 일관되게 나타남을 명확히 보여줍니다. 그러나 정상 M2와 수정된 M1 사이의 거리 측정치는 이 선택적 경로에서도 유의미한 변화를 보이지 않았습니다. M1 단계가 M2 단계로 전환되는 데 대한 저항성은 M1에 중요한 M2 반응이 결여되어 있기 때문일 수 있습니다. Gelbach 등23의 연구에서도 콘드로이틴/헤파란 설페이트 분해 및 뉴클레오티드 당 대사 경로의 활성이 M2에 특유하다는 점이 언급되었습니다. 따라서, 이전에 언급된 30개의 반응 수정 외에도 M2에만 존재하는 130개의 고유 반응을 추가하여 새로운 수정된 M1 모델을 생성했습니다. 이후 선택된 경로에 대한 Hausdorff 거리를 계산한 결과, 수정된 M1과 정상 M2 사이에서 4.3037 × 10⁻⁴에서 4.3009 × 10⁻³로 감소했습니다. Hausdroff 거리 계산 및 기타 거리 측정 방법에 대한 세부 사항은 보충 자료 5에 수록되어 있습니다. 위에서 언급된 선택된 경로의 유동 분포는 수정 전후를 t-SNE 플롯으로 시각화하여 그림 7에 표시되었습니다.

Fig. 7

t-SNE plot visualizing the M1 phenotype, M2 phenotype, and the modified M2 phenotype represented by blue, red, and green, respectively. (a) The flux distribution of M1 and M2 polarized state before any modifications (green represents the M2 and yellow represents M1 flux). (b) The t-SNE plot of Modified M2 (adjustment of 30 reactions) represented by red. (c) Understanding the metabolic shift in M1 select pathways with modifications on 30 reactions. (d) The new flux distribution after adding 130 unique M2 reactions that shift M1 to M2 more effectively. Modified M1 in both (c) and (d) is marked in blue.

Full size image

Discussion

Advancing respiratory medicine, improving disease management, and developing targeted therapies necessitate a deep understanding of alveolar macrophage (AM) polarization mechanisms54. In this study, we reconstructed genome-scale metabolic (GSM) models for healthy and activated states of AM to investigate the metabolic shifts that occur when AMs adopt either a pro-inflammatory (M1) phenotype or an anti-inflammatory (M2) phenotype.

Among the integration techniques applied, E-flux effectively captured the biological intricacies of all activation states of alveolar macrophages (AM), generating more comprehensive metabolic models. For instance, genome-scale metabolic (GSM) models derived through E-flux successfully incorporated key AM-related pathways, including nitric oxide (NO) production, succinate and itaconate synthesis, and complete fatty acid oxidation (FAO) and fatty acid synthesis (FAS) pathways30. Thus, E-flux preserves the entire metabolic framework, enabling a holistic exploration of metabolism and providing a more in-depth analysis of metabolic reprogramming compared to models generated using switch-based methods20,21,22,23,24,29.

E-flux integrates gene expression‎ data with gene-protein-reaction (GPR) associations to impose constraints on metabolic reactions, making it a robust tool for constructing biologically relevant metabolic networks when suitable transcriptomic datasets are available. However, due to limitations in dataset availability, the closest approximation of the M2 macrophage phenotype in the lung was derived from bronchial tissue biopsies (GSM41084)26, which served as a proxy for lung-specific data. Despite multiple studies reporting that this dataset as a good representation of M2 macrophage in lungs, an ideal case will be to get iL-4 stimulated AM. To the best of our knowledge, at the time of reconstruction of these models, no such datasets were available. Recent advances such as study by Travaglini et al. (2020)63and Sun et al. (2021)56 have led to scRNA and snRNA sequencing of AM to study certain diseased states; however, such advanced data are not available for the M1 and M2 phase specifically. Therefore, the most suitable transcriptomics data available for all the conditions were selected.

One of the most important aspects of GSM models is determining the objective function that does justice to the biology of the cell. In most cases the biomass equation that defines the growth of the relevant cellular system is used as an objective function. However, an immune cell’s metabolic objective cannot be confined to a growth function (especially not for a tissue resident AM) and/or a specific metabolic function such as ATP or protein synthesis. Hence in this study, we used objective-independent techniques such as FVA, FSA and metabolite centric analysis along with a basal level of biomass as a constraint, thus ensuring the maintenance to biological activities. Next, we conducted exhaustive analysis of the metabolic pathways of AM and its phenotypes. In doing so, our GSM models highlighted the pathways and reactions that can be linked to the inflammatory/anti-inflammatory responses. These pathways, specifically bile acid synthesis, arachidonic acid metabolism, chondroitin heparan sulphate biosynthesis can be linked to the M1 phenotype acquirement and the initiation of inflammation. Chen et al. (2016) investigated three bile acids and the effect of change in their concentration to cytotoxicity in the cell44. Additionally, oxylipins which are obtained via arachidonic acid metabolism play a very important role in the regulation of inflammation and the formation of other important leukotriene metabolites such as LTA459,60.

Our results add to the preexisting potential of chondroitin/heparan sulphate biosynthesis and degradation as a significant modulator for initiation and resolution of inflammation47,48,49. Chondroitin/heparan sulphate biosynthesis and degradation pathways are known to either contribute to the formation or degradation of an important metabolite called Heparan sulphate47. The formation of heparan sulfate is a crucial step for the recruitment, adhesion, crawling, and transmigration of leukocytes from the circulation to the site of inflammation48. Despite their important roles, these pathways have not been fully explored for their potential contribution in lung injury and inflammation. The enhanced HS biosynthesis in M1 and enhanced HS degradation in M2 phenotype as predicted by the models show the potential to be key regulator of initiation and resolution of inflammation in AMs.

Our bilevel optimization framework MetaShiftOptimizer is able to identify set of reactions originating from several metabolic pathways such as production and transport of glycogenin and production of Leukotriene B4, that can be linked to various important signaling network such as the PI3K/AKT, TLRs/NF-kB, JAK/STAT, Notch, and JNK signaling that are known to drive the macrophage polarization55,57. Figure 6shows such metabolic reactions/pathways identified by MetaShiftOptimizer and the possible link to the signaling mechanisms. While several optimization frameworks such as OptForce, and MOMA ae extensively used to identify the reactions to enable a metabolic network to enhance a certain metabolite production, MetaShiftOptimizer identifies reactions, which when modified, can shift one phenotype to the other64,65. To this end, the decrease in distance between the flux spaces of normal and modified conditions show that the modifications were able to shift the phenotype to the desired state. Even the small differences in fluxes can be statistically and biologically significant due to cascading effects on downstream pathways. For example, a slight increase in a bottleneck reaction flux could indicate metabolic shifts that are biologically meaningful, particularly in tightly regulated systems or pathways involved in complex immune cells such as AM. Thus, this framework can be used to identify targets points to shift polarization to desired phenotype to effectively fight pathogens and diseases.

In conclusion, our models are able to successfully captivate the biological essence of the healthy and activated state of AM. We have used available experimental data and phenomenon reported in literature to validate our models. Our results indicate several potential targets such as pyruvate, glycogen, chondroitin/heparan sulphate biosynthesis and degradation, bile acid synthesis, arachidonic acid metabolism as driving forces in macrophage polarization and/or significant markers of each activated phase. However, the activated AM models in this study are simplified version unlike the complex M1/M2 dynamics as typically observed in vivo. In reality, AM polarization is highly plastic, and context dependent and influenced by a complex array of signals in the lung microenvironment that results in a wide range of functional states beyond the M1/M2 dichotomy which unfortunately is very challenging to capture by any GSM model due to lack of suitable omics data. Additionally, our study is also limited to the in silico analysis and does not directly include experimental validation while substantiating the observed phenomena based on literature evidence. In future, with further advancement in study of AM and its activated phenotypes and when advanced omics data become available, these models will further be tested and refined for better prediction capability. Moving forward, experimental validation would be our focus (through establishing collaboration) so that some of the model findings can be tested.

논의

호흡기 의학의 발전, 질병 관리의 개선,

그리고 표적 치료법의 개발은

폐포 대식세포(AM)의 분화 메커니즘에 대한 깊은 이해를 요구합니다.

본 연구에서는 건강한 상태와 활성화된 상태의 AM에 대한 유전체 규모 대사(GSM) 모델을 재구성하여

AM이 염증 촉진형(M1) 또는 염증 억제형(M2) 표현형을 채택할 때 발생하는 대사 변화를 조사했습니다.

적용된 통합 기술 중 E-flux는 폐포 대식세포(AM)의 모든 활성화 상태의 생물학적 복잡성을 효과적으로 포착하여 더 포괄적인 대사 모델을 생성했습니다.

예를 들어, E-flux를 통해 도출된 유전체 규모 대사(GSM) 모델은 질산산화물(NO) 생성, 수산화물 및 이타콘산 합성, 완전 지방산 산화(FAO) 및 지방산 합성(FAS) 경로 등 AM과 관련된 주요 경로를 성공적으로 포함했습니다30. 따라서 E-flux는 전체 대사 프레임워크를 유지하여 대사 과정을 종합적으로 탐구할 수 있으며, 스위치 기반 방법20,21,22,23,24,29를 사용해 생성된 모델보다 대사 재프로그래밍에 대한 더 깊은 분석을 제공합니다.

E-flux는 유전자 발현 데이터와 유전자-단백질-반응(GPR) 연관성을 통합하여 대사 반응에 제약을 가함으로써, 적절한 전사체 데이터셋이 사용할 수 있을 때 생물학적으로 관련성 있는 대사 네트워크를 구축하는 강력한 도구입니다. 그러나 데이터셋의 제한으로 인해 폐의 M2 대식세포 표현형을 가장 잘 근사한 것은 기관지 조직 생검(GSM41084)26에서 유래했으며, 이는 폐 특이적 데이터의 대용으로 사용되었습니다. 이 데이터셋이 폐 내 M2 대식세포의 좋은 대표라고 보고된 다수의 연구가 있음에도 불구하고, 이상적인 경우는 iL-4 자극된 AM을 얻는 것입니다. 본 모델 재구성 시점까지 이러한 데이터셋은 존재하지 않았습니다. Travaglini 등(2020)63과 Sun 등(2021)56의 최근 연구는 특정 질환 상태를 연구하기 위해 AM의 단일 세포 RNA(scRNA) 및 단일 핵 RNA(snRNA) 시퀀싱을 수행했지만, M1 및 M2 단계에 특화된 이러한 고급 데이터는 없습니다. 따라서 모든 조건에 적합한 가장 적절한 전사체학 데이터가 선택되었습니다.

GSM 모델의 가장 중요한 측면 중 하나는 세포의 생물학을 적절히 반영하는 목적 함수를 결정하는 것입니다. 대부분의 경우 관련 세포 시스템의 성장을 정의하는 생물량 방정식이 목적 함수로 사용됩니다. 그러나 면역 세포의 대사적 목표는 성장 기능(특히 조직 거주형 AM의 경우)이나 ATP 또는 단백질 합성과 같은 특정 대사 기능으로 제한될 수 없습니다. 따라서 본 연구에서는 FVA, FSA 및 대사체 중심 분석과 같은 목표 독립적 기술을 사용했으며, 생물학적 활동을 유지하기 위해 기초 수준의 생물량을 제약 조건으로 설정했습니다. 다음으로, AM의 대사 경로와 그 표현형을 대상으로 포괄적인 분석을 수행했습니다. 이 과정에서 GSM 모델은 염증/항염증 반응과 연관된 경로와 반응을 강조했습니다. 특히 담즙산 합성, 아라키돈산 대사, 콘드로이틴 헤파란 설페이트 생합성 경로는 M1 표현형 획득과 염증 발현의 시작과 연관될 수 있습니다. Chen et al. (2016)은 세 가지 담즙산과 그 농도 변화가 세포 독성에 미치는 영향을 조사했습니다44. 또한 아라키돈산 대사로부터 생성되는 옥실리핀은 염증 조절과 LTA4와 같은 중요한 류코트리엔 대사물의 형성에 매우 중요한 역할을 합니다59,60.

우리의 결과는 콘드로이틴/헤파란 황산염 생합성 및 분해가 염증의 발현과 해결에 중요한 조절인자로서의 잠재력을 추가로 입증합니다47,48,49. 콘드로이틴/헤파란 황산염 생합성 및 분해 경로는 중요한 대사산물인 헤파란 황산염의 형성 또는 분해에 기여하는 것으로 알려져 있습니다47. 헤파란 황산염의 형성은 순환계에서 염증 부위로 백혈구의 모집, 부착, 이동, 및 전이 과정에 필수적인 단계입니다48. 이러한 경로는 중요한 역할을 함에도 불구하고 폐 손상 및 염증에서의 잠재적 기여도는 충분히 탐구되지 않았습니다. 모델 예측에 따라 M1 형질에서 HS 생합성 증가와 M2 형질에서 HS 분해 증가가 관찰된 것은 AM에서 염증의 발현과 해결을 조절하는 핵심 조절자로서의 잠재성을 보여줍니다.

우리의 이중 최적화 프레임워크 MetaShiftOptimizer는 글리코겐 생산 및 운반과 류코트리엔 B4 생산과 같은 여러 대사 경로에서 유래한 반응 집합을 식별할 수 있으며, 이는 PI3K/AKT, TLRs/NF-kB, JAK/STAT, Notch, JNK 신호전달 경로와 같은 중요한 신호전달 네트워크와 연결될 수 있습니다. 이러한 경로는 대식세포 분극화를 유도하는 것으로 알려져 있습니다55,57. 그림 6은 MetaShiftOptimizer가 식별한 대사 반응/경로와 신호전달 메커니즘과의 가능한 연결 관계를 보여줍니다. OptForce 및 MOMA와 같은 여러 최적화 프레임워크는 특정 대사체 생산을 향상시키기 위해 대사 네트워크를 활성화하기 위한 반응을 식별하는 데 널리 사용됩니다. 반면 MetaShiftOptimizer는 반응을 식별하여 이를 수정함으로써 한 표현형을 다른 표현형으로 전환할 수 있는 반응을 찾아냅니다64,65. 이 목표를 달성하기 위해 정상 조건과 수정된 조건 간의 유동 공간 간의 거리 감소는 수정된 조건이 표현형을 원하는 상태로 전환할 수 있음을 보여줍니다. 유동량의 작은 차이도 하류 경로에 대한 연쇄 효과로 인해 통계적 및 생물학적 의미가 있을 수 있습니다. 예를 들어, 병목 반응의 유동량 증가가 생물학적으로 의미 있는 대사 변화, 특히 AM과 같은 복잡한 면역 세포에 관여하는 엄격히 조절되는 시스템이나 경로에서 발생할 수 있습니다. 따라서 이 프레임워크는 병리체 및 질병과 효과적으로 싸우기 위해 극성 전환을 원하는 표현형으로 전환하기 위한 목표 지점을 식별하는 데 사용될 수 있습니다.

결론적으로, 우리 모델은 AM의 건강한 상태와 활성화된 상태의 생물학적 본질을 성공적으로 포착했습니다. 우리는 기존 실험 데이터와 문헌에 보고된 현상을 활용해 모델을 검증했습니다. 결과는 피루vate, 글리코겐, 콘드로이틴/헤파린 황산염 생합성 및 분해, 담즙산 합성, 아라키돈산 대사 등 대식세포 극성화의 주요 동인 또는 각 활성화 단계의 중요한 지표로 작용할 수 있는 여러 잠재적 표적을 제시했습니다. 그러나 본 연구의 활성화된 AM 모델은 체내에서 일반적으로 관찰되는 복잡한 M1/M2 동역학과 달리 단순화된 버전입니다. 실제로 AM 분극화는 매우 가소성이 높으며, 맥락에 의존적이며 폐 미세환경의 복잡한 신호 체계에 의해 영향을 받아 M1/M2 이분법 beyond를 넘어 다양한 기능적 상태를 보입니다. 이는 적합한 오믹스 데이터의 부족으로 인해 어떤 GSM 모델로도 포착하기 매우 어려운 특징입니다. 또한 본 연구는 in silico 분석에 한정되어 있으며, 문헌 증거를 바탕으로 관찰된 현상을 뒷받침하기 위해 실험적 검증은 직접 포함하지 않았습니다. 향후 AM 및 그 활성화된 표현형에 대한 연구가 진전되고 고급 오믹스 데이터가 확보되면, 이러한 모델은 예측 능력을 향상시키기 위해 추가로 테스트 및 정교화될 것입니다. 앞으로 실험적 검증(협력을 통해)이 우리의 주요 초점이 될 것이며, 이를 통해 모델의 일부 결과를 검증할 수 있을 것입니다.

Methods

Transcriptomics data processing

An exhaustive literature search was conducted to identify the appropriate set of transcriptomics data which included the transcriptomic profiles of healthy non-smokers24AM induced by Lipopolysaccharides (LPS) and bronchial tissues induced by interleukin-4 (iL-4) resembling M1 phase25and M2 phase26, respectively. The data obtained were used as input for Gene Set Enrichment Analysis (GSEA) tool28. GSEA is a tool that is used for pathway analysis based on the transcriptomic state of the cells. It was used to compare the pathway activity of healthy AM with the M1 phase, healthy AM with the M2 phase, and the M1 with the M2 phase. In this process, genes are ranked based on the correlation between their expression‎ and the class distinction using any suitable metric. GSEA calculates the enrichment score (ES) and its significance level using p-values28. The default metrices such as meandiv was used in the desktop GSEA for this study. The output from the GSEA run generated lists of enriched pathways for the M1 and M2 phase that mainly focused on signaling pathways and are available in the GitHub repository. The obtained results are impar with the evidence reported in literature for M1/M2 and healthy phenotype of alveolar macrophages.

Using the raw data set, we deduced a list of genes that were also present in the Human1 metabolic model. Quantile normalization technique was used for all the datasets to obtain the gene expression‎ values at all conditions. After normalization, the average was taken for each condition from all the samples and then was scaled to −1000to 1000. The list of genes for healthy AM, M1 phase, and M2 phase were 2,173, 2,951, and 2,390 respectively. The expression‎ values of these genes were integrated into Human123model by using different integration techniques (iMAT21, mCADRE20, INIT29and E-flux22) as described below, to reconstruct models of healthy AM and its activated phases. The details on the datasets and normalization are available in Supplementary Data 6.

GSM model reconstruction

The transcriptomics data obtained for each of the phenotypes of AM was integrated into Human1, a global human metabolic reconstruction consisting of 13,417 reactions, 10,138 metabolites (4,164 unique), and 3625 genes23. Three context-specific AM metabolic reconstructions were obtained by implementing both switch and valve approaches of omics integration. Among various methods available in both the categories of switch and valve approach, iMAT21, INIT29, mCADRE20, and E-flux22were used in our study. iMAT (integrative metabolic analysis tool) is an optimization-based program that can be used to integrate the available omics data with GSM models for the prediction of metabolic fluxes20. The modified version of iMAT was used in which instead of classifying the overall reactions into three categories (highly expressed, lowly expressed, and moderately expressed), the reactions were divided as either highly expressed or lowly expressed with the biomass precursors always included in the highly expressed set. The cutoff point of a p-value less than 0.05 was used to include reactions in the highly expressed set (the same threshold was used for all conditions). The biomass reaction for the healthy and the activated state was adopted from Bordbar et al16. The formulation was constructed in such a way that all the reactions from the highly expressed set were always made active and the minimum number from the lowly expressed reaction set was added to obtain the specified objective. This resulted in a pruned model which is significantly smaller than the original human model with the reactions, and metabolites specific to AM and its activated stages. On the other hand, E-flux only requires the change of the upper bound and lower bound on each reaction depending on the gene expression‎ level22. The forward reactions consisted of a lower bound of 0 and a unique upper bound according to the gene expression‎ levels. The backward reactions ranged from a unique lower limit to 0 as an upper bound. And for the reversible reactions, the range spans from the negative value specific to each reaction to its corresponding positive value based on the gene expression‎ levels. Additionally, GSM models via INIT (Integrative Network Inference for Tissues) and mCADRE (Context-specificity Assessed by Deterministic Reaction Eval‎uation) methods were reconstructed20,21,22,29. Both the approaches integrate the expression‎ data into the global Human metabolic model while ensuring the production of certain important metabolites. However, INIT, mCADRE, and iMAT prune the human model specific to a context by turning off a number of reactions. This could lead to the omission of certain reactions that are crucial for the specific phase (healthy or M1 or M2), while E-flux does not deactivate any reaction. By not pruning reactions based on gene expression‎, E-flux provides a more inclusive representation of the metabolic models, allowing for a broader exploration of pathway activities. This approach is particularly useful to capture the overall functionality of the metabolic models under different conditions without biasing these towards specific gene expression‎ patterns.

Hence, the GSM models for healthy AM, M1, and M2 phases were obtained by implication of the mentioned approaches. We compared the distribution of metabolites, reactions, and pathways in all the GSM models and the details are discussed in the Results and Discussion section. To ensure the biological relevance of these GSM models, we used techniques such as Flux Balance Analysis (FBA)33and Flux Variability Analysis (FVA)34 to analyze and improve model connectivity.

Flux balance, flux sum, flux variability, metabolite centric analysis

Flux Balance Analysis (FBA) is used in this study to analyze the flow of metabolites in different conditions. FBA is a widely used approach to study biochemical networks, namely, genome-scale metabolic models that contain the known metabolic reactions in a biological system and the associated genes and enzymes66. The GSM model is represented by a stoichiometric matrix which contains the metabolites as columns and the reactions as rows. The upper and lower bound act as a constraint on each of the reactions based on nutrient availability and other (microenvironment or genetic) conditions. FBA generates a flux value for each reaction. Flux Variability Analysis (FVA) is an extension of FBA which calculates the maximum and minimum possible flux for each of the reactions in the model at a specific condition33. Additionally, Flux Sum Analysis39 is used to obtain the metabolite pool size in different conditions (M1 and M2). FSA is conducted for important metabolites such as ATP, NADPH, NO, succinate, and itaconate in normal (without turning off any pathways), glycolysis off and OXPHOS/TCA off condition. The metabolite centric analysis was conducted by using RPAm function in COBRA toolbox in MATLAB (the details are available in the Supplementary Data 4).

Model curation

The three metabolic models were curated by using the classic design-build-test-refine cycle to be able to ensure proper network connectivity and accurate reflection of the metabolic capabilities of the AM cell. To ensure all the important AM pathways such as NO cycle, FAS, FAO, were activated, thermodynamically infeasible cycles (TICs) need to be properly handled. TICs are cycles created by reactions that carry fluxes even in the absence of nutrients essential for cellular growth and functionality. The TICs can cause the metabolic model to produce metabolites higher/lower than expected, by activating reactions that would be off in a biological scenario. However, if essential reactions are eliminated or the directionality of these reactions are changed without proper review, the behavior of the metabolic model might shift away from the known biological phenomena of the cell. Hence, it is extremely important to refine metabolic models by using efficient and effective methods.

We used OptRecon (an inhouse tool, currently unpublished), that has been developed as an expansion of OptFill​​, a tool previously developed by our group with different functionalities67. The initial function of OptRecon is to refine GSMs by removing TICs; however, the process of removing TICs from GSMs was found to be much more difficult than the process of incorporating reactions without creating TICs, and thus the method was upgraded to be able to expand a minimal model (i.e., minimum number of biochemical reactions required to satisfy the objective, in our models the number was found to be 143) by adding reactions from a database (the database consisted of all but these 143 reactions from Human1). OptRecon generated three possible solutions to avoid formation of any TICs and ensure optimal connectivity. These solutions consisted of either blocking a reaction completely or changing the direction of the reaction. Before incorporating any changes, an exhaustive literature search was conducted to ensure that none of the biologically relevant pathways were fully omitted or partially affected due to these changes. Special attention was given to novel AM pathways such as production of NO from arginine in healthy AM, production of succinate, itaconate and citrate in TCA cycle during M1 phase and citrulline and urea production in NO cycle during M2 phase. FBA and FVA techniques were used to check the fluxes of the metabolic models ensuring proper network connectivity. All the fluxes from FBA and FVA in the absence of nutrients were found to be zero as expected in healthy and activated AM GSM models, while in other conditions the fluxes were found to be in accordance with the biological nature of AM.

Identification of minimal modification for polarization shift

The GSM model for each activated phase comprises of numerous reactions. When comparing the flux range distribution across the M1 and M2 phases, we examined whether any reactions exhibited a unique distribution. However, manually observing and comparing each reaction, which may have different flux ranges and potentially contribute to polarization shifting when limited, was a highly time-consuming approach to address the problem.Hence, in this work we propose a bilevel optimization formulation which identifies a set of reactions capable to shift the one phenotype to another. Previous works such as OptForce ensures the overproduction of desired metabolite while considering the flux values of wildtype strain by identifying a force set65,67. Here, our focus lies in identifying the specific reactions that, when manipulated through upregulation, downregulation, or deletion, can induce a shift in the metabolic network from one activated state (M1 or M2) to the other. To achieve this, the minimization of metabolic interventions was optimized while ensuring the minimal discrepancy in reaction fluxes between the M1 and M2 states.

Initially, we began by creating a collection of reactions that fall into three categories: those that must be upregulated (MUSTU), those that must be downregulated (MUSTL), and those that must be deleted (MUSTdel). In order to distinguish these three sets of reactions, we conducted a comparison of the flux ranges acquired via Flux Variability Analysis (FVA). These three sets have unique reactions i.e. a reaction included in a MUSTU set only belongs to MUSTU and is not present in other sets. The minimal alterations required for the transition of both M1 and M2 to the alternate phase encompassed 30 reactions, originating from diverse pathways, and metabolizing distinct compounds.

The formulation is following:

min∑jyu+yl+yx

Subject to: min∑j|M1(j)−M2(j)|

Subject to:∑jSij.vj=0∀i∈I

(1)

vj≥vjmax.yu+LBj(1−yu)∀j∈MUSTU

(2)

vj≤vjmin+UBj(1−yL)∀j∈MUSTL

(3)

LBj(1−yx)≤vj≤+UBj(1−yx)∀j∈MUSTdel

(4)

LBj≤vj≤UBj∀j∈J

(5)

vbiomass≥0.03

(6)

In the above formulation, the yu, yL and yx are the binary variables associated with the reactions in the MUSTU, MUSTL and MUSTdel sets. The binary variable is 1 if the reaction is to be modified (For example yu is 1 when the reaction is to be upregulated and 0 otherwise). M1(j) and M2(j) are the sets of all the reaction fluxes in the M1 GSM model and the M2 GSM model, respectively, under similar nutrient conditions (obtained from FBA). Sij is the stoichiometric coefficient for the metabolite i in the reaction j and vj is the flux of each reaction. Sets I and J include all metabolites and reactions known to occur within AM and its activated phenotypes, respectively. Each reaction is constrained by a upper bound (UBj) and a lower bound (LBj) that specify the direction based on thermodynamic constraints. The pseudo-steady state assumed in Eq. 1 allows the assumption that all internal metabolites are consumed and produced at an equal rate. The Eqs. 2, 3, and 4 re-adjust the upper bound and lower bound of each of the reactions that is to be upregulated, downregulated and deleted accordingly. Finally, as shown in Eq. 6, in order to ensure a minimum level of biomass can be produced, a flux of 0.03 mol/g cell DW/h is forced, which corresponds to the biomass flux under healthy condition. From this algorithm, a set of reactions is identified which when modified as suggested shift the polarized state of the activated AM. The reactions identified by the modified algorithm can be linked to the signaling pathways that are known to be directly responsible for macrophage polarization.

Measure of similarity/dissimilarity

In order to assess the overall impact and extent of change resulting from the introduction of constraints in the GSM models, we employed techniques such as Hausdorff distance, Jaccard distance, and Jaccard index61. The Jaccard index is computed by considering the intersection and union of a single data point. Flux Balance Analysis (FBA) is employed to derive a permissible distribution of metabolic fluxes in a steady-state system within a Genome-Scale Metabolic (GSM) model. However, it is important to note that the resulting fluxes are not exclusive solutions. The GSM models are typically characterized by being underdetermined, context-specific, and physiologically relevant flux solutions66. These solutions can be further refined to a single solution by imposing additional restrictions62. Eflux2 is an enhanced version of FBA that deduces a metabolic flux distribution from transcriptomics data and addresses the limitation of E-flux by offering a singular solution62. The Jaccard similarity index is a quantitative measure that assesses the degree of similarity between two sets of data. This index was computed using this distinct solution. It is expressed as a percentage range from 0 to 100%, with a larger percentage indicating a greater level of similarity61,68. However, the solution of E-flux2 is not constant and varies depending on the specified value of the goal function. For instance, the solution set derived from FBA’s maximum biomass differs from the solution set derived from FVA’s maximum flux set as biomass. To the best of our knowledge, there is no specific growth rate documented for alveolar macrophages. Therefore, the flux ranges obtained via Flux Variability Analysis (FVA) were utilized to produce 50,000 samples using the Flux Sampling method39 for both the regular activated state and the modified activated states. These samples were then used to calculate the Jaccard Index, Jaccard distance and Hausdorff distance, and to construct t-SNE plots. Supplementary Data 5 contains detailed information on the concept and calculation of the Jaccard index, Jaccard distance, and Hausdorff distance.

Data availability

Data Availability: All the generated and used data information is provided in this paper. Additionally, the GSM models (xml versions), and the GSEA results can be found in this GitHub repository: https://github.com/ssbio/Alveolar-Macrophage. Code Availability: The code for MetaShiftOptimizer in both GAMS and python scripts are also available in the same repository (https://github.com/ssbio/Alveolar-Macrophage).

Code availability

The code for MetaShiftOptimizer in both GAMS and python scripts are also available in the same repository (https://github.com/ssbio/Alveolar-Macrophage).

References

1. Blussé Van Oud Alblas, A. & Van Furth, R. The origin of pulmonary macrophages. Immunobiology 161, 186–192 (1982).

   Article PubMed Google Scholar 

2. Kulikauskaite, J. & Wack, A. Teaching old dogs new tricks? The plasticity of lung alveolar macrophage subsets. Trends Immunol. 41, 864–877 (2020).

   Article CAS PubMed PubMed Central Google Scholar 

3. Van Furth, R. & Blussé Van Oud Alblas, A. The current view on the origin of pulmonary macrophages. Pathol. Res. Pract. 175, 38–49 (1982).

   Article PubMed Google Scholar 

4. Green, G. M., Jakab, G. J., Low, R. B. & Davis, G. S. Defense mechanisms of the respiratory membrane. Am. Rev. Respir. Dis. 115, 479–514 (1977).

   CAS PubMed Google Scholar 

5. Fall, F. et al. Metabolic reprograming of LPS-stimulated human lung macrop