Systems biology reveals how altered TGFβ signalling with age reduces protection against pro-inflammatory stimuli

Osteoarthritis (OA) is a degenerative condition caused by dysregulation of multiple molecular signalling pathways. Such dysregulation results in damage to cartilage, a smooth and protective tissue that enables low friction articulation of synovial joints. Matrix metalloproteinases (MMPs), especially MMP-13, are key enzymes in the cleavage of type II collagen which is a vital component for cartilage integrity. Transforming growth factor beta (TGFβ) can protect against pro-inflammatory cytokine-mediated MMP expression. With age there is a change in the ratio of two TGFβ type I receptors (Alk1/Alk5), a shift that results in TGFβ losing its protective role in cartilage homeostasis. Instead, TGFβ promotes cartilage degradation which correlates with the spontaneous development of OA in murine models. However, the mechanism by which TGFβ protects against pro-inflammatory responses and how this changes with age has not been extensively studied. As TGFβ signalling is complex, we used systems biology to combine experimental and computational outputs to examine how the system changes with age. Experiments showed that the repressive effect of TGFβ on chondrocytes treated with a pro-inflammatory stimulus required Alk5. Computational modelling revealed two independent mechanisms were needed to explain the crosstalk between TGFβ and pro-inflammatory signalling pathways. A novel meta-analysis of microarray data from OA patient tissue was used to create a Cytoscape network representative of human OA and revealed the importance of inflammation. Combining the modelled genes with the microarray network provided a global overview into the crosstalk between the different signalling pathways involved in OA development. Our results provide further insights into the mechanisms that cause TGFβ signalling to change from a protective to a detrimental pathway in cartilage with ageing. Moreover, such a systems biology approach may enable restoration of the protective role of TGFβ as a potential therapy to prevent age-related loss of cartilage and the development of OA.


Supplementary details of model construction
The complete model schematic can be seen in S5 Fig, and is a modified expansion of two previously published models [1,2]. These were simplified in ordered to reduce computational burden whilst still fitting the same profiles shown in the original models . In order to achieve this, some assumptions were made which marginally changed the structures. The original IL-1+OSM model was presented in three sections, OSM, IL-1 and pro-MMP activation [1]. As we were only looking at MMP13 mRNA, and did not explore the transition to active MMP, we removed the pro-MMP activation component completely. The remaining two sections were simplified ensuring that they still fitted the same profiles. To achieve this some approximations were made. Three proteins were originally present in the model: Dusp16, Mkp1 and Pp4. As these were upregulated by the same source and only relevant for blocking, we replaced all three of these species with one "Block" species. This blocked all three reactions that the original proteins did. However, each reaction was blocked at a different rate in order to have the same effect that the original proteins would have had. In the original model Traf6 was internalised by Irak2 binding to receptor bound IL-1. We changed the receptor profile for IL-1 to include an internalisation step that replaced the need for Irak2 and Traf6. The phosphorylation and subsequent nuclear translocation of Stat3 was also altered so that it could be represented by a single self-limiting reaction. Finally, the original models included mRNAs for many of the proteins. Apart from cFos and MMP13, these were removed so that protein synthesis could be more simply modelled by one reaction.
The Alk1/Alk5 sections of the model are based on the TGFβ section described in [2]. The most significant change was the removal of the Sox9, aggrecan and collagen components from the original model as these were not of interest in the current study. Mmp2 was completely removed as sensitivity analysis showed it had no effect on the model in this form. The role of active Runx2 was also altered, no longer inducing the production of pro-MMP-13 but instead producing MMP13 mRNA. This change was reasonable as Runx2-induced upregulation of MMP13 has been shown to be at the mRNA level [3]. We also included an additional reaction so that Smad1/5/8 signalling leads to an increase in p38 phosphorylation [4].
The original models had not been created to be combined and as a result they did not share any interactions. Therefore, the combined model was extended to include reactions that link the TGFβ and IL-1+OSM signalling pathways, based on knowledge from the literature. We considered two different mechanisms: inhibition of AP-1 (cJun homodimers or cJun/cFos heterodimers) by JunB; and increased degradation of MMP13 mRNA via TGFβ signalling (see below).

