ECHO, the executable CHOndrocyte: A computational model to study articular chondrocytes in health and disease

Computational modeling can be used to investigate complex signaling networks in biology. However, most modeling tools are not suitable for molecular cell biologists with little background in mathematics. We have built a visual-based modeling tool for the investigation of dynamic networks. Here, we describe the development of computational models of cartilage development and osteoarthritis, in which a panel of relevant signaling pathways are integrated. In silico experiments give insight in the role of each of the pathway components and reveal which perturbations may deregulate the basal healthy state of cells and tissues. We used a previously developed computational modeling tool Analysis of Networks with Interactive Modeling (ANIMO) to generate an activity network integrating 7 signal transduction pathways resulting in a network containing over 50 nodes and 200 interactions. We performed in silico experiments to characterize molecular mechanisms of cell fate decisions. The model was used to mimic biological scenarios during cell differentiation using RNA-sequencing data of a variety of stem cell sources as input. In a case-study, we wet-lab-tested the model-derived hypothesis that expression of DKK1 (Dickkopf-1) and FRZB (Frizzled related protein, WNT antagonists) and GREM1 (gremlin 1, BMP antagonist) prevents IL1 β (Interleukin 1 beta)-induced MMP (matrix metalloproteinase) expression, thereby preventing cartilage degeneration, at least in the short term. We found that a combination of DKK1, FRZB and GREM1 may play a role in modulating the effects of IL1 β induced inflammation in human primary chondrocytes.

observed [10,15,16]. In addition, articular cartilage conditioned medium is able to inhibit hypertrophic differentiation in growth plate cartilage and chondrogenically differentiated MSCs [12][13][14]. We have previously found that this is at least partially due to the presence of the natural BMP and WNT antagonists GREM1, DKK1 and FRZB, that are prominent in healthy articular cartilage, but not in OA [2,17].
Chondrogenic differentiation and hypertrophy are directly and tightly regulated by the activity of two main transcription factors. SRY-Box 9 (SOX9) is the master transcription factor for chondrogenic development and a key inhibitor of hypertrophic differentiation [18]. RUNX2 (Runt-related transcription factor 2, also known as core-binding factor subunit alpha-1 (CBF-alpha-1)) is a transcription factor that facilitates hypertrophic differentiation that occurs in the growth plate, and a key factor for osteoblastogenesis during subsequent bone formation [19,20]. The exact activity of these factors seems key in determining whether permanent articular or transient growth plate cartilage is formed. However, the complexity of the signaling network determining the activity of SOX9 or RUNX2 prevents a thorough understanding of the mechanisms underlying the formation of transient or permanent cartilage (Fig. 1b). This lack of understanding hampers the development of successful therapies for OA.

ANIMO in biology
Signaling networks are traditionally represented as static graphs. However, in the past years it has become obvious that the temporal and spatial information in these networks confers important dynamic behavior. As static networks do not allow quick modifications to test hypotheses or to include novel findings, a more widespread use of interactive exploration of biological networks and their dynamics could cause a paradigm shift in our understanding of biological networks.
Computer models can be used to study the dynamics of biological processes in the cell. As such, computational models enhance understanding of the biological systems and can be used to store data, interpret experimental data, help in the design of new experiments and even to identify possible diagnostic biomarkers or therapeutic targets. There has been a rise of systems biology in the field of regenerative medicine and tissue engineering [21][22][23][24][25][26]. It is important to note that the main goal of a computational model of cellular pathways is not to describe the biological system in the greatest detail, but to learn about the overall functional outcomes of the cell in response to diverse biological inputs in this model. In this sense, the model merely requires a reasonable description of the minimal number of pathways and quantitative parameters that accurately predict its molecular behavior and generates plausible outcomes that are testable by experimental observation.
We developed a visual-based modeling tool, ANIMO (Analysis of Networks with interactive Modeling), for the biological community that have no prior modeling experience and/or a solid mathematical background. ANIMO is a powerful tool to formalize knowledge on molecular interactions [27,28]. This formalization entails giving a precise mathematical (formal) description of molecular states and of interactions between molecules. Such a model can be simulated, thereby in silico mimicking the processes that take place in the cell.
ANIMO is a computational modeling tool that enables executable modeling of network dynamics in order to mimic biological phenomena in silico. ANIMO is based on the intuitive graphic interface offered by Cytoscape [29]. In addition, ANIMO has the ability to predict biological responses, both by manually testing hypotheses, as well as by using the model-checking capabilities offered by the underlying formal tool UPPAAL [30][31][32].
In this manuscript, we validated an articular cartilage model (generated as described in an accompanying MethodsX article) by performing straight forward experiments, such as changing the activities of the 7 cytokine/growth factor ligands and performing in silico knock-out and overexpression experiments and comparing this to previous data that was not used in the model-generation. We then used the model to investigate if DKK1, FRZB and GREM1 play a role in preventing IL1β induced cartilage degradation, which was validated in primary human chondrocytes.

Modeling articular cartilage using ANIMO
Growth plate cartilage and articular cartilage share a common lineage in development (reviewed in [33,34]). During development, the master transcription factor for cartilage, SOX9 is active. In endochondral bone formation, a process normally found in the growth plate, the transcription factor RUNX2 becomes active (Fig. 1A). However, the exact way signals are integrated to regulate activity of RUNX2 and/or SOX9 is still a black box (Fig. 1B). The transcription factors SOX9 and RUNX2 are key factors in cartilage development and homeostasis by regulating the gene expression of chondrogenic and cartilage hypertrophic genes, respectively (reviewed in [35]). During cartilage development, SOX9 complexes with SOX5 and SOX6 to regulate expression of COL2A1 (collagen type II), COL1A1 (collagen type I), ACAN (aggrecan) and PtHrP, thereby preventing hypertrophic differentiation [36]. In contrast, RUNX2, which can complex with RUNX3, is essential for the hypertrophic differentiation of chondrocytes and upregulates expression of COL10A1 (collagen type X) and IHH [37].
To describe the processes regulating cell-fate in chondrocytes, we developed a computational model, which we named the Executable CHOndrocyte, or ECHO. ECHO consists of seven signaling pathways known to be important in cartilage development and maintenance: Insulin growth factor (IGF), Parathyroid Hormone Related Protein (PTHrP), Bone Morphogenetic protein (BMP), Fibroblast growth factor (FGF), transforming growth factor β (TGFβ), Wingless-type MMTV integration site family (WNT), Indian Hedgehog (IHH). ECHO was initially based on a Boolean model of the growth plate [38,39], which we converted to an ANIMO model as described in the accompanying MethodsX manuscript. This resulted in a growth plate model (Fig. 1C).
The growth plate model (GP, Fig. 1C) was adapted to an articular cartilage model (AC) based on global gene expression microarrays of growth plate and articular cartilage, see MethodsX Section 1.3 and [2]. Specifically, we added the WNT antagonists DKK1 and FRZB and the BMP antagonist GREM1 downstream of SOX9 and included known cross-talk of these factors as shown in Fig. 1D. In addition, we found that the expression of 4 genes, whose proteins were already represented in our model: ERK2, p38γ, GSK3β and Smad3 was significantly changed. We therefore multiplied downstream interaction parameters with the relative expression levels to take differences between tissues into account. For example, the microarray data show that an articular chondrocyte expresses less p38 than a growth plate cell ( Fig. X6: AC/GP = 0.68); as a consequence, we multiplied the parameters of the interactions downstream of p38 by 0.68 (Fig. 1E).
Because of their key roles in cartilage SOX9 and RUNX2 are implemented in ECHO as outputs of the cell fate. Since SOX9 activity is necessary for the permanent state of (articular) cartilage, and RUNX2 is essential for cartilage hypertrophy, we assume the SOX9 + state to be an articular cartilage phenotype while the RUNX2 + state is a hypertrophic phenotype in our models.

