Large-scale computational modelling of the M1 and M2 synovial macrophages in rheumatoid arthritis

Rheumatoid arthritis (RA) is a chronic inflammatory autoimmune disease lacking a cure. Current treatments primarily manage symptoms, often with adverse effects. The innate immune system, particularly macrophages, plays a significant role in RA pathogenesis. Macrophages exist in various phenotypes, with M1 macrophages driving inflammation through cytokine production and matrix degradation, and M2 macrophages promoting resolution through anti-inflammatory cytokine production and tissue repair. In RA, an imbalance exists with a higher M1/M2 ratio, highlighting the potential of targeting macrophages therapeutically. While molecular interaction maps provide valuable information, their static nature limits predictive capabilities. Dynamic modeling, such as Boolean modeling, offers a solution to analyze emergent behaviors in complex systems. However, analyzing large-scale Boolean models is computationally demanding. This study aimed to develop an efficient computational framework to overcome this challenge, allowing for comprehensive analysis of large-scale Boolean models of M1 and M2 macrophages in RA.

The literature highlights the critical role of macrophages in RA pathogenesis, with an imbalance between pro-inflammatory (M1) and anti-inflammatory (M2) phenotypes contributing to disease severity. Existing approaches target macrophage modulation, but lack of specific clinical therapies targeting macrophages necessitates further investigation. High-throughput experimental techniques generate extensive data on molecular interactions, often represented as networks. While molecular interaction maps offer a rich knowledge source, their static nature limits their predictive power. Dynamic modeling, specifically qualitative approaches like Boolean models, is better suited for systems lacking detailed kinetic data and allows exploration of emergent properties. Previous work has demonstrated the application of Boolean modeling to biological systems, but scaling to large-scale models remains a challenge.

This study developed a computational framework to build, analyze, and validate large-scale Boolean models of RA macrophages. The framework utilizes publicly available molecular interaction maps from the RA-Atlas, which were converted into executable Boolean models using the CaSQ tool. A modified version of the BMA tool, deployed on a high-performance computing cluster, enabled efficient identification of model attractors (steady states). The identified steady states were validated against gene expression data (GSE97779 dataset) and existing literature. Differentially expressed genes (DEGs) between RA and healthy control synovial macrophages were identified using both literature review and transcriptomic data analysis. The expression levels were discretized into binary values (0 for underexpressed, 1 for overexpressed). Similarity scores were computed to compare model steady states with the experimentally observed DEG expression profiles. Steady states with the highest similarity scores were selected, and their average vector represented the calibrated model state. In silico knockout (KO) simulations were then performed to assess the effects of single and double knockouts of potential therapeutic targets on macrophage phenotypes (apoptosis and proliferation). The Therapeutic Target Database (TTD) was used to identify potential therapeutic targets already experimentally modulated, focusing on those with at least one inhibitor. The effects of these KO simulations were evaluated by comparing the perturbed model states to the calibrated states. Receptor double knockouts were also simulated to investigate their synergistic effects.

The framework successfully generated and calibrated Boolean models for both RA M1 and M2 macrophages. The RA M1 macrophage model comprised 233 nodes and 64 inputs, while the RA M2 macrophage model had 169 nodes and 39 inputs. Model validation showed high similarity between the model's predicted states and experimentally observed values (99% for M1, 96.5% for M2). In silico KO simulations identified several potential therapeutic targets: * **NFkB:** Inhibition led to M1 macrophage apoptosis and suppressed proliferation, demonstrating its potential as a therapeutic target. * **ERK1:** Inhibition suppressed M1 macrophage proliferation. * **GSK3B:** Inhibition promoted M2 macrophage proliferation and suppressed apoptosis. * **JAK1/JAK2:** Combined inhibition suppressed M1 macrophage proliferation and induced apoptosis. * **Notch1/ERK1:** Combined inhibition induced M1 macrophage apoptosis and suppressed proliferation. Double knockout simulations revealed synergistic effects for some target combinations, particularly in promoting M1 macrophage apoptosis. However, double knockouts of receptors did not significantly alter macrophage phenotypes, suggesting complex crosstalk between signaling pathways.

The study successfully demonstrated the feasibility of using a large-scale Boolean modeling approach to analyze complex biological systems, specifically macrophage behavior in RA. The high-performance computing framework allowed for efficient analysis of large models, overcoming computational limitations. The identification of NFkB, JAK1/JAK2, ERK1/Notch1, and GSK3B as promising therapeutic targets aligns with existing literature and provides potential new drug combinations. However, the limitations of the M1/M2 dichotomy as a simplified representation of macrophage activation and the lack of phenotype-specific data for model calibration need to be considered. Future work should incorporate a wider range of macrophage phenotypes and integrate additional datasets to improve model robustness and accuracy.

This study provides a robust computational framework for analyzing large-scale Boolean models of complex biological systems. The application to RA macrophages identified several potential therapeutic targets and drug combinations, offering insights into mechanisms underlying RA pathogenesis. Future research should focus on expanding the model to incorporate a wider range of macrophage phenotypes and integrating additional data for improved model accuracy and prediction of therapeutic efficacy.

The study's limitations include the simplification of macrophage activation into the M1/M2 dichotomy, which might not fully capture the complexity of macrophage heterogeneity. The model's calibration relied on available data, and the inclusion of additional, particularly phenotype-specific, data would improve model accuracy. The Boolean modeling approach, while efficient, is qualitative and may not fully capture the quantitative aspects of biological processes.