Nucleus pulposus cell network modelling in the intervertebral disc

In 2020, the global prevalence of low back pain (LBP) exceeded half a billion cases. This condition was responsible for 7.7% of all Years Lived with Disability (YLDs), making it the leading cause of disability worldwide¹. IDD is one of the main causes of LBP, as it explains between 26% and 42% of cases². The IVD is the largest avascular organ in the human body, located between the vertebral bodies and consists of three main tissues: the nucleus pulposus (NP); the annulus fibrosus (AF); and the cartilage endplate (CEP). Due to the avascular nature of the organ, the nutrient supply to NPC is made via diffusion through the interstitial fluid of the ECM, resulting in a slow tissue turnover. IDD is multifactorial with several inter-relating risk factors including genetic inheritance, age, metabolic factors, mechanical loads, and other environmental factors, leading to altered cellular phenotypes. These alterations translate into a disruption of the balance between the anabolic and catabolic processes that regulate the ECM of the disc, tilting the balance towards catabolism (see Fig. 1), ultimately contributing to the onset of LBP.

Structure and main ECM components of a non-degenerate, IVD (left) and a degenerate IVD (right). Within a non-degenerate human IVD the anabolic and catabolic components that regulate the ECM of the disc are kept in balance. During disc degeneration the balance is dysregulated, resulting in decreasing matrix synthesis of COL2A and ACAN and promotion of degrading enzymes (MMPs, ADAMTs). Regarding the disc morphology clear boundaries between NP and AF are difficult to distinguish with degeneration. Created with BioRender.com (2022).

Full size image

During disc degeneration the native NP and AF cells of the IVD can acquire pro-inflammatory phenotypes³ that might trigger pro-catabolic events with a strong inertia over time. The early stage of degeneration involves a cascade of cellular changes resulting in loss of proteoglycans and thus associated dehydration, in the NP.

Several experimental studies have been conducted to explore the biochemical environment of the IVD, which is responsible for the onset and the progression of its degeneration. Some have focused on the catabolic effect of pro-inflammatory cytokines on the ECM^4,5,6,7 and how the inhibition of these cytokines could help in IVD treatment^8,9,10. Other studies have explored the role of anti-inflammatory cytokines, such as interleukin-4 (IL-4) and interleukin-10 (IL-10)^11,12,13, and the potential role of growth factors, including growth differentiation factor 5 (GDF5)¹⁴ in the alleviation of IDD through promotion of matrix synthesis, cell proliferation and differentiation, reduced cell apoptosis, and tissue regeneration.

Despite the avail of anti-inflammatory and pro-anabolic factors, it is difficult, however, to understand the initiating factors (disc microenvironment, biomechanics, genetics and epigenetics, bacterial infection of the disc and the gut microbiome) involved in the shift from anabolism to catabolism, by using only in-vivo or in-vitro techniques¹⁵. Computational methods, like Finite Element Models (FEM), offer valuable insights by simulating the mechanical loading of the IVD structure in both healthy and degenerated states. These models can predict changes in disc tissue structure, composition, and properties, including cellular activities such as nutrition^16,17,18,19. However, FEM requires the use of constitutive models and laws with well-defined quantitative parameters, posing challenges when modelling targets at cellular and/or sub-cellular levels. Given the necessity of capturing regulatory mechanisms at these finer scales to understand IDD, modelling approaches must grapple with the intricacies of these interactions, often relying on knowledge or semi-qualitative biological measurements.

According to such a need, biological networks offer graph-based knowledge representations in a holistic and intuitive way. In these models, proteins are depicted as nodes interacting through activation and/or inhibition links. Static interaction maps are transformed into computable models to predict cellular responses, either through logical relationships among nodes and links, or via ordinary differential equations (ODEs). When knowledge needs to be mapped to cope with the lack of specific laws, ODE can be derived through Boolean rules or fuzzy-logic interpolations thereof^20,21. The latter can provide information about relative (or pseudo levels of) protein concentrations or expression, on a continuous normalized scale, and effectively recognises key cell regulators responsible for the ECM disruption. Indeed, numerous studies employing cell network models have aimed to pinpoint key molecules in IDD, with potential applications in diagnosis, disease monitoring, and therapeutic interventions^22,23,24,25.

Xu et al.²² provided fresh insights into the regulatory network of IDD, by integrating proteomics and transcriptomic profiling data from degenerated human NP cells²². Their analysis reveiled six key regulators (Tumor Necrosis Factor (TNF) Alpha Induced Protein 6 (TNFAIP6), Chitinase-3-like protein 1 (CHI3L1), Keratin 19 (KRT19), Dermatopontin (DPT), Collagen Type VI Alpha 2 (COL6A2) and Collagen Type XI Alpha 2 (COL11A2)) with potential functional roles in the IDD process, and they identified two regulators already reported (KRT19, COL6A2) as important IDD markers in independent studies. Li et al. ²³ conducted a comprehensive genome-wide study on RNA expression profiles, constructing protein and disease-gene interaction networks, thereby identifying three novel IDD-specific genes²³. Moreover, they identified entrectinib and larotrectinib as potential treatments for IDD via the NTRK2 gene²³. Huang et al.²⁴ constructed a regulatory mechanism network of lncRNA/circRNA– miRNA–mRNA in IDD²⁴. They proposed potential interaction axes, such as NEAT1–miR-5100–COL10A1 and miR663AHG/HEIH/hsa-circ-0003600–miR-4741–HAS2/HYAL1/LYVE1, shedding light on molecular mechanisms in IDD, notably implicating type X collagen synthesis and hyaluronic acid metabolism imbalance in degeneration. Finally, a top-down network modelling approach was proposed to estimate NPC responses in multifactorial environments^21,26, by systematically translating experimental data into feed-forward model parameters, facilitating high-level simulations of cell activity based on factors such as cell nutrition and pro-inflammatory activity.