Changes in activities of receptors ligands and their antagonists influence the outcome of the network
Not all nodes in the network influence the outcome of the stable states equally. To determine the influence of nodes on the outcome of the network, we used Monte Carlo simulations as described in MethodsX 1.4, in which all nodes were initially assigned a random, uniformly distributed activity level over the entire range of possible values (i.e. the interval [0, 100]). Each initialized model was then simulated until a stable state was reached. Fig. 2A shows the 15 nodes that have the most prominent effect on cell fate. For detailed description of the experiment and the roles of these nodes in cartilage, we refer to the MethodsX article.
The top nodes from Fig. 2A are either ligands/antagonists, that can be pro-SOX9, pro-RUNX2 or both, depending on their activity and the cross-talk with other pathways (Fig. 2B), or they can be intracellular signaling molecules that have a great influence on the network. The extracellular growth factors that activate cellular signaling pathways, together with antagonists GREM1 and FRZB (Fig. 2B), are expected to have a great influence on the network, since the presence and activity of these nodes initiate responses downstream of these signaling pathways, resulting in changes in cell fate. DKK1 is also shown, but plays a lesser role in the model, primarily because of redundancy of its effect with respect to FRZB, in our model. Based on the results in Fig. 2A, FGF and WNT can be classified as pro-RUNX2, whereas PTHrP, TGFβ and IGF1 have a pro-SOX9 effect. BMP and IHH have a more complex, mixed downstream effect and activate both RUNX2 and SOX9. In addition, we can find that a few single nodes by itself play an important role in determining cell fate: SMAD7, MSX2, PPR (PTHrP receptor), DLX5, and RUNX2 and SOX9.
We validated the likelihood of these findings with literature before proceeding with the next experiments. DLX5, distal-less homeobox 5, is a transcription factor that is not bone or cartilage specific but is expressed in the early stages of bone formation and mediates RUNX2 expression [40]. DLX5 is a positive regulator of chondrocyte maturation (read: in GP chondrocytes) during endochondral ossification [41]. MSX2 (Msh Homeobox 2) is a transcription factor important for bone formation, and controls IHH expression to stimulate chondrocyte maturation [42]. Both MSX2 and DLX5 are regulated by BMP signaling (reviewed in [43,44]). As such, it is clear why these factors play such an important role in regulating the bi-stable character of our models.
The second subnetwork of nodes is a small network of transcription factors that gives insight in the origin of ECHO's bi-stable behaviour (Fig. 2C). In this network, SMAD3 is shown as the key signal transduction molecule that promotes SOX9 activity, and the powerful role of MSX2 is stressed by showing its role as an inhibitor of WNT expression. It can be concluded that these nodes are key to regulating the cartilage phenotype in our model.

The effects of extracellular ligands on cell fate are dose-dependent and different in the growth plate and articular cartilage models
ECHO summarizes the effects of multiple pathways, allowing us to study their cross-talk. However, the activation of each pathway depends also on the dose of the related extracellular ligand. To investigate in which measure the dose of ligands influences the cell fate, we performed a series of in silico dosage experiments. We investigated the dose dependency for all extracellular ligands (Fig. 3). For example, to better understand the effects of BMP, its activity was artificially fixed at different levels in the range 0-100, in steps of 10 levels. All other nodes were randomly initialized and the resulting cell fate distribution was analyzed over the course of 10,000 simulations for each BMP initial Top 15 nodes with the most prominent effect on the SOX9 and RUNX2 stable states in both GP and AC ECHO-models. A. The values in the "SOX9 + " and "RUNX2 + " columns give the median initial activity for the simulations that resulted in SOX9 + and RUNX2 + cell fates, respectively. The simulation results were divided into categories corresponding to the resulting cell fate (SOX9 + and RUNX2 + ), and the median initial activities of all nodes were computed. The median initial activity of a node with no influence on a cell fate is 50, which would correspond to the node's value being randomly chosen on the interval [0,100]. Conversely, if a node activity needs to be above/below a certain threshold in order to obtain the desired fate, the corresponding median value deviates from 50. The higher a node's absolute deviation from 50, the more a cell fate depends on that node being active/inactive. The values for the "Influence" column were computed as the average absolute deviation from 50 of each node. B, C. Identifying the subnetworks involving the top nodes allowed us to highlight two groups of nodes that have the largest influence on cell fate. B. All extracellular signaling molecules in ECHO (central row) are major factors in the determination of cell fate in both the GP and AC model. IHH exerts its effect by activating Gli2, which in turn plays a role in the regulation of WNT, BMP, PTHrP and TGFβ. DKK, FRZB and GREM are antagonists of WNT signaling and BMP signaling respectively, and together favor the SOX9 + state. DKK1, FRZB and GREM1 are only present in the AC model C. SOX9, RUNX2, DLX5 and MSX2 are at the core of the bi-stable behavior of ECHO. SOX9 and RUNX2 stimulate their own activity and inhibit each other. DLX5 and MSX2 also inhibit each other and have opposite effects on RUNX2. SMAD3 is the main effector of TGFβ and acts as a strong pro-SOX9 node in this subnetwork.
value. This experiment, for the GP model, shows that there is an optimum dosage of BMP to reach the SOX9 + fate, with a peak at BMP = 50. At higher BMP levels, there is an increasing tendency towards RUNX2 + . This dual role of BMP2 has been partially shown in other studies. For example, knock-out of BMP2 has been shown to be essential for both chondrocyte proliferation and chondrocyte hypertrophy in the growth plate [45]. In addition, overexpression of BMP2 results in osteophyte formation in an experimental model of OA [46]. These data indicate that BMP2 at low levels is necessary for chondrogenesis, but high BMP2 expression leads to cartilage hypertrophy and OA. This is nicely captured in our models. However, our models indicate that the levels at which BMP2 functions may be different between articular cartilage and growth plate cartilage, and that this may be regulated by the activity of the IHH/Gli2 pathway.
It is interesting to note that some ligands have an effect that is dosage-dependent: for example, the highest chance of obtaining a SOX9 + fate is achieved when the initial dosage of BMP is at intermediate levels.
In most other cases, the percentage of cell fates linearly depend on the dosage of the ligands: for example, higher initial activities of PTHrP always lead to more SOX9 + , while lower activity of TGFβ give more RUNX2 + . Higher activity of TGFβ and of IHH elicits a strong power to induce the SOX9 + fate, most notably in the AC model. Only at high IGF1 activity did we observe a stable SOX9 + state, while lower activity results in a Null state. Interestingly, most labs use ITS-Premix (insulintransferrin-selenous acid-Premix) as a supplement in the chondrogenic differentiation medium. The insulin content is 100-fold above physiologic concentrations and, as compared to IGF1, has a two order of magnitude lower affinity for the IGF1R. This possibly results in nonspecific activation of the IGF1 signaling pathway [47]. It has been shown that IGF1 activation is additive to TGFβ in chondrogenically differentiating MSCs [47]. The authors showed that IGF1 modulated chondrogenesis of MSCs by stimulating proliferation, inducing expression of chondrogenic markers and regulating apoptosis. This stimulatory role of IGF1 in cartilage is reflected in our model.

