A large-scale Boolean model of the rheumatoid arthritis fibroblast-like synoviocytes predicts drug synergies in the arthritic joint

Abstract

Rheumatoid arthritis (RA) is a complex autoimmune disease with an unknown aetiology. However, rheumatoid arthritis fibroblast-like synoviocytes (RA-FLS) play a significant role in initiating and perpetuating destructive joint inflammation by expressing immuno-modulating cytokines, adhesion molecules, and matrix remodelling enzymes. In addition, RA-FLS are primary drivers of inflammation, displaying high proliferative rates and an apoptosis-resistant phenotype. Thus, RA-FLS-directed therapies could become a complementary approach to immune-directed therapies by predicting the optimal conditions that would favour RA-FLS apoptosis, limit inflammation, slow the proliferation rate and minimise bone erosion and cartilage destruction. In this paper, we present a large-scale Boolean model for RA-FLS that consists of five submodels focusing on apoptosis, cell proliferation, matrix degradation, bone erosion and inflammation. The five-phenotype-specific submodels can be simulated independently or as a global model. In silico simulations and perturbations reproduced the expected biological behaviour of the system under defined initial conditions and input values. The model was then used to mimic the effect of mono or combined therapeutic treatments and predict novel targets and drug candidates through drug repurposing analysis.

Fig. 1. Construction and simulation of a large-scale, modular Boolean model of RA-FLS for evaluating novel drug combinations.

The RA map was converted into an executable Boolean model using the map-to-model framework described in ref. . Using single-cell omic datasets and literature studies, the RA generic model was subsequently enriched in RA-FLS-specific data. The RA-FLS model focuses on five phenotypes (apoptosis, cell proliferation, inflammation, matrix degradation, and bone erosion) characteristic of RA’s fibroblasts. Individual phenotype-specific submodels and a five-phenotype global model were created. Biological scenarios extracted from the literature were used to evaluate and validate the models’ behaviour leading to some modifications of the original models. The modified RA-FLS model was then used to test mono and combined RA therapies. Drug repurposing analysis and further drug combination simulations led to a panel of suggestions of drug combinations that are predicted to have a favourable outcome (apoptosis active, cell proliferation, inflammation, bone erosion and matrix degradation inactive).

Fig. 2. Number of nodes per module/model and shared components among modules.

a Number of nodes of the phenotype-specific modules and the global model. b Venn diagram of all the five-phenotype-specific components. The core of 191 nodes is shared among all five modules, and only a few are characteristic of the corresponding phenotype-specific module.

Fig. 3. Trap spaces of the modified global model.

Heatmap displaying the trap spaces for the tested biological scenarios in the modified global model (a) and the comparison between expected and obtained values shown with colour codes (b). The y axis shows all the tested scenarios’ names, as mentioned in Table 2, regarding all the five phenotypes as outcomes on the x axis with DEFAULT inputs set to one. ^†Represent scenarios where additional conditions were given as a known biological behaviour* or as model conditions**. Trap spaces colour codes: -1 (unfixed) , 0 (OFF) , 1 (ON) . Expectation graph colour codes: Score: expected value, obtained value; 1: Yes [OFF, OFF & ON, ON] ; 0: No [ON, OFF & OFF, ON (conflict)] , -1: Undefined .

Fig. 4. Calculating continuous time phenotypic probabilities of the selected initial conditions.

a Simulation with IL6_Extracellular_space active and Inflammation as output. The Inflammation phenotype gets activated in the presence of IL-6. b Simulation with PDGFA active and Cell proliferation as output. The Cell proliferation phenotype gets activated in the presence of PDGFA. c Simulation with TNFSF11 (RANKL) inactive, SFRP5 active and Bone erosion as output. The Bone erosion gets deactivated in the presence of SFRP5 and the absence of TNFSF11. d Simulation with MMP1 inactive and Matrix degradation as output. Matrix degradation phenotype gets deactivated in the absence of MMP1. e Simulation with FASLG active, AKT2 active, BID inactive and Apoptosis as output (RA original global model). Apoptosis gets activated in the presence of FASLG, AKT2 and the absence of BID. f Simulation with FASLG active, AKT2 active, BID inactive and Apoptosis as output (RA modified global model). Apoptosis gets inactive in the presence of FASLG, AKT2 and the absence of BID due to the dominant-negative regulator CAV1.

Fig. 5. Oscillatory behaviour of TP53 and the mitochondrial proteins.

a Trap spaces show that along with TP53 and MDM2, mitochondrial proteins like BAX and cytoplasmic protein CASP9 remain unfixed. b Simulation with CC shows TP53 and MDM2 oscillations while all inputs remain active. c Simulation with CC shows oscillatory behaviour of the mitochondrial proteins when all inputs remain active.

Fig. 6

Heatmap showing identified trap spaces regarding the five phenotypes while testing different drug targets with DEFAULT inputs conditions as 1 with the modified version of the RA-FLS model.