Despite notable progress, there remains a lack of models to capture the regulation of NPC activity or phenotype, which is complex due to the multitude of pro-anabolic and pro-catabolic factors, which impact on the complex biochemical environments within the NP in non-degenerate IVDs and during IDD. Within articular cartilage, Segarra-Queralt, M. et al.²⁷ has developed a knowledge-driven network, with a verified balance between pro-anabolic and pro-catabolic processes, which was calibrated and validated against relative protein concentrations²⁷. However, to date no such model has been developed to simulate the complex interactions and regulation within the IVD.

Hence, this study aimed to develop an extensive knowledge driven network model to simulate IVD cell activity, augmenting current approaches in IVD systems biology. Our approach prioritizes high-level modelling, focusing solely on the effects of proposed interactions among soluble proteins, without delving into the intricacies of intracellular signalling mechanisms. This facilitates assessments against falsifying experiments and retains the bare essential descriptors of normal and perturbed cell activity, while preserving model interpretability²⁷. Furthermore, the present network modelling approach provides a curated corpus, shared hereby, that uniquely reflects the activity of NPC in terms of regulated soluble factors, according to current knowledge about non-degenerate IVD and IDD.

Overview

An initial literature-based network model was meticulously crafted, capturing the complex anabolic and catabolic protein interactions within NPC, critical for ECM regulation in the IVD. After first simulation assessments, this model underwent enrichment by integrating general protein–protein(p-p) interactions sourced from the Homo sapiens dataset within the String database (Szklarczyk, D., Gable, A. L., Nastou, K. C., Lyon, D., Kirsch, R., Pyysalo, S., Doncheva, N. T., Legeay, M., Fang, T., Bork, P.‡, Jensen, L. J.‡, & von Mering, C.‡ (2021). The STRING database in 2021: customizable protein–protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Research, 49(D1), D605-D612. doi:10.1093/nar/gkaa1074) along with additional information about closely related cells such as chondrocytes. Following this enrichment process, the model successfully portrayed the anticipated NPC baseline behaviour within a non-degenerate IVD. To further validate its qualitative accuracy, the model underwent thorough testing against independent experimental data concerning the activity of human non-degenerate IVD NPC in in vitro cell culture treatments that favoured either anabolic or catabolic responses. Moreover, a sensitivity analysis was conducted through a full factorial design of experiments to pinpoint the soluble regulators exerting the most significant impact on structural proteins and degrading enzymes within the network. Specifically, emphasis was placed on identifying cytokine and growth factor nodes with the greatest influence (Fig. 2). In essence, this approach yielded a unique and knowledge-driven NPC model capable of replicating expected cellular activity within the context of a non-degenerate IVD. The model was shared via the BioModels.net and the underlying literature corpus was shared through Zenodo repository^28,29,30.

A literature review was firstly conducted to collect all the information about biochemical stimulation and activity of NPC in IVD regulation. Then, a knowledge-based RNM was built, and its baseline was calculated and analysed. The following enrichment method fed continuously the network with new connections from general protein–protein interactions in Homo Sapiens. The RNM was re-assessed, and a sensitivity analysis identified the cytokines and growth factors with the highest impact in the system, which was again compared to independent literature knowledge.

Full size image

Knowledge-based RNM

From the literature review, we gathered information about soluble proteins that were shown to affect NPC activity. The initial topology consisted of 31 different nodes (proteins) and 59 edges that represented node interactions in terms of activations or inhibitions (Fig. 3). Full network detail can be found in the following repository²⁸ and in Supplementary Table 1. Pro-inflammatory cytokines and degrading enzymes had the highest degree of connectivity (DoC), with Interleukin-1 beta (IL-1β) being dominant with 28 DoC. Conversely, only a few connections linked the anti-inflammatory/anti-catabolic and pro-anabolic regulators such as IL-10 (2 DoC), IL-4 (5 DoC), TGF-β (1 DoC), Insulin-like growth factor (IGF) (1 DoC), GDF5 (4 DoC) and Tissue Inhibitor of Metalloproteinase (TIMP) (2 DoC), to the rest of the network.

The network consists of high-level p-p interactions, including the most important soluble proteins in IVD regulation. These are connected through activation (green lines) and inhibition (red lines) links. The size of the node is proportional to the number of connections, i.e. nodes with bigger size have more connectivity.

Full size image

Node activations were calculated using Mendoza’s method, which applies fuzzy logic to bridge the gap between Boolean logic and ordinary differential equations²⁰. This semi-quantitative approach allows for more detailed representations of biological regulation, capturing continuous levels of (normalised) activity in the network (see the Methods for additional details). The exact code used can be found in the corresponding repository²⁹ while the model was deposited in BioModels³⁰ and assigned the identifier MODEL2411040001 with which can be accessed. After conducting 100 simulations to ensure stability, we examined the distribution of the resulting steady states (SS) for each protein. These distributions were skewed, with small variance of the SS values. To enhance clarity for readers, the median of the baseline was represented using bar plots (see Fig. 4), while the detailed boxplot analyses are available in the Supplementary Fig. 1. The activation of pro-catabolic or abnormal matrix components such as Collagen Type X Alpha 1 (COL10A1), Collagen Type I Alpha (COL1A), Interferon gamma (IFN-γ), TNF, Interleukin-1 alpha (IL-1α), IL-1β, Interleukin-6 (IL-6), C-C Motif Chemokine Ligand (CCL), Matrix Metalloproteinase-3 (MMP3), Matrix Metalloproteinase-9 (MMP9), Matrix Metalloproteinase-13 (MMP13), A Disintegrin and Metalloproteinase with Thrombospondin motifs-4/5 (ADAMTS4/5) and Vascular Endothelial Growth Factor (VEGF) was notably abundant. Conversely, only a few anabolic / anti-catabolic components, including IL-4, IL-10, IGF1 and GDF5 were observed to be activated. Interestingly, the activation levels of IVD ECM proteins nodes, Collagen Type II Alpha (COL2A) and Aggrecan (ACAN), were calculated to be very low.