Combinatorial dosage experiments give insight in the effect of competing extracellular ligands
To study the effects of cross-talk between the pathways in ECHO, we selected two extracellular ligands with opposite effects and performed a series of in silico combinatorial dosage experiments. The results for a selection of this experiment are shown in Fig. 4. Other tested  combinations did not show a clear difference in response between the GP and AC model.
Again, we note that, particularly in the GP model, intermediate levels of BMP give the best chances of obtaining a SOX9 + cell fate. When the level of BMP is too low, the intrinsic predisposition of ECHO models to reach a Null fate can sometimes be counterbalanced by high levels of TGFβ, IHH or IGF1, especially in the AC model (Fig. 4, left panel). For comparison, when the initial BMP activity is higher, the chance to encounter RUNX2 + cell fate rises. In all cases, the differences between GP and AC models are significant, with a much higher likelihood to reach the SOX9 + state in the AC model. The comparison between WNT and BMP again highlights the importance of intermediate BMP levels to reach SOX9 positivity.
Clear opposite effects are shown by WNT and TGFβ: only when one is low and the other high can the model reach a non-Null fate, with WNT forcing equilibria towards RUNX2 + states and TGFβ favoring SOX9 + states. In general, we note that only relatively high levels of WNT can lead the model towards a RUNX2 + cell fate in the presence of pro-SOX9 ligands. We also note that different pro-SOX9 ligands can counteract the effect of high levels of BMP to different degrees (Fig. 4). IGF1 is less powerful than IHH, which in turn is less powerful than TGFβ in preventing the RUNX2 + fate resulting from high levels of BMP in the GP model. In line with the previous results, the AC model shows an increased dominance of the SOX9 + cell fate.

Changing cell fate through perturbations of extracellular ligands
There is evidence that hypertrophy-like changes in chondrocyte terminal differentiation are also found in OA (reviewed in [16]). Therefore, we hypothesized that the switch from a SOX9 + to a RUNX2 + cell state would be easier than a switch from RUNX2 + to a SOX9 + stable state, and that the GP model would be more prone to the RUNX2 state than the AC model. We performed experiments to test transitions of cell states in silico to aid in our understanding of the reversibility of growth plate and articular cartilage tissue development. These studies essentially model cell programming (i.e., trans-differentiation) that could permit deliberate switching of cellular states in tissue engineering applications.
For this experiment, the starting situation was either the RUNX2 + or the SOX9 + stable state. Each combination of two extracellular ligands was perturbed from their steady state values by changing their activity levels to fixed levels over the complete range [0,100]. The model behavior is then evaluated until a stable state is reached and the resulting state is recorded (Fig. 5). The figure shows the nonlinear range of activity levels that results in each cell fate. Perturbations of one or two ligands reprograms a RUNX2 + state to switch into a SOX9 + state in both the GP and the AC models. Note that in the AC model it is possible to have a switch from RUNX2 + to SOX9 + using high levels of TGFβ alone: this can be seen in the first row of the table, where the effect of TGFβ dosage is plotted in vertical stripes (Fig. 5A). Some combinations lead to a SOX9 + state only with extreme values: for example, the combination BMP and TGFβ leads to a SOX9 + state only with 0 BMP and 100 TGFβ, indicated in row 4 of the table as BMP (− ) and TGFb (+).
Perturbations of one or two ligands can make a SOX9 + state switch to a RUNX2 + state in both GP and AC models (Fig. 4B). Note that BMP and WNT can be used alone (first two rows) to cause a switch from SOX9 + to RUNX2 + if their activity is set at high enough levels in the GP model, while this does not hold in the AC model. Also, starting from the stable SOX9 + state in row 4, increasing BMP activity (X-axis) with decreasing PTHrP activity (Y-axis) causes a switch to the RUNX2 + state, indicated by BMP (+) and PTHrP (− ) (Fig. 5B). In this example, we also note that intermediate values of BMP favor the SOX9 + state. It can be seen that activation of the WNT or FGF signaling pathways results in a switch from SOX9 + to RUNX2 + (Fig. 5A). This can be explained as both WNT and FGF signaling have been related to induction of hypertrophy in cartilage (reviewed in [10]).
In general, the number of ways to reprogram cells to achieve a SOX9 + ➔ RUNX2 + switch decreases from the GP model to the AC model, while the reverse switch becomes more accessible (cf. the number of "No switch" between the two models). These observations are in line with our hypothesis and indicates that the AC model describes the fate of a stable articular chondrocyte. In particular, the changes from the GP to the AC model increase SOX9 + stability to the expense of RUNX2 + . Moreover, the range of activity levels that cause a switch to SOX9 + is broader in the AC model than in the GP model.

Modeling cartilage development: not all stem cells have the same chondrogenic capacity
Our previous models captured events representing growth plate and articular cartilage. To enable research into potential combinatorial treatments for OA, it is necessary to capture the different steps of chondrogenesis in our models. To investigate the chondrogenic capacity of various stem cells, we performed in silico stem cell differentiation experiments in ECHO based on RNASeq data from embryonic stem cells, mesenchymal stem cells and articular chondrocytes. From a tissue engineering perspective, it is necessary to define the exact circumstances that dictate cell fate decision. We initialized the top 15 most important nodes (from Fig. 2) in ECHO with different configurations, based on RNASeq data of articular cartilage, bone marrow Mesenchymal Stem Cells (BMSC), Adipose-derived mesenchymal stem cells (AMSC) or human embryonic stem cells (ESC) (data sets: [48][49][50], MethodsX table X5), then used the model to simulate the reaction of the cell to a differentiation stimulus. The stimuli were chosen based on BMSC chondrogenic differentiation protocols commonly used in our lab [8,9,51]. At the end of each experiment, the resulting cell state was recorded. The differentiation stimuli used for the different cell types and the resulting cell fate percentages after the simulations are shown in Fig. 6.
We observe differences in chondrogenic capacity between the different cell types, as also observed in literature. In silico chondrogenesis of the AC model yields efficiencies that are between articular chondrocytes and BMSCs under the same stimuli, indicating that our AC model may be more representative for an early articular chondrocyte rather than a mature articular chondrocyte. Interestingly, the model predicts that it is possible to generate SOX9 + cells from ESC by using a 'standard' MSC differentiating medium containing 40 TGF and 80 BMP. However, we know from biological experiments that multiple stages of differentiation, each containing different stimuli, are necessary to obtain cartilage [52,53].
Another observation from our in silico experiments is that BMSC and AMSC respond similarly to the stimuli, with a slightly better chondrogenic capacity of the AMSC. This corresponds to previous reports that state that BMSC and AMSC are equally effective in cartilage formation in co-culture conditions with articular chondrocytes [51,54], and that AMSC may have a slightly better chondrogenic capacity, at least in OA conditions [55]. When BMP and WNT are both active in our model, RUNX2 + or Null states are reached in all models, except in the AC model, where DKK1, FRZB and GREM1 function as BMP and WNT antagonists in at least a percentage of the simulations. This also corresponds to experimental data, showing that both WNT and BMP activity leads to cartilage hypertrophy and RUNX2 activation [10,56].