Inhibition of AP-1 by JunB
There is a substantial body of evidence showing that JunB is a direct downstream target of Smad3 [5,6] and serves to antagonise rapid gene transcription [7][8][9]. Studies have shown it can displace cFos in the AP-1 complex and bind to cJun [5,7]. It has also been demonstrated to limit the effects of IL-1, in particular its effect on MMP synthesis [10]. Finally, it has been shown to repress gene expression in epithelialmesenchymal transition in cancer due to upregulation by TGFβ. This demonstrates TGFβ can induce significantly high JunB expression to limit gene transcription [6]. We modelled these effects by adding reactions for JunB production via Smad2/3 and JunB binding to cJun to inhibit AP-1 transcription (see yellow highlighted cells in Table 2).

Increased degradation of MMP13 mRNA via Smad2/3
As the protein or microRNA responsible for increased MMP13 mRNA degradation via Smad2/3 signalling is currently unknown we added a species with ID 'Dummy'. We assumed that transcription of the 'Dummy' gene requires an activator protein 'A' and that 'A' needs to be phosphorylated by Smad2/3 for its activity. The 'Dummy' protein then interacts with MMP13 mRNA to induce its degradation. The reactions for this mechanism are shown by the highlighted green cells in Table 2.

Model simplifications
There were some assumptions still in place from the original models that should be mentioned. These were justified in the supplementary data in [2] and the main text of [1]. The first model simplification is that TGFβ turnover is not included, and so it can only be added to the model in the initial conditions or via an event. This simplification was left in place as TGFβ can enter the system for a variety of reasons and incorporating these into the model would have added significant complexity. Instead, TGFβ is inactivated and activated in the model but the total pool remains constant. Any active TGFβ that enters the system is inactivated in a process that degrades Alk5 and Smad7. Alk1 can also be degraded if TGFβ is bound to an Alk1/Alk5 heterodimer. Inactive TGFβ can then be activated by integrin, which was included as a species in the model. It is known that TGFβ can be activated by many biological components such as plasmin, F-spondin or MMP-2/-9 [11,12] as well natural detergents such as lysophosphatidylethanolamine (LPE) and lysophosphatidylcholine (LPC) [13]. For simplification, Integrin is representative of all mechanisms of TGFβ activation.
The TGFβ-driven increase in MMP13 mRNA is by Runx2, which is activated by phospho-Smad1/5/8 and inactivated by phospho-Smad2/3. This is a simplification as Runx2 can be inactivated by other mechanisms [14]. However, as we are specifically interested in how TGFβ affects Runx2, having just phospho-Smads affecting its activation provides insight into how their effect on Runx2 changes with age.
The model works through the assumption that IL-1+OSM-driven MMP13 expression is a result of cFos/cJun heterodimers and cJun homodimers. Our group has recently reported that this may be a simplistic view since ATF3 appears to be important in this upregulation [15]. However, it is also evident that cFos/cJun heterodimers have a significant effect on ATF3 expression [15], such that it is reasonable to assume that blocking cFos/cJun in the way predicted in the model would have the same overall effect on MMP13 expression.
The combined model uses mass action kinetics but there are differences in the kinetic rate laws for homodimer formation between the deterministic and stochastic versions. In the deterministic model, the rate law for a homodimerisation reaction, 2X → Y is k*X 2 . However, in stochastic models, the rate is proportional to the number of ways that two molecules of type X can interact, namely X*(X-1)/2. This means that the reaction can only take place if there are at least two molecules of X. An example is for cJun dimerization which has the equation: k1*(cJun_P)*(cJun_P-1)*0.5. Only one compartment was used for the model entitled "cell". The complete lists of species and reactions in the combined model are shown in Tables S1 and S2.

Further details of parameter estimation
To parametrise the model we initially ran a global parameter scan using the genetic algorithm, followed immediately by a local algorithm Hooke and Jeeves. To fit MMP13 data we used the percentage of repression in the IL-1+OSM+TGFβ samples, compared to IL-1+OSM alone. This was accomplished by simulating the model with only IL-1+OSM active, then using the generated data to determine the level of repression necessary at each time point. For example, at 24 hours the particle number for MMP13 with simulated IL-1+OSM treatment was 109. We knew from the data that when TGFβ was present this should be reduced by 42%, so at 24 hours the value was changed to 63. This was repeated for all known time points and the IL-1+OSM+TGFβ model was fitted to those data.