The baseline of the knowledge-based regulatory network is represented through bar plots and reflects a degenerate state of the disc. Activation level (0 for no activation and 1 for the highest activation), of each regulator is given in the SS of the system and are colour-labelled depending on their biochemical properties in the IVD regulation.

Full size image

Enriched RNM

The first objective was to build a regulatory network that reaches a SS representative of a non-degenerate IVD, i.e., with a SS characterized by high pro-anabolic node activity and low pro-catabolic activations. However, the initial topology converged to a clearly pro-catabolic SS (Fig. 3). This motivated the enrichment process through the STRING database, aimed at transitioning the network from representing a pro-catabolic state to a pro-anabolic state of the IVD. The SS in our network are influenced by three main factors: (1) the network topology, which serves as the primary determinant; (2) the initial node activations that were randomly varied over 100 sets of initial conditions, to test the robustness of the topology effect, and to the same attractor in 90% of the cases, i.e. reflecting the deterministic importance of the topology; (3) sustained activation (clamping) of specific nodes, which can shift the network to alternative SS, representing perturbed phenotypes. Importantly, these perturbed phenotypes have been found to correspond to stable attractors within the network, indicating their robustness against perturbations. More detailed information about these methods can be found in the Methods section.

As a result, general p-p interactions in Homo sapiens, along with relevant interactions in chondrocyte regulation, were incorporated (Supplementary Table 2). The topology of the enriched network (found in the corresponding repository²⁸), highlighted the importance of the potential missing links, including anti-inflammatory cytokines such as IL-10 (18 DoC) and IL-4 (19 DoC), growth factors like TGF-β (15 DoC), IGF1 (8 DoC), Cellular Communication Network Factor 2 (CCN2) (3 DoC) and GDF5 (6 DoC), as well as inhibitors of matrix degradation: TIMP1/2 (9 DoC) and TIMP3 (3 DoC) (Fig. 4). The enriched model included 33 different proteins and 153 interactions. Compared to the initial RNM, we included five new proteins (Fig. 5, denoted by bold borders) and 94 new links (Fig. 5, represented by bold lines). Some commonly known or expected interactions such as IL-4 inhibits IL-1β, IL-10 inhibits IL-1β, TGF-β inhibits ADAMTS4/5 and TGF-β inhibits MMP2 could not be found in the IVD-specific literature but were reflected in the STRING database. The updated enriched network was finally able to represent a pro-anabolic activity with high expression level of structural proteins and growth factors and low expression level of (pro)-inflammatory cytokines and degrading enzymes (Fig. 6). The detailed boxplot analyses are available in the Supplementary Fig. 2.

This network encompasses p-p interactions among pivotal soluble proteins involved in IVD regulation, along with general p-p interactions in Homo sapiens. Connections between proteins are depicted through activation (green lines) and inhibition (red lines) links. The size of each node corresponds to the number of its connections, with larger nodes indicating higher connectivity. Newly incorporated nodes are highlighted with bold borders, while bold lines denote newly added links resulting from the enrichment process.

Full size image

Activation level (0 for no activation and 1 for the highest activation) of each regulator is given in the SS of the system and are colour-labelled depending on their biochemical properties in the IVD regulation.

Full size image

Overall corpus results after enrichment

Considering the sources of knowledge out of the literature review, and after the enrichment process, 103 peer-reviewed papers were used for the corpus. From these papers, 38 included information about outcomes of in vitro culture of IVD NPC, yielding a subtotal of 98 p-p interaction links within the network (Table 1). Upon categorizing these experimental sources by the pathological state of the IVD, 56 links were derived from studies examining NPC behaviour in non-degenerate human and animal IVD tissues (Table 2), while 50 links were informed by investigations focusing on NPC behaviour on NPCs extracted from degenerate IVDs.

Assessment of the network

To assess the functionality of the final system and provide an initial independent validation of the model, we tested qualitatively the enriched network against the respective results of two independent experimental treatments of human non-degenerate NPC, reported in the literature^31,32.

Experiment 1: Network perturbation with the pro-inflammatory cytokine TNF. In accordance with findings by Tekari et al.³¹, the anticipated outcomes involved increased expressions of MMP3 and IL-6, along with decreased expressions of ACAN, COL2A, IGF1, and TGF-β. Concurrently, as demonstrated by Millward-Sadler et al.³², the experiments yielded a downregulation of COL2A gene expression, coupled with an upregulation of MMP3, MMP9, and MM13 gene expression. In Fig. 7A, the baseline of the enriched network is depicted by blue bars, while yellow bars represent the baseline following TNF stimulation (for boxplots see Supplementary Fig. 3A). The observed behaviours align qualitatively with the anticipated outcomes from the literature. Furthermore, a Mann–Whitney U test was conducted to assess the significance of the two non-normal distributions. The analysis revealed that all changes were statistically significant (p < 0.05), except for the protein Progranulin (PGRN).

Testing of the RNM by replicating two independent experimental studies made in human healthy IVD NPCs. Initial SS (blue bars) and final SS after corresponding perturbations (yellow bars) with the pro-inflammatory cytokine A TNF and B IL-1β.

Full size image

Experiment 2: Network perturbation with the pro-inflammatory cytokine IL-1β. In line with observations by Millward-Sadler et al.³², the anticipated responses included the upregulation of MMP3, MMP9, and MMP13 gene expression. These expected outcomes are visually represented in our system by the yellow bars in Fig. 7B (for boxplots see Supplementary Fig. 3B), demonstrating the reproduction of protease responses to IL-1β stimulation. Additionally, the Mann–Whitney U statistical test confirmed the significance of all observed changes (p < 0.05).