Case-study: using ECHO to predict the effect of IL1β on the activity of SOX9 and RUNX2 in primary human chondrocytes
Our studies on switching cell fates from SOX9 + to RUNX2 + identified the GP model, especially in the RUNX2 + state, as the model closely describing a hypertrophic phenotype. However, a switch from a SOX9 + state to a RUNX2 + state in the AC model can be obtained by changing inputs, leading to an (in silico) hypertrophic state. As described before, hypertrophy is seen in some patients with OA, and this hypertrophy is thought to contribute to the development of osteophytes [15]. One of (caption on next page) the biggest challenges in validating our model generated hypotheses, is that work on human cartilage development and OA is performed in vitro in primary cells. These primary cells are often isolated from patients undergoing total joint replacement therapy for OA management. We have observed that the cartilage in the OA joints is not homogeneously damaged throughout the joint (non-published observations and [57]). Much of the research in primary chondrocytes with regards to the response of cells to external stimuli has been performed in cells isolated from an OA joint. There is evidence that the overall inflammation status affects all tissues in the joint, including cartilage [58].
We have previously shown that there is a difference in the amount of IL1β expression in healthy looking cartilage from an OA joint and its OA paired sample [59]. There is a positive feedback loop between IL1β and its gene, suggesting that once the cells are exposed to IL1β, joint inflammation is imminent. We hypothesized that OA and healthy chondrocytes respond differently to exposure to IL1β signaling, and in particular that healthy cells are somehow protected against inflammatory signals such as IL1β.
IL1β signaling was not initially included in our GP and AC models, so we started by adding this pathway to our model (Fig. 7A, Supplemental  Fig. 1). We have previously generated a computational model of IL1β signaling in cartilage [59]. Using this model, we identified that WNT activity is upregulated in human OA chondrocytes by reducing the expression of the WNT antagonists DKK1 and FRZB via iNOS [59].
With the addition of the IL1b pathway to the AC model, we want to investigate the influence of that important pathway on healthy AC cells. Because of the robustness of the healthy AC cell (represented by the AC model in its SOX9 + state), simply setting IL1β to its maximum activity with no inhibition (i.e. setting the activity of the IL1b node to 100% and disabling the IL1b -| IL1b interaction) would not be enough to switch from SOX9 + to RUNX2 + . In Fig. 7A, B we show the result of this attempt: 7A is the initial state of the simulation, and 7B is the final state. The numerical activity values for some nodes in the network in the final state of the simulation are shown also in the first row in Fig. 8B.
Based on the finding that IL1β induces cartilage degradation [60], we assumed that addition of IL1β would reinforce the RUNX2 + phenotype and weaken the SOX9 + one. We tested this by setting the initial activity of IL1β to 100%, disabling its self-inhibition (to keep it artificially high) and performing two simulation runs in the AC version of ECHO, one using the RUNX2 + initial activity configuration and the other using SOX9 + as initial condition.
The activity of the RUNX2 and SOX9 nodes does not change in either simulation, but the simulation starting from SOX9 + ends in a state where some nodes have significantly changed their activities (see Fig. 7B). Interestingly, production of RUNX2 is observed (node "Runx2 prot" activity at 40%); however, the transcription factor seems to lack the needed post-translational modifications (node "Runx2 PTM" activity at 0), which makes the overall state of the network "quasi-SOX9 + ".

Validation of ECHO predictions with wet-lab experiments
In ECHO, we initially try to investigate general mechanisms, which can be highly influenced by donor variation, disease state [57], and time of cells in culture [4]. Testing our hypotheses without the use of computational models requires (too) many cells, which is especially hard when using primary cells and tissues. As described in the accompanying MethodsX paper, ECHO was generated based on literature data of cartilage development in a variety of species. In addition, ECHO was built based on a variety of inputs: gene expression, data on protein interaction and signal transduction pathways, as well as protein activity data are incorporated into this large network. Most of what we know about cell's responses comes from gene expression (qPCR) data. To validate ECHO predictions, we need to look at the influence of IL1β on the activity of SOX9 and RUNX.

