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

Rheumatoid arthritis is a chronic autoimmune inflammatory disease characterized by persistent synovial inflammation, joint destruction, and systemic comorbidities. Synovial macrophages are abundant and correlate with disease severity, with an imbalanced M1/M2 ratio skewed toward proinflammatory M1 cells. The study addresses the need for macrophage-focused therapeutic strategies by asking whether large-scale logic-based models can capture disease- and cell-type-specific macrophage phenotypes and predict perturbations that selectively suppress M1 while supporting M2 in RA. The purpose is to build executable Boolean models from curated RA macrophage interaction maps, validate them against transcriptomic data and prior knowledge, and use them to prioritize mono- and combination therapies that modulate apoptosis and proliferation phenotypes in M1 and M2 macrophages. This is important because existing RA therapies do not target macrophages specifically, and network-level, dynamic understanding may reveal selective intervention points.

The paper reviews systems biology approaches for complex diseases, contrasting quantitative models—which require kinetic parameters and are limited in scale—with qualitative methods scalable to hundreds of components. Petri nets can capture multiple abstractions but are challenging for large biological networks lacking stoichiometry. Logic-based models, especially Boolean models, are parameter-free and intuitive for signalling and gene regulatory networks, representing species as ON/OFF with dynamics defined by logical rules (AND/OR/NOT). Attractors (steady states or cycles) correspond to long-term phenotypes, making their computation central but challenging due to exponential state spaces. Prior work has produced curated, static molecular interaction maps for diseases (e.g., RA-Atlas) but their static nature limits prediction. Automatic inference tools (e.g., CaSQ) can convert maps into executable models, yet large input spaces and many inputs complicate attractor identification and biological contextualization.

- Maps to Boolean models: Used RA-Atlas CellDesigner SBML maps for RA M1 and M2 synovial macrophages and converted them to executable Boolean models with CaSQ v1.1.4, exporting BMA JSON format (granularity 1). Initial network sizes: M1 map 601 components/405 reactions → Boolean model 309 nodes, 75 inputs, 562 interactions; M2 map 513 components/323 reactions → 254 nodes, 57 inputs, 430 interactions. - Focus on cell-specific phenotypes: Selected nodes upstream of phenotypes to focus analysis. For M1: upstream of Apoptosis, Proliferation, Osteoclastogenesis → reduced to 233 nodes (64 inputs). For M2: upstream of Apoptosis and Proliferation → reduced to 169 nodes (39 inputs). Inputs to excluded nodes set to 1 (BMA default). - Differential expression evidence: Curated literature for RA synovial macrophage-specific expression changes and analyzed GEO microarray GSE97779 (9 RA synovial macrophage samples, 5 healthy PB monocyte-derived macrophage samples). Quantile normalization (preprocessCore), DEA with limma, adj. p<0.05. Identified differentially expressed molecules mapped to model nodes and discretized: overexpressed=1, underexpressed=0. Totals: 105 DE components in M1 model; 88 in M2 model. - Attractor computation: Fixed inputs corresponding to DE evidence to reduce combinations. M1: fixed 43 of 64 inputs → 2^21 input combinations. M2: fixed 24 of 39 inputs → 2^15 combinations. Used a Linux-deployable BMA console tool on an HPC node (96 single-core CPUs, 768 GB RAM) with parallelization (joblib) to compute attractors under synchronous update; all found attractors were steady states. - Biological filtering: Filtered steady states by phenotype coherence reflecting RA context: for M1, Proliferation=ON, Apoptosis=OFF, Osteoclastogenesis=ON; for M2, Proliferation=OFF, Apoptosis=ON (capturing increased M1/M2 ratio and bone resorption). All M1 steady states passed; 8192 M2 steady states passed. - Calibration using similarity scores: Compared each filtered steady state to the discretized DE vector using the simple matching coefficient S=(N00+N11)/(N00+N11+N10+N01). Selected steady states with the highest score and computed their mean vector to define calibrated model states. • M1: 384 top-scoring steady states; 222/233 nodes fixed; 11 unfixed; matched 104/105 DE nodes (99%); mismatch: CASP7 predicted 0 vs observed 1. • M2: 96 top-scoring steady states; 158/169 nodes fixed; 11 unfixed; matched 85/88 DE nodes (96.5%); mismatches: BCL2L1, MCL1 (observed up, model 0), and CASP3 (observed down, model 1). - In silico perturbations: Queried Therapeutic Target Database (TTD) for model-present targets with inhibitors (1,643 candidates screened). Identified 71 targets in M1 model and 60 in M2 model. Simulated single KOs and pairwise KOs from calibrated states; compared perturbed phenotype states to calibrated states. Also tested all double KOs of receptors present in the models (M1: 406 pairs; M2: 300 pairs).