To further assess our system’s response to key anabolic regulators and explore potential rescue strategies, we simulated two distinct degenerate environments for NPC: one with pro-inflammatory cytokine TNF (Fig. 8A) and another with IL-1β (Fig. 8D), both set at their maximum activation levels. Under the TNF baseline, IL-4 (Fig. 8B) and IL-10 (Fig. 8C) significantly increased ACAN and COL2A levels (p < 0.05) and decreased levels of IL-1β, IL-6, Interleukin-8 (IL-8), Interleukin-12 Alpha (IL-12A), Interleukin-17 Alpha (IL-17A), MMPs, ADAMTS4/5 and VEGF (p < 0.05). Additionally, IL-4 significantly increased IGF1 levels (p < 0.05). When using the IL-1β baseline, TGF-β significantly reduced the expression of MMPs, ADAMTS4/5, and VEGF (p < 0.05) while significantly increased ACAN and COL2A levels (p < 0.05) (Fig. 8E). Notably, GDF5 activation in the IL-1β baseline significantly increased ACAN and COL2A levels (p < 0.05), although it did not decrease catabolic factors (Fig. 8F). When it comes to CCN2, its stimulation under IL-1β treatment, did not have any significant change while its combination with TGF-β did not promote further anabolism that TGF-β alone.

A degenerate baseline was produced to promote catabolism by clamping TNF (A) or IL-1β (D) to 1. Rescue strategies were simulated by stimulating the TNF degenerate baseline with B IL-4, C IL-10, and the IL-1β degenerate baseline with E TGF-β, F GDF5, G CCN2 and H TGF-β and CCN2 (yellow bars).

Full size image

Sensitivity analysis

The final objective was to conduct a sensitivity analysis to identify the mediators with the most significant impact on ECM synthesis proteins (ACAN, COL1A, and COL2A) and degrading enzymes (ADAMTS4/5, MMP1, MMP2, MMP3, MMP9, and MMP13) among the 15 cytokines and 5 growth factors in the network. Simulation results revealed that among all cytokine combinations in the RNM, IL-1Ra had the highest impact on both examined groups (Fig. 9A, B). Additionally, IL-6 and IL-1α showed a smaller direct effect on COL1A and COL2A, while IL-8 affected only COL2A (Fig. 9A). Conversely, MMP1 and MMP2 were influenced by the inflammatory cytokines IL-6 and IL-1α (Fig. 9B). Regarding the growth factors, TGF-β had the most substantial impact on the selected groups of interest (Fig. 10A, B), while GDF5 significantly affected COL1A, IGF1 affected COL1A, MMP3 and MMP13, VEGF affected COL1A and CCN2 affected MMP3 and MMP13 (Fig. 10A, B).

ANOVA analysis results showing the most significant cytokines (green bars) that affected A the structural proteins and B the degrading enzymes are represented by green bars. Non-significant combinations are indicated by pink bars, while combinations that fall outside the range of significance are not depicted. +: positive impact, −: negative impact.

Full size image

Significant combinations of growth factors that impact A structural proteins and B degrading enzymes are represented by green bars. Non-significant combinations are indicated by pink bars, while combinations that fall outside the range of significance are not depicted. +: positive impact, −: negative impact.

Full size image

The disruption of a balanced regulation of the ECM by IVD cells might cause or accelerate the degeneration of the disc. Biochemical factors need to be explored carefully, not only to understand the mechanisms that contribute to ECM degradation, but also to predict and/or prevent the triggering of these mechanisms. Multiple experimental measurements have generated valuable knowledge about IVD protein response against pro-catabolic, anti-catabolic and growth factor stimulations, but the integration of such knowledge is still missing, as targeting only a single protein or few proteins at once is not sufficient to test the effect of the many proteins that can shape the biochemical microenvironment of NPC. Network modelling can provide a holistic and intuitive vision of the complex behaviour and interactions among proteins, and complement in-vitro studies, as it allows to explore perturbations /interactions that would require further experimental exploration.

Different approaches can be adopted to build a network model: either gather literature information about the biochemical IVD regulation (knowledge-driven network) or exploit directly experimental measurements (data-driven network), as performed by Baumgartner, L. et al.^21,26. While the latter can predict IVD cell responses in a pro-inflammatory environment, in terms of protein expression, it is limited on TNF and IL-1β stimulations and it did not consider important pro-anabolic and pro-catabolic factors in IVD research. On the other hand, our model provides a unique corpus of knowledge able to represent a relatively large spectrum of NP cell activity predictions upon different perturbations. In addition, it can test virtually perturbations still not experimentally explored, which facilitates, and identifies priorities for future experimental study designs.

The initial system exhibited an anticipated catabolic response, given that 52% of the experimental samples sourced from the literature either originated from degenerated discs or were artificially induced into a degenerative state prior to testing. However, our primary aim was to construct a network that accurately represented a non-degenerate IVD, enabling us to evaluate how various stimuli influence metabolic shifts. Consequently, enriching the network with relevant p-p interactions, particularly those aligned with our objectives, was deemed a logical progression. By incorporating key missing interactions (depicted as bold lines in Fig. 5), we successfully attained an anabolic state. This state is distinguished by heightened activation levels of essential ECM components such as ACAN and COL2A, along with pivotal growth factors including TGF-β, CCN2, IGF1 and GDF5. Conversely, activation levels of angiogenic factors (e.g., VEGF), catabolic markers (MMPs, ADAMTS4/5), and pro-inflammatory cytokines (such as interleukins, TNF, and IFN-γ) were notably suppressed. This iterative process underscored not only the efficacy of our approach but also the imperative for further experimental studies specifically focused on IVD biology. Out of the 103 papers analysed, 38 pertained to IVD NP cells, 22 focused on cartilage chondrocyte cells, while the remaining 43 addressed other cell types, highlighting the pressing need for more comprehensive investigations within the field.