IL1β affects SOX9 and RUNX2 activity
In ECHO, addition of IL1β resulted in a quasi-SOX9 + network state, where the steady states did not change from SOX9 + to RUNX2 + , but the downstream effects changed. For validation of our model-predictions, we want to both investigate the effect of network changes on the activity of SOX9 and RUNX2, as well as on downstream effects. To investigate the direct effect of IL1β on the activity of SOX9 and RUNX2, we performed transcription factorfluorescent recovery after photobleaching (TF-FRAP) [61]. We have previously shown that TF-FRAP is a very fast method that shows the direct effect of external influences on transcription factors [61]. We have shown that the mobility of SOX9 correlates to its transcriptional activity. The higher protein mobility correlated to a loss of DNA binding and a lower gene expression of its target genes COL2A1 and ACAN [61]. TF-FRAP is a single cell method, enabling us to investigate the range of responses in the heterogeneous  1 (0, 1, 2, … 100), while all other nodes were initialized as in a RUNX2 + state (A), or a SOX9 + state (B). After one simulation, the resulting stable state was recorded. The X-axis of all graphs corresponds to the first perturbed node, and the Y-axis to the second. For example, in B4, the X-axis indicates BMP, the Y-axis PTHrP activities. The '+' sign means addition of the factor is necessary to obtain a switch in cell fate, '-'means absence of this factor is a requirement to obtain a switch. A. Perturbations of one or two ligands that can make a RUNX2 + state switch to a SOX9 + state in GP and AC. B. Perturbations of one or two ligands that can make a SOX9 + state switch to a RUNX2 + state in the GP and AC models.   population. Parameters give information on the mobility of the protein and are measured in terms of half-time of recovery (t ½, indicative of the diffusion rate of the protein) and immobile fraction (IF, indicative of the amount of protein that is tightly bound to DNA), see methods and [61].
We used primary chondrocytes from donors undergoing total joint replacement therapy due to OA. From these donors, chondrocytes were isolated from the affected area (OA chondrocytes) as well as from the macroscopically healthy site (Healthy looking/HL chondrocytes). The cells were transfected with either SOX9-mGFP or RUNX2-eGFP (Fig. 7C) and treated with 10 ng/ml IL1β for 20 min before performing TF-FRAP (Fig. 7D,E,F). The changes in mobility are compared to untreated control cells (Supplemental Fig. 1). SOX9 in OA cells shows a lower Immobile Fraction (45%) than in the HL cells (49%), indicating that SOX9 is more active in the HL cells. The HL cells display a larger change in Immobile Fraction in the presence of IL1β than the OA cells, albeit not significant due to the large heterogeneity. The recovery t ½ of SOX9 after IL1 treatment in both HL (12.8 s vs 13.4 s in the control) and OA (12.6 s vs 13.6 s in the control) cells is significantly lower, indicating a lower residence time of SOX9 on the DNA. In both HL and OA chondrocytes cells we observe little binding of RUNX2, as indicated by the low immobile fraction (<30%). There is a significant increase in RUNX2 IF in chondrocytes treated with IL1β as compared to the control, although the IF is slightly increased in the OA cells after IL1 treatment, and the recovery t ½ is higher in the HL cells after treatment, indicating a longer residence time of RUNX2 on the DNA (Fig. 7F, Supplemental Fig. 1).
To show that treatment of cells can change the SOX9 mobility, we performed a positive control where we treated the cells with 100 ng BMP7 and performed SOX9-FRAP. This resulted in a significant increase in the immobile fraction (52% for HL, 50% for OA) and in the half-time of recovery (14.4 s for HL, 15.5 s for OA), indicating that SOX9 is more active (Supplemental Fig. 1). This corresponds to the ECHO prediction on the influence of BMP on SOX9 in Fig. 3. Together, these data indicate that the effects of IL1β treatment on SOX9 and RUNX2 activity result in only small effects at the level of transcription factor activity. This corresponds to what we observe in ECHO: a quasi SOX9 + state is obtained.

IL1β treatment reduces expression of SOX9 targets and increases expression of MMP
To test the effects of IL1β treatment on gene expression of SOX9 target genes and MMP expression, HL and OA cells were treated with 10 ng/ml IL1β for 24 h and the expression of COL2A1, ACAN, and MMP1, MMP3 and MMP13 was tested and compared to untreated controls (Fig. 7G). Expression of COL2A1 and ACAN was decreased in OA cells when compared to the healthy cells (Fig. 7G). The expression of MMP1, 3 and 13 is significantly increased in both cell types after IL1β treatment (Figure7G). Moreover, the expression of all MMPs was significantly higher in OA cells treated with IL1β as compared to the expression in healthy cells. The opposite was true for ACAN and COL2A1 expression, which was significantly decreased in the OA cells, but not in the healthy cells. It appeared that loss of expression of the matrix proteins ACAN and COL2A1 is partially prevented in healthy cells, whereas in OA cells the presence of IL1β results in an even higher loss of matrix.

DKK1, FRZB and GREM1 may play a role in protecting healthy cartilage from IL1β
To elucidate a possible mechanism by which healthy cells are, at least partially, protected against IL1β, we turned to our AC-model + IL1 ECHO model. We repeated the simulation starting from SOX9 + disabling each one of DKK1, FRZB and GREM1 in turn: this resulted in a RUNX2 + final state only when FRZB was disabled, while the two other cases resulted in "quasi-SOX9 + ", with a relatively high amount of RUNX2protein activity. This did not lead to a shift in cell state from SOX9 + to RUNX2 + . In addition, we observed an increase in MMP activity in our model, since the expression of MMP was directly downstream of iNOS [59]. We note that disabling GREM1 allowed the BMP pathway to become active, but this was not enough to change the global network state. Indeed, these nodes were identified as the most significantly different in expression between GP and AC cartilage and the expression of DKK1 and FRZB is lost in OA [2,17,57,62]. If we reduce the activities of FRZB, and to a lesser extent DKK1 and GREM1, in the AC model when it is in the SOX9 + state, and we then add IL1β, we observe a large increase in MMP expression and a rescue of the IL1β induced loss in collagen 2 expression (Fig. 8A,B).
Since ECHO predicts that a variety of combinations of DKK1, FRZB and GREM1 affect Collagen 2 and MMP expression, we only tested the combination of DKK1, FRZB and GREM1 in presence or absence of IL1β in OA cells by TF-FRAP. All treatment conditions lowered the immobile fraction of SOX9, but only addition of DKK1, FRZB and GREM1 to the cells in presence of IL1β slightly, but significantly, lowered the immobile fraction of SOX9 compared to the control, but not compared to the other conditions. The t ½ of A2 was not significantly changed (Supplemental Fig. 2), indicating that the residence time of SOX9 on the DNA is unaltered. Together, these data show that a possible influence of IL1β + DKK1, FRZB, and GREM1 is not detectable at the level of SOX9 activity.
All treatments increased the IF of RUNX2, but only IL1β treatment resulted in a significant increase in IF, indicating that in OA cells RUNX2 becomes more active in the presence of IL1β. The residence time on the DNA for RUNX2 was decreased in presence of DKK1, FRZB and GREM1, as indicated by a significant decrease in t ½ (Supplemental Fig. 2). Treatment with IL1 increased the recovery half-time, while co-treatment of IL1β with DKK1, FRZB and GREM1 brought the halftime back to control values (Supplemental Fig. 2). This indicates that IL1β slightly increases RUNX2 activity in OA chondrocytes and that addition of DKK1, FRZB and GREM1 may be able to modulate this effect. However, the heterogeneity of the cell response is too large to get a clear image of how DKK1, FRZB and GREM1 may modulate SOX9 and RUNX2 at the level of transcription factor activity using FRAP.

OA cells respond differently to IL1β stimulation than healthy chondrocytes due to the absence of DKK1, FRZB and GREM1
To investigate the possible role of DKK1, FRZB and GREM1 in modulating the response of cells to IL1β, we performed qPCR in HL and OA chondrocytes, which we supplemented with recombinant DKK1, FRZB, and GREM1, in presence or absence of IL1β. We observed a slightly higher expression of the SOX9 target genes COL2A1 and ACAN in the presence of IL1, DKK1, FRZB and GREM1 as compared to the cells treated with IL1β alone. This effect was seen in all 3 donors tested, although the effect is most prominent in the OA donors (Fig. 8C). Due to the large donor variation and the limited number of donors tested, we did not perform statistical analysis. We observed that the combination of DKK1, FRZB and GREM1 was not able to prevent IL1β induced MMP expression (Fig. 8D). These data are in accordance with the predictions of ECHO, where 100% activity of MMP expression and lower expression of COL2A1 are predicted. In addition, we observe a reduction in the number of apoptotic cells in the presence of IL1β with DKK1, FRZB and GREM1, as compared to IL1β alone (Fig. 8E). These data confirm our ECHO generated hypothesis that DKK1, FRZB and GREM1 play a small role in preventing IL1β induced cartilage damage in healthy cartilage. This has not been shown before.
Together, this data shows that ECHO, our executable chondrocyte model, is excellently suited for storing data. It enables the dynamic and time-dependent visualisation of molecular mechanisms involved in the gradual transformation of a healthy chondrocyte in an osteoarthritic cell. Executable models can furthermore be used to predict new mechanisms for cellular responses in cartilage health and disease. This information can be used for the rational design of the most promising new wet-lab experiments with respect to the desired outcome. We believe that executable computer models of molecular mechanisms of disease may become a valuable asset in the development of new treatments.