- Large-scale executable models and calibration: Successfully generated RA M1 (233 nodes, 64 inputs post-focus) and RA M2 (169 nodes, 39 inputs) Boolean models from curated maps; computed all steady-state attractors for 2^21 (M1) and 2^15 (M2) input combinations under HPC parallelization. After phenotype-based filtering and calibration, models reproduced 99% (M1: 104/105) and 96.5% (M2: 85/88) of differentially expressed molecules, with limited mismatches (M1: CASP7; M2: BCL2L1, MCL1, CASP3). - Disease-coherent steady states: All M1 steady states matched RA phenotype criteria; 8,192 M2 steady states passed filtering. - Predicted therapeutic interventions (single targets): • NFkB inhibition: selectively induced M1 apoptosis and suppressed M1 proliferation; no effect on M2 phenotypes. • ERK1 inhibition: suppressed M1 proliferation; reduced proinflammatory cytokines (CCL2, CSF2, IFNG, IL-18, IL-1, IL-6, TNF); blocked VEGFA synthesis in M2 model. • GSK3B inhibition: in M2, induced proliferation and suppressed apoptosis, promoting anti-inflammatory macrophages. - Predicted synergistic combinations (M1 model): Among 2,485 combinations, two pairs showed synergy on M1 phenotypes: • JAK1 + JAK2: together induced M1 apoptosis and suppressed proliferation (single KOs ineffective alone). • ERK1 + Notch1: combined KO induced M1 apoptosis and suppressed proliferation (ERK1 alone suppressed proliferation only). - Receptor double KOs: Among 406 (M1) and 300 (M2) receptor KO pairs, none perturbed apoptosis or proliferation phenotypes, highlighting pathway redundancy/crosstalk. - Computational performance: Utilized BMA on HPC (96 CPUs, 768 GB RAM). All identified attractors were steady states under synchronous updates.

The study demonstrates that automatically inferred, large-scale Boolean models of RA synovial macrophages can be calibrated to disease- and cell-specific behaviours using curated literature and transcriptomic data, enabling phenotype-aligned predictions. The models recapitulate the elevated M1/M2 imbalance and osteoclastogenic context of RA and reproduce most observed differential expressions. Key predictions identify NFkB as a selective suppressor of proinflammatory M1 macrophages, consistent with its known central role in M1 cytokine programs and with therapeutic effects of NFkB pathway modulation. GSK3B inhibition is predicted to promote M2 macrophage survival and proliferation, aligning with literature on anti-inflammatory macrophage polarization and CREB1-mediated survival signaling. Synergistic effects of ERK1+Notch1 and JAK1+JAK2 on inducing M1 apoptosis and suppressing proliferation suggest combination strategies that go beyond single cytokine blockade and could translate to improved macrophage-targeted therapy in RA. The absence of phenotype changes upon receptor double KOs underscores robust intracellular pathway redundancy and crosstalk buffering extracellular receptor perturbations. Methodologically, leveraging HPC-enabled BMA allows exhaustive steady-state exploration despite large input spaces. While focus on steady states facilitates comparison with expression data and stable phenotypes, it may miss oscillatory dynamics. Discrepancies in apoptosis-related components reflect the use of macrophage- and RA-specific but not phenotype-specific datasets; future phenotype-specific omics should refine calibration. Overall, the framework supports hypothesis generation and prioritization of macrophage-selective interventions in RA.

This work presents a scalable framework to convert curated disease maps into executable large-scale Boolean models, compute and filter attractors, calibrate against disease- and cell-type-specific evidence, and simulate therapeutic perturbations. Applied to RA synovial macrophages, the calibrated M1 and M2 models reproduced most differential expression signals and RA-relevant phenotypes, and predicted interventions: NFkB inhibition to selectively suppress proinflammatory M1 macrophages; GSK3B inhibition to promote anti-inflammatory M2 macrophages; and synergistic combinations ERK1+Notch1 and JAK1+JAK2 to induce M1 apoptosis and inhibit proliferation. Future work includes integrating additional and phenotype-specific datasets to improve calibration, representing macrophage activation on a continuum beyond the M1/M2 dichotomy, and constructing a multicellular RA model by coupling macrophage maps with fibroblast and Th1 cell maps from the RA-Atlas to assess system-level therapies and scalability.

- Data specificity: Calibration relied on RA and macrophage-specific datasets that were not phenotype-specific (M1 vs M2), contributing to mismatches in apoptosis-related nodes (e.g., CASP7, BCL2L1, MCL1, CASP3). Access to RA M1/M2 phenotype-specific omics would improve accuracy. - Modeling granularity: Boolean formalism necessitated discretization of expression data and excluded non-differentially expressed components; multi-state/fuzzy approaches could add granularity but increase computational complexity. - Attractor scope: Only steady states were analyzed; excluding cyclic attractors may miss relevant oscillatory/periodic dynamics. - Map and signature coverage: Limited overlap with recently proposed macrophage signatures and potential gaps in RA-specific annotations can affect comprehensiveness. - Redundancy and receptor effects: Lack of phenotype changes upon receptor double KOs may reflect strong pathway redundancy; experimental validation is needed to confirm predictions. - Computational demands: Exhaustive input combination analysis scales exponentially; requires HPC resources and careful optimization to remain tractable for larger multicellular models.