The significance of enrichment and subsequent network predictions finds robust support in the experiments conducted by Li et al.³³. The authors investigated the impact of treating degenerate NP canine cells with IL-10 and TGF-β, revealing a notable suppression in the release of inflammatory cytokines IL-1β and TNF³³. Such effect shall further limit the additional catabolic expression of proteases and deterioration of the ECM^34,35. Ge et al.¹³ showed that IL-10 treatment on NP degenerate cells could significantly increase the protein expression of type II collagen and aggrecan while decreasing protein levels of type X collagen and endogenous mRNA expression of TNF and IL-1β¹². A study that shows the importance of our enriched regulators (IL-4, TIMP1/2, and TIMP3), was demonstrated by Kedong et al.¹¹, who determined the anabolic role of the anti-inflammatory IL-4 against the inflammatory markers IL-6 and IL-8¹¹ while TIMPs act as a balance factor to the catabolic behaviour of MMPs and ADAMTs^36,37. Notably, TIMP3 emerges as a potent inhibitor of ADAMTS4 and ADAMTS5, playing a pivotal role in safeguarding against ECM matrix-degrading enzymes³⁸. Moreover, despite the concurrent increase of TIMP1, a major MMP inhibitor, alongside MMPs during disc degeneration, TIMP3 exhibits a distinct behaviour, not following the trajectory of cell immunopositivity for ADAMTS during degeneration³⁹. Another inclusion was CCN2, which have been reported as possible regenerative molecule in the IDD when combined with its downstream mediator TGF-β⁴⁰. Lastly, the positive outcomes of enrichment are further exemplified by the addition of IL-1Ra, the natural inhibitor of IL-1⁴¹. When delivered through gene therapy, this receptor antagonist could significantly diminish the expression of MMP1, MMP3, MMP13 and ADAMTS4, as well as the activity of these enzymes in the NP and the AF, by more than 95%⁴¹.

In our pursuit of independently assessing the enriched network, we conducted a thorough literature search to identify experimental studies not initially incorporated into the knowledge representation that formed the network topology. The literature was filtered to solely target stimulated NP cells derived from non-degenerate human IVD samples and build a set of tests as relevant to human IVD biology, and coherent, as possible. First, we stimulated the network with the pro-inflammatory cytokine TNF. This stimulation prompted a significant catabolic response within our system, evidenced by the upregulation of IL-6 (p < 0.05) and MMP3 (p < 0.05), accompanied by the downregulation of ACAN (p < 0.05), COL2A (p < 0.05), IGF1 (p < 0.05), and TGF-β (p < 0.05). Notably, this trend aligned with findings reported by Tekari et al.³¹ as well³¹. When comparing the nodal activation level under TNF stimulation with the second experimental study from Millward-Sadler et al.³², we observed consistent upregulation of MMP3 (p < 0.05), MMP13 (p < 0.05), IL-1β (p < 0.05), and IL-1Ra (p < 0.05), alongside downregulation of COL2A (p < 0.05) and COL1A (p < 0.05). Secondly, we perturbated the network with the pro-inflammatory cytokine IL-1β, to confirm the upregulation of MMP3 (p < 0.05), MMP13 (p < 0.05), MMP9 (p < 0.05), and TNF (p < 0.05) as measured from Millward-Sadler et al.³². Essentially, our in-silico simulations corroborated these experimental results, providing further validation of our enriched networks’ predictive capabilities.

In another aspect of our assessment, we sought to explore the potential effects of IL-4, IL-10, TGF-β, GDF5 and CCN2, all recognized for their reported anabolic responses in a degenerative environment^{11,12,33,42,43,44,45,46,47,48}. To establish a degenerate baseline, we used two approaches: clamping TNF at a level of 1 and clamping IL-1β at a level of 1. These cytokines are proposed to play a key role in intervertebral disc degeneration (IDD), as indicated in previous studies^4,32,49,50. Baseline characterized by elevated expression levels of MMPs, ADAMTS4/5, and VEGF, all associated with ECM degradation and angiogenesis^32,39,51,52, were thus created. We then tested the effects of anti-inflammatory cytokines and growth factors on these baselines. We saw that anti-inflammatory cytokines could rescue the IDD when applied in an environment like the TNF baseline, while growth factors (TGF-β better than GDF5) would promote anabolism when applied in environment similar to IL-1β baseline.

More specifically, Ge et al.¹³ has previously demonstrated that IL-10 could be used to delay IDD through its anti-inflammatory response by inhibiting p38 MAPK pathway activation¹², while recently Kedong et al.¹¹ has proposed that the recombinant IL-4 delivery into the IVD might have a beneficial therapeutic effect by reducing disc inflammation following herniation¹¹. In our in silico simulations, IL-4 decreased MMP3 and IL-17A, consistent with in vitro results from Chambers et al.⁵³ and Sandoghchian et al.⁵⁴. Furthermore, IL-4 increased ACAN, COL2A, and IGF1, aligning with experimental results in chondrocyte cells^55,56. IL-10 also increased COL2A, a result observed experimentally in NPC¹² and decreased pro-inflammatory cytokines (IL-6, IL-8, IL-12A, IL-17A), corroborating findings in other cell types^57,58. However limited studies have investigated the effects of IL-4 and IL-10 on human NP cells from non-herniated material which could be infiltrated with immune cells. TGF- β stimulation led to depletion of MMP2 expression, a result reported previously in bovine IVD NP cells by Pattison and colleagues⁵⁹. It also eliminated the catabolic markers MMP3, MMP9, MMP13, and ADAMTS4/5, whose activity is closely related to ECM degradation⁵. Although the downregulation of these markers by TGF-β has been observed in general human cell types^27,60, such findings are lacking in NP cells, underscoring the need for further experimental studies specific to IVDs. Notably, TGF-β induced a high expression level of TIMP3, a result not previously reported in NP cells but observed in fibroblasts and chondrocytes^61,62, which share similar functional and cellular characteristics with NP cells^63,64,65. Additionally, TGF-β significantly increased structural proteins ACAN and COL2A, consistent with findings in NPCs by Masuda et al.⁶⁶. GDF5, while not significantly altering the degenerate state, did decrease TNF expression, consistent with findings by Guo et al.^43,67, suggesting GDF5’s efficacy in suppressing ECM degradation. Our in silico model also promoted anabolism by increasing ACAN and COL2A, as found by Chujo et al.¹⁴ and Le Maitre et al.⁴⁴. Whilst CCN2, could not reverse the catabolic effects of IL-1β in this model, it has been reported to have regenerative potential when combined with TGF-β40. However, further testing in human NP cells is necessary, as the current findings are based on experiments in rat discs. CCN2 is a matricellular protein present in the extracellular matrix (ECM) that plays a critical role in signalling mechanisms⁶⁸. Therefore, the involvement of signalling pathways must be considered, particularly since TGF-β is known to positively regulate CCN2 through the Smad signalling pathway⁴⁶.