Discussion
We successfully used ANIMO to generate models of growth plate and articular cartilage based on an earlier model of the growth plate [39,63] and microarray data of articular cartilage [2]. These models were used to obtain insight into the signaling network regulating cartilage development and OA. We found a possible mechanism by which healthy cartilage is protected against low levels of inflammatory signals, via suppression of the WNT and BMP pathways.
Both the BMP and the WNT pathways are implicated in cartilage hypertrophy and OA (reviewed by [10,64]). In ECHO, BMP is a node with much influence on the determination of cell fate in the AC and GP models. Much research in in vitro and in animal models has shown that BMP2 is involved in cartilage development and endochondral ossification by regulating both SOX9 and RUNX2 [43]. In addition, SOX9, ACAN, collagen type II were reduced at the mRNA levels in BMP2/4 double KO embryos or BMP2 conditional KO mice, while RUNX2 expression was reduced at the protein level in the proliferating and prehypertrophic areas BMP2/4 double KO embryos [45]. Moreover, BMP2 induced RUNX2 expression at both transcriptional and posttranscriptional levels through inhibition of CDK4 and prevented RUNX2 degradation [45].
ECHO shows that the BMP concentration has to be tightly controlled in order to induce a SOX9 + cell state in the GP model, so that little BMP is not enough to induce any differentiation, but too much BMP leads to RUNX2 activity and thus hypertrophy and collagen 10 expression. This is verified in ECHO by both the co-activity experiment in Fig. 4 and by the perturbation experiment in Fig. 5B (SOX9 > RUNX2) in which constitutive activation of BMP in combination with the absence or KO of FGF, IGF, TGFβ or IHH never yields a SOX9 + state in both the GP and AC model. This is corroborated by the fact that different BMP antagonists, such as GREM1 and Noggin, are expressed in specific layers of the articular and growth plate cartilage [2,17,65].
Other research highlights the importance of BMP in the cartilage network: for example, it has been shown that the lack of BMPR1A leads to significant chondrodysplasia and almost eliminated the chondrocyte phenotype with decreased SOX9, collagen II, and proteoglycan expression in conditional BMPR1a KO mice [66][67][68].
The WNT signaling pathways, and specifically the WNT/beta-catenin canonical pathway, has been described to be important for both correct cartilage development and hypertrophy (reviewed in [10,64]). We have previously found the WNT antagonists DKK1 and FRZB to be important for the prevention of hypertrophy and we found that the expression of DKK1 and FRZB was reduced in OA [2,57,17,69]. In addition, we found that IL1β was able to reduce DKK1 and FRZB expression in human chondrocytes that were isolated from healthy looking cartilage regions in an OA knee via upregulation of iNOS [59]. These relatively healthy cells had a basal level of DKK1 and FRZB expression and had a reduced IL1β and MMP expression as compared to the OA cells from the same joint [59]. Prolonged exposure to inflammatory factors, such as IL1β, has been implicated in the hypertrophic differentiation via both NFκB and MAPK signaling [70,71], and our data suggest that the loss of DKK1, FRZB and GREM1 plays a role in the loss of the protective properties of healthy cells against IL1β exposure. This is, to our knowledge, a new finding.
Another important player in ECHO is IGF1. KO of IGF1 in both GP and AC models reduced the percentage of SOX9 + cell fates, while overexpression increased the percentage of SOX9 + cell fates, see Fig. X7 in the accompanying MethodsX manuscript. In contrast, the percentage of RUNX2 + cell fates were increased in the GP model when IGF1 was absent. This is not surprising as IGF or insulin are used in most standard protocols for MSC chondrogenesis [72,73].
Both GP and AC model indicate an important role of TGFβ in the development of SOX9 + state. This is supported by literature: Loss of TGFβ expression has been correlated with OA in mouse models [74]. In addition, the correlation between genetic alterations of TGFβ and OA have been reviewed in [75]. For example, it has been shown that impaired chondrocyte terminal differentiation was observed in ERK1 and ERK2 KO mice [76]. This is also observed in our model: ERK1/2 knockout prevents a transition from a SOX9 + to a RUNX2 + state, whereas overexpression induces a RUNX2 + state (see Fig. X7). In the GP model, ERK1/2 knock-out results in a transition to a Null-state, indicating the importance of ERK1/2 in the regulation of hypertrophy and endochondral ossification via RUNX2. This is corroborated by literature [76][77][78].
We determined the top 15 nodes with the most prominent effect on the SOX9 + and RUNX2 + stable states in ECHO. We showed how computational models can aid in predicting the chondrogenic differentiation capacity based on RNASeq data of only these 15 genes. In addition, these experiments predict the cellular response on a variety of external signals that are commonly used in differentiation experiments. In our top 15 most important nodes, we also identified PTHrP, IHH and GLI2/3. The roles of IHH and PTHrP in regulating bone and cartilage formation have been well established [79]. In ECHO, we find that both overexpression and KO of PTHrP or IHH leads to a transition from SOX9 to RUNX2 to Null (data not shown). Expression of SOX9 and RUNX2 as well as PTHrP was low and growth was inhibited in temporomandibular joint in IHH KO mice, indicating that IHH is indispensable for proliferation and expression of transcriptional regulators such as RUNX2 and SOX9, however, these defects were partially restored in double IHH/ Gli3 mutants, suggesting that IHH function is modulated and restricted by Gli3 and Gli3(R) [80,81]. Growth inhibition and less expression of RUNX2, RANKL (Receptor Activator Of Nuclear Factor Kappa B Ligand), and MMP13 was found in TGFα (Transforming Growth Factor Alpha) KO mice [82]. PTHrP dependent and independent effect were observed in TGFβRII KO mice [83]. The expression of collagen type X, MMP13, ADAMTS4 and ADAMTS5 (A Disintegrin-Like And Metalloprotease (Reprolysin Type) With Thrombospondin Type 1 Motif, 5 (Aggrecanase-2)) was up-regulated in TGFβRII (Transforming Growth Factor Beta Receptor 2) KO mice, indicating cell differentiation process was accelerated and matrix degradation activated. PTHrP KO mice show increased endochondral bone formation [83]. Limb explant treated with PTHrP showed inhibition of collagen type X expression and enhanced chondrocyte proliferation. This effect was potentiated in GLI2 KO limbs, while it was blocked in GLI3 KO limbs and was only partially inhibited by hedgehog ligand blockade. PTHrP negatively regulated Gli mediated transcription in cell cultures and regulated the level of the repressor form of Gli3 in a PKA dependent manner. These results show that PTHrP regulates growth plate chondrocyte proliferation and differentiation in part through the activity of Gli3, suggesting a crucial role for Gli3 in growth plate chondrocyte development [84].
From a therapeutic standpoint, it is interesting to know if and how it is possible to switch cell fate through perturbation of the network by extracellular ligands. Indeed, it was possible in ECHO to define various factors and especially combinations of factors able to switch cell fate (Fig. 5). It is experimentally challenging to test these predictions in human primary cells. This is especially the case for the switch of a RUNX2 + to SOX9 + cell fate, since epigenetic regulation has likely occurred in the process of osteoarthritis (reviewed in [85]). The resulting methylation of cartilage specific genes, such as SOX9, will therefore prevent the actual switch from a RUNX2 + to a SOX9 + phenotype in OA cells. However, this strengthens the argument of using computational modeling, since it allows us to simulate osteoarthritis development, and as such provides insight into the molecular boundaries that define therapeutic efficacy. In addition, the combination of factors tested in our in silico experiments have individually been described to have a role in cartilage development and in OA. However, the combined effects of these factors have not yet been conclusively shown. It is therefore likely that when designing therapies for treatment of cartilage defects multiple factors will have to be targeted in order to get the desired response. Computational modeling is very well suited to investigate combinatorial treatments for disease.
Our case-study shows that ECHO can be used to investigate cellular responses to changes in network activities, both by manipulating nodes in so-called knock-down and overexpression experiments, as well by changing activities of extracellular ligands. In our case-study, we prioritized the experiments to be performed, thereby limiting the number of experimental conditions from 9 to 4 conditions.
Our study shows a possible role of DKK1, FRZB and GREM1 in protecting healthy cartilage from the effects of inflammation and indicates that this protective effect is only partially dependent on the transcriptional activity of SOX9 and RUNX2. However, the exact mechanism by which this happens in healthy chondrocytes needs to be investigated further both by adapting ECHO and by extensive testing in many completely healthy chondrocytes from various donors. This is outside the scope of this paper.
In this manuscript, we present the application of GP and AC models, with adaptation to an model in which IL1β is active, which may represent OA-like states. In these models we are able to investigate the network response to changes in the cellular environment or changes in the intracellular signaling network. We used these models to obtain insight into the mechanisms of stem cell differentiation and cell fate determination in disease. We validated our model topology and activity by using RNASeq data as well as by performing FRAP and qPCR analysis. The model hypothesis that DKK1, FRZB and GREM1 could function as the gate-keepers of cartilage homeostasis has been validated in a wet-lab experiment that showed that the IL1β induced loss of Collagen 2 and Aggrecan expression was partially prevented by, and apoptosis was reduced in the presence of DKK1, FRZB and GREM1 in human chondrocytes.