Here, we also investigated the potential rescue effects of IL-4, IL-10, TGF-β and GDF5 in a degenerate baseline generated from IL-1β and TNF co-stimulation. Only TGF-β was able to decrease the catabolic effects on the MMPs, while no significant changes were observed in structural proteins to promote anabolism (Supplementary Fig. 4). The expression of TNF and IL-1β in degenerate IVDs is complex and can vary independently. It is unlikely that IVDs would express only one cytokine (IL-1β/TNF), however depending on the stage of degeneration and infiltration of inflammatory cells following AF and CEP fissures, it is possible that some discs could have differential levels^32,69. Understanding these dynamics is crucial for developing targeted therapies for disc degeneration.

The sensitivity analysis aimed to determine key determinants of the ECM regulation in the disc. Among the cytokines the anti-catabolic factor IL-1Ra had the most impact, showing the importance of the natural IL-1 inhibitor, whose deficiency might develop IDD⁷⁰. During disc degeneration COL2A and ACAN are decreased⁷¹, which is associated with dehydration of the disc and reduced ability to withstand loads. Thus, if IL-1Ra via the inhibition of native IL-1 agonists (α and β) could prevent the decrease in matrix synthesis, which accompanied with growth factor stimulation, either native or applied could promote new matrix synthesis, IDD could be halted, and potential regeneration induced. However, the realization of such therapeutic potential hinges upon a deeper understanding of IL-1Ra’s regenerative properties, necessitating further experimental investigations. Moreover, sensitivity analysis underscored the profound impact of transforming growth factor-beta (TGF-β) on both structural proteins and matrix-degrading enzymes. This finding emphasizes the pivotal role of TGF-β in orchestrating ECM homeostasis and suggests its potential as a therapeutic target for IDD management. Moreover, in-silico co-stimulation of TGF-β with IL-1β and TNF was able to reduce the presence of pro-inflammatory cytokines and increase ACAN, explaining its strong impact. These insights not only deepen our understanding of the molecular mechanisms underlying disc degeneration but also highlight promising avenues for therapeutic intervention aimed at promoting disc regeneration.

To conclude, the proof-of-concept development of an in silico network model demonstrates the importance of a baseline in network modelling. We have achieved a unique corpus of IVD biochemical regulators and a baseline which is proposed to represent the non-degenerate state of the disc. Our network showed that TGF-β is promising in the regeneration of the disc, while more exploration of the IL-4 and IL-10 anti-inflammatory effects is needed. Considering rescue strategy simulations, the combination of autologous anti-inflammatory serum containing TGF-β, GDF5, IL-10, and IL-4 emerges as compelling potential biological treatment for IDD. In addition, the initial state of the nodes plays a very important role in the response of the system against stimulations and for that more IVD specific experimental studies should be conducted. Current treatments for IDD are mostly conservative and can only alleviate the pain, but not regenerate the disc⁷². It is essential to acknowledge that while our model provides valuable insights, it does not capture intracellular changes or include pathways, which are crucial for comprehensive analysis. Future iterations of our regulatory network model should incorporate these elements to serve as a versatile tool for guiding experimental studies, predicting outcomes, and testing IDD drugs.

Corpus

As a first step, we searched in the open database PubMed for IDD chondrocyte p-p interactions. Specifically, our search terms included the following combinations:

• (Intervertebral disc degeneration) AND (nucleus pulposus cell)

• (Intervertebral disc degeneration) AND (cytokine)

• (Intervertebral disc degeneration) AND (anti-inflammatory cytokine)

• (Intervertebral disc degeneration) AND (growth factor)

• (Intervertebral disc degeneration) AND ((Matrix metalloproteinases) OR(MMP))

• (Intervertebral disc degeneration) AND ((a disintegrin and metalloproteinase with thrombospondin motifs) OR(ADAMTS))

• (Intervertebral disc degeneration) AND (chemokine)

From the plethora of articles identified, a meticulous review yielded 34 peer-reviewed studies that formed the foundation for constructing our initial IVD RNM. This compilation encompassed data pertaining to direct activation, inhibition, upregulation, and downregulation effects among a diverse array of molecular entities crucial in regulating IVD function. These entities spanned structural proteins, degrading enzymes, angiogenic factors, pro-inflammatory, inflammatory, and anti-inflammatory cytokines, as well as chemokines.

To clarify how the literature review informed the development of our network, we provide the following example. A well-supported interaction between IL-1β and COL2A was consistently identified, with multiple studies reporting that IL-1β downregulates COL2A expression, thus facilitating its straightforward inclusion in the model. Though, we identified two conflicting interactions. The first was between IL-1β and IL-10, with studies showing both upregulation and downregulation. We prioritized data from human nucleus pulposus cells cultured in 3D, reflecting the physiological context of IVD biology. The second conflict involved IL-1β and ACAN; while low doses of IL-1β were linked to increased ACAN expression, we emphasized significant effects at higher doses to clarify the transition from a normal disc to a degenerative state. This systematic approach ensured that our network was constructed using the most reliable and contextually relevant data, enhancing its robustness and relevance to IVD pathology. In refining our dataset, we opted to streamline the representation of CCLs by merging them into a single category denoted as CCL. However, we noted a distinct behaviour exhibited by CCL22 compared to others within this category. While IL-1β activated other CCLs, it paradoxically inhibited CCL22⁷³, warranting its separation from the merged group. This final dataset is accessible in the corresponding repository²⁸. Leveraging the Cytoscape_v3.8.2 platform, we employed a node-edge visualization approach to depict proteins as nodes and illustrate their activation and/or inhibition effects through connecting edges. This strategy offers an intuitive and comprehensive representation of the intricate regulatory mechanisms governing IVD biology, facilitating in-depth analysis and interpretation of IDD-related protein interactions.

Semi-quantitative system resolution

To transition the knowledge-based regulatory network (RN) from a static snapshot to a dynamic representation, we employed a method proposed by Mendoza, which utilizes fuzzy logic semi-quantitative interpolation of boolean solutions²⁰. This approach enables the conversion of boolean logic into ordinary differential equations, facilitating the modelling of dynamic behaviour within the network.

The calculation of the final expression for each node involves a set of ordinary differential equations, with the dynamics described by the following equation:

Here, the left-hand side represents the activation function, while the right-hand side depicts a decay function proportional to the node x_n, weighted by a factor γ > 0. The activation term incorporates the sigmoid function ω, which computes the total input of a node x_n.

Generally, each node has different connectivities, according to the number of activators and/or inhibitors (Fig. 11). Thus, ω is given by the following equations, each of them describing the possibility of each node having solely activators, solely inhibitors or a combination of both. The equation is as follows:

(2)

In a regulatory network model, proteins are represented by nodes that interact among each other through activation (red lines) and/or inhibition (green lines) edges.

Full size image

0 ≤ x_n ≤ 1

0 ≤ ω_n ≤ 1

h_n, α_nk, β_nl, γ_n > 0

{({x}_{{nk}}^{{alpha }})} is the set of activators of x_n

{({x}_{{nl}}^{i})} is the set of inhibitors of x_n

is used if x_n has only activators

is used if x_n has only inhibitors

is used if x_n has both activators and inhibitors

In these equations α > 0 and β > 0 are the weights of activators and inhibitors, respectively and take the values of α = 1, β = 1 for simplicity as in spite of their values, they do not affect the monotonically behaviour of the ω function between its intervals. Finally, we choose h = 10 and γ = 1 for all nodes, as proposed by Mendoza et al. (2006), a choice serving the convenient and smooth biological behaviour of the system.

The system was solved using the fourth-order Runge-Kutta method in Matlab 2018b (The MathWorks, Inc, Massachusetts, USA) with the ode45 solver to determine the global attractor or SS of the system (Fig. 12). This SS represents a practical approximation of an artificial time period within a cellular process, likely confined to a localized region of interest. Consequently, we designated this SS as the baseline of each protein within the network. To ensure the robustness of our baseline calculations, we conducted 100 simulations for each protein, varying initial conditions based on Mersenne Twister 0 seeds. This rigorous approach allowed us to confirm the consistency of the baseline across various initial conditions. Additionally, it ensured that each node converged to a SS, reflecting the system’s stable behaviour. While the SS values did not conform to a normal distribution, indicating variability in protein levels, their distribution remained unimodal, being simply right- of left-skewed, depending on the node. Hence, this does not undermine the existence of a global attractor, which signifies a stable state of the system. In our simulations, approximately 90% of the runs converged to this same global attractor, reinforcing the consistency of the system’s behaviour, and the deterministic capacity of the defined topologies despite the observed variability in protein expression. In response to this observation, we opted to employ a boxplot visualization technique to illustrate all SS values comprehensively. This approach offers a detailed representation of the baseline, capturing the distribution of SS values along with key statistical metrics such as median, quartiles, and outliers. However, owing to the limited variance in the values and the skewed distribution, boxplots were not sufficiently sensitive in highlighting significant differences. To enhance clarity and facilitate better interpretation for readers, we chose to present the median of the baseline using barplots. Meanwhile, the detailed boxplot analyses, encompassing additional statistical insights, can be accessed in the Supplementary Figs. 1 and 2. This approach ensures that readers can easily discern the baseline trends while still having access to comprehensive complete datasets provided by the boxplots. The code can be found in the following repository²⁹.

The dynamics of the RNM were captured through a system of ordinary differential equations (ODEs), where the activation of each node is determined by the number of inhibitors and activators. The ODEs were solved using the fourth-order Runge-Kutta method (RK4), initialized with random system conditions. After 30 s, the system converged to a SS, serving as the baseline for each protein within the network. Notably, 100 simulations were conducted, and all resulting SS values were utilized for subsequent analysis, as illustrated by the boxplots.

Full size image

Enrichment

The enrichment of the regulatory network was done firstly through the STRING (Version 11.5) public database. This procedure retrieved information about protein–protein interactions which we could not find in literature when looking specifically for IVD NPC regulation. Text mining, experiments and databases were selected as interaction sources while the interaction score with highest confidence of 0.900 was applied in the direction of having strong and reliable biological information. We utilized the existing nodes of our network to identify new p-p interaction links. For each potential new link, we reviewed the relevant sources (literature or databases) provided by STRING to determine the directionality of the interactions, identifying which proteins acted as activators or inhibitors of others. If clear directionality was not established for a link, we chose not to add that information to the model, ensuring that only well-supported interactions were represented. Additionally, a manual enrichment process was conducted by scouring PubMed for information on growth factors such as TGF-β, CCN2 and IGF, as there were still gaps in the network regarding these factors. Furthermore, to ensure comprehensive coverage of key proteins implicated in IVD pathogenesis, additional proteins were manually added including MMP1, MMP2, TIMP1, CCN2 and IL-1Ra based on recent investigations documented in the PubMed literature. These proteins play pivotal roles in IVD pathophysiology, and their inclusion was deemed essential for a thorough understanding of regulatory network dynamics.