Computational cluster
The in silico experiments on ECHO were executed as computational jobs on a cluster composed of 32 Hadoop/HPC 16-core compute nodes and managed with SLURM (Simple Linux Utility for Resource Management). Scripts used for parallelization of simulations are provided as Supplementary Material 'scripts'. The scripts are organized in subfolders with number corresponding to the images in the manuscript, that are the results of the scripts.
As each of the simulations in an in silico experiment is independent from the others, these experiments belong to the so-called group of embarrassingly parallel problems [19]. In this setting, the absence of inter-process communications allows for an arbitrary division of the jobs among computational units (i.e. CPU cores), gaining a performance boost linear with the number of computational units. This means that all available computational units could be exploited optimally: for example, if 512 cores are available, 512 simulations could be run in parallel at the same time. In practice, this significantly increases the feasibility of large in silico experiments: using 512 cores to perform the experiment on which Table X2 and Fig. 2A are based took about 17 h; the same experiment would have taken about one year if all its simulations would have been executed sequentially on a single core.

In silico experiments
In silico experiments described in this paper are based on simulations performed on ECHO. The typical experiment performed on ECHO consisted in making minor changes to the configuration of the network, for example setting the initial activity of some nodes to fixed values in a given range, then simulating the evolution of the resulting system on a period sufficiently long to have the network stabilize at a final state. The activity levels of all nodes in the final state were then collected in tables and analyzed in accordance to the aim of the experiment. In some cases, the final state could not be directly recognized as one of the three possible stable states (SOX9 + , RUNX2 + , Null): that happens because keeping one or more nodes fixed at a certain value may block the evolution of the network in case those nodes need to change before proceeding towards a stable state. In order to resolve those undetermined cases, an additional analysis was performed on the results of the stalled simulations: the nodes artificially held fixed to the given value were let evolve normally (i.e. the artificial block was removed) and an additional simulation was computed starting from the final state of the undetermined simulation. This always led to one of the three possible final states SOX9 + , RUNX2 + , Null, allowing us to classify the outcome of the given modification.

In silico cell differentiation
The in-silico experiments consisted of four main steps: 1. set the initial activity of nodes in ECHO to represent a given cell type, based on our RNASeq data, 2. add the proper factors to represent known differentiation stimuli, based on the protocols that were previously established in our lab (see main manuscript). 3. simulate the model until a steady state is reached, 4. compare the resulting cell type represented by the model with known differentiation results.
The initial configurations we chose represent four cell types: adiposederived stem cells (AMSC), bone marrow-derived stem cells (BMSC), healthy cartilage and embryonic cells. RNAseq data collected in laboratory experiments on these four cell types were used to define the initial configurations of ECHO. In particular, from the RNAseq data we derived the activity levels of nodes in ECHO. As the expression of different genes can have extreme quantitative differences (up to an order of 100,000), the average RNAseq value for each node was translated into activity levels by first normalizing it on a logarithmic scale, and then mapping the resulting values to a [0, 100] linear scale. This was used to represent gene expression values as activity levels in ECHO. According to the in silico experiments previously described in this paper, choosing the activity for the 15 most influential nodes (the "Top 15 nodes", see Fig. 2) in ECHO should be enough to determine the cell fate. We restricted the initialization step of ECHO to only the Top 15 nodes, identified in the experiment for Fig. 2.
For each of the four considered cell types (AMSC, BMSC, cartilage, embryonic), we performed 100,000 random initialization experiments (thus 400,000 experiments in total) with the following initialization settings: -all nodes in ECHO apart from the Top 15 start at a random initial activity and are free to change their value according to the interactions present in ECHO; -the initial activity of the Top 15 nodes is set to the corresponding values derived from the RNAseq data for the chosen cell type.
In each of the 100,000 simulations per cell type, the evolution of the network was simulated in three stages: 1. keep the activity level of the Top 15 nodes fixed at their initial value and perform a simulation starting from the initial configuration for an amount of time, T, sufficient for all interactions to exert their full effect; 2. set a new value for the activity of specific nodes in order to simulate a differentiation stimulus and keep their value fixed, while letting all other nodes (including the Top 15) free to change their value, then continue the simulation for T more time units; 3. let every node in the network free to change its value according to the interactions influencing it and compute the simulation results after T more time units.
Each experiment simulated the network evolution for 3 T time units.