For the manual enrichments the keywords in PubMed were:

• ((intervertebral disc degeneration)) AND ((transforming growth factor beta) OR (TGF-β)) (20/10/2021)

• ((intervertebral disc degeneration)) AND ((insulin-like growth factor) OR (IGF)) (20/10/2021)

• ((intervertebral disc degeneration)) AND ((Matrix Metallopeptidase 1) OR (MMP1)) (20/10/2021)

• ((intervertebral disc degeneration)) AND ((Matrix Metallopeptidase 2) OR (MMP2)) (20/10/2021)

• ((intervertebral disc degeneration)) AND ((tissue inhibitor of metalloproteinases) OR (TIMP1)) (20/10/2021)

• ((intervertebral disc degeneration)) AND ((Interleukin-1 receptor antagonist) OR (IL-1Ra)) (20/10/2021)

• ((intervertebral disc degeneration)) AND ((CCN2) OR (Connective tissue growth factor)) OR (CTGF) (25/09/2024)

As a last step of the enrichment, we used missing links found in osteoarthritis chondrocytes regarding TIMP1, TIMP2, TGF-β, IL-10 and IL-4²⁷, as osteoarthritis is a disease with common mechanisms as IDD this was a similar pathway and cell type to enhance the network. Detailed information for the enriched network used can be found in the corresponding repository²⁸. In the final enriched network ADAMTS4 and ADAMTS5 were merged as ADAMTS4/5, TIMP1 and TIMP2 as TIMP1/2, while IL-23A and IL-32 were removed to simplify our network as they showed limited connections. Eventually, the stable SS of the enriched network was solved by using the method explained in the semi-quantitative system resolution (Fig. 12) as done for the knowledge-based network initially.

Assessment of the network

After building the network we evaluated its performance by exploring the anabolic or catabolic effect of each protein in IVD regulation (database). A widely accepted principle dictates that a high activation level of structural proteins such as COL2A and ACAN, as well as growth factors, coupled with a low activation level of proteases, ADAMTS, and pro-inflammatory cytokines, promotes anabolism within the IVD. Conversely, a low activation level of structural proteins and growth factors, along with a high activation level of MMPs, ADAMTS, and pro-inflammatory cytokines, indicates a catabolic state.

The experimental system was further tested by comparison to independent experiments from the literature, by replicating theoretically two in-vitro experiments. More specifically, in the first test according to Tekari et al.³¹, in which they used non-degenerate human IVD NP cells collected from trauma patients undergoing spinal fusion surgery without any history of disc degeneration before the operation, the current in silico network was stimulated with the pro-inflammatory cytokine TNF. In the second test, according to Millward-Sadler et al.³², in which they used non-degenerate and degenerate human NP cells obtained from post-mortem samples or from surgery, and graded histologically according to the method of Sive and colleagues⁷¹, we stimulated the network with the pro-inflammatory cytokines TNF and IL-1β. In both tests, we analysed the behaviour of our network post-stimulation and compared it with the results obtained from the corresponding experiments. This comparison allowed us to evaluate the functionality and accuracy of our network model in replicating experimental biological responses.

Sensitivity analysis

A sensitivity analysis was conducted to identify the mediators that most significantly affected the system (Fig. 13). Specifically, the analysis aimed to determine which cytokines/chemokines (CSF2, IFN-γ, IL-1α, IL-1β, Il-1Ra, IL-4, IL-6, IL-8, IL-10, IL-12α, IL-17A, IL-18, TNF, CCL, CCL22,), growth factors (GDF5, IGF1, TGF-β, CCN2, VEGF), and their second-order interactions significantly impact the groups of structural proteins (ACAN, COL1A, COL2A) and degrading enzymes (ADAMTS4/5, MMP1, MMP2, MMP3, MMP9, MMP13).

As a first step of the sensitivity analysis a full factorial design of experiment was made to generate all the possible combinations among the 15 cytokines and the 5 growth factors, with two treatment levels (0 and 1). These combinations were then used as initial conditions to solve the system semi-quantitatively and finally get the corresponding SS for each combination. Finally, by using an ANOVA test, the combinations that altered the initial baseline the most was identified.

Full size image

As a first step, a full factorial design of experiment was made and described by n^p total interactions, when n is the number of levels (0 and 1 in our case) and p is the number of factors/proteins tested. In that way all the possible combinations among the analysed proteins are generated and represent the initial states. More specifically, 2¹⁵ = 32.768 simulated combinations among the 15 cytokines and 2⁵ = 32 simulated combinations among the growth factors were tested. For each combination, the direct effect and the second-order effect of the tested factors to the interest groups was investigated. The second step was to solve again our system semi-quantitatively so many times as the total combinations of each tested group, by using each combination as initial state, to obtain the relevant SS. This led to 32.768 different SS for the cytokines and 32 different SS for the growth factors.

The last step included the analysis of the experiments by using the analysis of variance F statistic (ANOVA F Statistic). Such a method compares the variability between the groups to the variability within the groups. The formula used for the ANOVA F Statistic for each node of the network-based model is:

where MSB is the mean sum of squares between the groups

• cytokines-structural proteins

• cytokines-degrading enzymes

• growth factors-structural proteins

• growth factors-degrading enzymes

and MSW is the mean sum of squares within the groups (within the 15 cytokines and within the 5 growth factors). Thereof, the larger the F value, the greater is the evidence that there is a difference between the group means. To compare if the difference among the group means is statistically different, we computed the p-value. The significant F values (p ≤ 0.05) for the interested groups were plotted in Pareto charts in a logarithmic scale as the distribution for the proteins was very large. In the results only the direct effects of the cytokines can be seen for the sake of simplicity, but the analysis included the second-order interactions. Regarding the growth factors all the effects (direct and second order) can be visualized as they were few.