Interleukin addition to the AC model
Using the topology presented in [59], we extended ECHO with the nodes and interactions listed in MethodsX Supplementary table 6. The kvalues representing the interaction speeds were determined using the rules already adopted for the rest of the interactions in ECHO: "slow" reactions such as gene expression have a k-value of 0.1, while "fast" reactions such as phosphorylation have a k-value of 1.0. Note that the interactions iNOS -| DKK1 and iNOS -| GREM occur via gene expression (with iNOS inhibiting mRNA for both targets): for this reason, their kvalue is 0.1.
For each of 1 million simulations, we randomly initialized all nodes in the network and simulated the model until it reached a stable state. After all simulations were done, we computed the percentage of cell fates (see Section 1.4 and Table X2 in MethodsX). Table 1, below, shows the comparison of cell fate percentages between the AC model and the same model after the addition of the IL1β pathway.

Microarray and RNASeq assays
The microarray assays of growth plate and articular cartilage were performed previously [2]. In these microarray experiments, we used paired growth plate (5 samples in total) and articular cartilage (8 samples) from 4 patients. Raw and normalized data were deposited in the Gene Expression Omnibus (GSE-32398).
RNASeq data of articular cartilage, bone marrow Mesenchymal Stem Cells (BMSC), Adipose-derived mesenchymal stem cells (AMSC) or human embryonic stem cells (ESC), were previously performed in the lab of Prof. van Wijnen and the data sets are found here: [48][49][50].

Cell culture and expansion
Human primary chondrocytes (hChs) were obtained from OA regions or relatively healthy looking full thickness cartilage, dissected from knee or hip biopsies of eight patients undergoing total joint replacement (for qPCR: OA: donor # 103 OA knee male-72yo, and donor # 104 OA hip female-73yo; healthy: donor # 105 HL knee male-62, and donor # 44 HL knee female-57yo; for TF-FRAP: OA donor #143 OA female 72yo, #114 female 65yo, HL donor #144 Male of unknown age) according to a previously published method [86]. To isolate cells, the cartilage was digested in chondrocyte proliferation medium containing collagenase type II (0.15% Worthington, NJ, US) for 20-22 h. Subsequently the hChs were expanded at a density of 3000 cells/cm 2 in chondrocyte proliferation medium until the monolayer reached 80% confluency. The chondrocyte proliferation medium consisted of DMEM supplemented with 10% fetal bovine serum (FBS), 1 × non-essential amino acids, 0.2 mM ascorbic acid 2-phosphate (AsAP), 0.4 mM proline, 100 U/ml penicillin and 100 μg/ml streptomycin.
Cells were washed with Imaging buffer, and IL1β (R&D systems, USA or PeproTech, USA) was supplemented to the imaging buffer at a final concentration of 10 ng/ml, unless otherwise noted. FRAP was performed starting from 20 min after treatment, unless otherwise noted. DKK1, FRZB and GREM1 (100 ng/ml each, R&D Systems, Minneapolis, USA) were added 40 min before IL1β addition and the combination was incubated at 37 • C for 1 h. In the treated cells, IL1β (10 ng/ml in single treatment, 5 ng/ml in the co-treatments) and/or DKK1, FRZB and GREM1 (100 ng/ml each) were present in the imaging buffer during the FRAP measurement.
Ten pre-bleach images were taken and the fluorescence intensity values were averaged to normalize the FRAP curve. We measured n ≥ 35 cells per donor, per condition. Matlab™ was used to analyze the FRAP data and the script is available upon request. A diffusion uncoupled, twocomponent method was used to interpret our FRAP data as described previously [61]. Individual FRAP measurement curves were averaged to get a single FRAP curve. To assess the statistical significance between the conditions Mann-Whitney U tests were applied using Origin software (OriginLab ®).
FRAP rates were calculated using: Single component fit : where F 0 is the value of the fluorescent intensity at the first post-bleach frame, A 1 is the amplitude of the fast diffusing population and τ is the time constant.
Two − component fit : where A 2 is the amplitude of slow diffusing population, τ 1 and τ 2 are the time constants of A 1 and A 2 respectively.
Half − time to recover : t ½ = ln(2)*τ Immobile fraction : where F I is the initial intensity and F E is the end value of the recovered intensity.

RNA isolation and qPCR
The hChs were used in passage two unless otherwise stated. When the monolayer reached 60% confluency, IL1β (BioLegend, USA or R&D system, USA) was added at a final concentration of 10 ng/ml for 24 h in the presence of normal expansion medium (DMEM with 10% FBS). DKK1, FRZB and GREM1 (R&D systems) were supplied at a concentration of 200 ng/ml for 24 h, simultaneous with the IL1β). Total RNA was isolated from cells using the NucleoSpin RNA II kit (Macherey-Nagel), according to the manufacturers protocol. The precipitated RNA was dissolved in RNase-free water and subsequently treated with RNase-free DNase I (Invitrogen life technologies, USA). cDNA was obtained from 1 μg of RNA (determined using Nanodrop 2000) with a cDNA synthesis kit Table 1 Cell fate percentage of the steady states in the AC model and in the AC model after addition of IL1β after one million simulations with random initialization of all nodes in the model. (BIO-RAD, Hercules, CA). QPCR was performed using SYBR Green sensimix (Bioline, London, UK) in the Bio-Rad CFX96 (Bio-Rad, Hercules, CA). For each reaction a melting curve was generated to test primer dimer formation and non-specific priming. GAPDH was used for gene expression normalization. All qPCRs were performed in technical triplicates and in 2 patients for each condition. Statistical analysis was performed using student's t-test or one-way ANOVA with Tukey's posthoc. p < .05 was considered statistically significant. Primer sequences are listed below. Primers used for qPCR. Supplementary data to this article can be found online at https://doi. org/10.1016/j.cellsig.2019.109471.

Ethics approval and consent to participate
All fresh cartilage specimens used in this investigation were collected under protocols approved by the Institutional Review Board of the MAYO Clinic or from the Gift of Hope Organ and Tissue Donor Network (Elmhurst, IL) (ORA#: L03090306). The RNASeq data and the gene expression microarray data used for this manuscript was previously published [48][49][50]. Human chondrocytes were isolated from joints that were left over from total joint replacement therapies/arthroplasties performed at the Orthopedisch Centrum Oost Nederland (OCON), Hengelo, The Netherlands. Patients donated these joints with informed consent and after approval of the local hospital medical ethical committees of the ZGT Hengelo and MST Enschede, The Netherlands.

Availability of data and material
All models presented in this manuscript are freely available as a supplement to this manuscript. All models and model updates are freely available at our website https://fmt.cs.utwente.nl/tools/animo. The building of ECHO is described in the accompanying MethodsX manuscript "Building of ECHO: the executable chondrocyte generated with ANIMO".