Novel application of multi-stimuli network inference to synovial fibroblasts of rheumatoid arthritis patients
© Kupfer et al.; licensee BioMed Central Ltd. 2014
Received: 9 July 2013
Accepted: 25 June 2014
Published: 3 July 2014
Network inference of gene expression data is an important challenge in systems biology. Novel algorithms may provide more detailed gene regulatory networks (GRN) for complex, chronic inflammatory diseases such as rheumatoid arthritis (RA), in which activated synovial fibroblasts (SFBs) play a major role. Since the detailed mechanisms underlying this activation are still unclear, simultaneous investigation of multi-stimuli activation of SFBs offers the possibility to elucidate the regulatory effects of multiple mediators and to gain new insights into disease pathogenesis.
A GRN was therefore inferred from RA-SFBs treated with 4 different stimuli (IL-1 β, TNF- α, TGF- β, and PDGF-D). Data from time series microarray experiments (0, 1, 2, 4, 12 h; Affymetrix HG-U133 Plus 2.0) were batch-corrected applying ‘ComBat’, analyzed for differentially expressed genes over time with ‘Limma’, and used for the inference of a robust GRN with NetGenerator V2.0, a heuristic ordinary differential equation-based method with soft integration of prior knowledge.
Using all genes differentially expressed over time in RA-SFBs for any stimulus, and selecting the genes belonging to the most significant gene ontology (GO) term, i.e., ‘cartilage development’, a dynamic, robust, moderately complex multi-stimuli GRN was generated with 24 genes and 57 edges in total, 31 of which were gene-to-gene edges. Prior literature-based knowledge derived from Pathway Studio or manual searches was reflected in the final network by 25/57 confirmed edges (44%). The model contained known network motifs crucial for dynamic cellular behavior, e.g., cross-talk among pathways, positive feed-back loops, and positive feed-forward motifs (including suppression of the transcriptional repressor OSR2 by all 4 stimuli.
A multi-stimuli GRN highly concordant with literature data was successfully generated by network inference from the gene expression of stimulated RA-SFBs. The GRN showed high reliability, since 10 predicted edges were independently validated by literature findings post network inference. The selected GO term ‘cartilage development’ contained a number of differentiation markers, growth factors, and transcription factors with potential relevance for RA. Finally, the model provided new insight into the response of RA-SFBs to multiple stimuli implicated in the pathogenesis of RA, in particular to the ‘novel’ potent growth factor PDGF-D.
KeywordsNetwork modeling Reverse engineering Rheumatoid arthritis Synovial fibroblasts Cytokines Growth factors Cartilage development Multi-stimuli modeling
Networks of molecular components such as genes, proteins, and metabolites play a crucial role in systems biology. Since high-throughput gene expression data are now more easily affordable, inference and modeling of gene regulatory networks (GRNs) has become more important over the last years. The regulation of gene expression can be visualized by networks and used to obtain new insights into the underlying biological mechanisms.
GRN inference from gene expression data is a widely used and accepted approach to reconstruct networks in systems biology. In this context, GRNs summarize gene regulatory interactions, including the regulation of gene expression by extra- and intracellular stimuli. The inferred networks open the possibility to better understand the underlying processes and cellular responses to manipulation and represent a starting point for the extraction of biological hypotheses.
Commonly, a GRN describing a biological network is a graph G = (V,E) where V represents the components (nodes) and E the relationships (edges) between the components. In the case of a GRNs, nodes represent genes and edges stand for transcriptional regulation [1, 2].
There are several inference methods, each using different sources and modeling assumptions that may lead to different results and visualizations. To address GRN inference from time series data, several methods and approaches have been used. For example there are vector autoregressive models [3–6], linear Bayesian networks [7, 8] and ordinary differential equation (ODE)-based approaches [9–11]. Regarding the fact that multi-stimuli experiments often lead to complex networks, especially if the data are time-resolved, heuristic network inference approaches are appropriate to handle the high number of possible structural connection parameters. Heuristic approaches possess the ability to reduce the computation time for network construction and still provide satisfactory inference results.
To our knowledge, there are only few heuristic methods for the inference of multi-stimuli experiments [12–14]. This type of experiments aims at investigating the relative importance of different stimuli for physiological and pathological processes, which may depend on more than one stimulus. In this case, the term multi-stimuli experiments is commonly used in the literature [12–15].
To address the challenge of GRN inference from multi-stimuli, time-resolved gene expression data, the heuristic inference algorithm NetGenerator V2.0 was chosen in the present study . The main reason to select this method is its ability to integrate prior knowledge obtained from different sources. This leads to a network that combines both expression data and prior knowledge, thus showing the capability of generating meaningful results in various biological and medical fields.
In the present study, the transcriptional regulation in synovial fibroblasts (SFBs) isolated from rheumatoid arthritis (RA) patients was studied by modeling the response to 4 external stimuli (IL-1 β, TNF- α, TGF- β, and PDGF-D). RA is a chronic inflammatory and destructive joint disease perpetuated by an invasive synovial membrane, the so-called pannus tissue. Activated and semi-transformed SFBs are key components of the inflamed synovial membrane [16, 17]. In the normal joint, SFBs are characterized by a balanced expression of proteins regulating the formation and degradation of the extracellular matrix (ECM). In RA, however, SFBs are permanently activated by cytokines, e.g., TNF- α and IL-1 β, which are potent pro-inflammatory cytokines especially produced by macrophages [16, 18]. Activated SFBs excessively express and secrete tissue-degrading enzymes such as matrix-metalloproteases (MMPs) or soft matrix components (e.g., collagens), thus both maintaining the degradation of ECM, cartilage, and bone, and inducing fibrosis of the affected joints [18, 19]. Moreover, SFB contribute to joint inflammation by increased expression of additional growth factors such as TGF- β or PDGF-D [17, 20, 21]. As a consequence, autocrine mechanisms are assumed to play a key role in synovial hyperplasia and the enduring activation of SFB . For instance, TGF- β enhances its own expression  and that of PDGF family proteins [23, 24]. TGF- β, although also exhibiting deactivating features , is known to support matrix remodeling/fibrosis and initial activation of inflammatory processes. In contrast, PDGF proteins act as potent growth factors for several cell types in the synovial membrane, including SFBs , and may serve as inducers of MMPs . In addition, TGF- β and PDGF-D are able to amplify the effects of other cytokines. When combined, both cytokines augment the secretion of pro-inflammatory and pro-destructive proteins by SFB ; also, TGF- β and PDGF-D have been independently shown to enhance the effects of IL-1 β[23, 25].
Although several reports have covered the effects of single or a combination of selected cytokines on SFB and their influence on RA pathogenesis (reviewed in ), SFB gene expression in response to 4 disease-relevant stimuli (IL-1 β, TNF- α, TGF- β, and PDGF-D) and the subsequent inference of GRNs have not been addressed to date. Therefore, this study provides new hypotheses for the interdependent regulation of SFB-derived gene expression profiles under the influence of different cytokines and growth factors.
Clinical data of patients
Disease duration (years)
NSAIDs, leflunomide, prednisone
NSAIDy, prednisone, MTX
Celecoxib, prednisone, MTX, leflunomide
Prednisone, MTX, alendronate
NSAIDs, prednisone, azathioprin
NSAIDs, MTX, etanercept (TNF-blocker)
Stimulation of SFBs
Cell stimulation and isolation of total RNA
At the end of the fourth passage, the SFB were washed in serum-free Dulbeccos modified Eagle’s medium (DMEM) and then stimulated with 10 ng/ml of either IL-1 β, TNF- α, TGF- β, or PDGF-D (PeproTech, Hamburg, Germany) in serum-free DMEM for 0, 1, 2, 4, or 12 hours (see Table 2). At each time point, the medium was removed and the cells were harvested after treatment with trypsin (0.25% in Versene; Invitrogen, Karlsruhe, Germany). After washing with phosphate-buffered saline, cells were lysed with RLT buffer (Qiagen, Hilden, Germany) and frozen at -70°C. Total RNA was isolated using the RNeasy Kit (Qiagen) according to the supplier’s recommendations.
The analysis of gene expression was performed using HG-U133 Plus 2.0 RNA microarrays (Affymetrix, Santa Clara, CA, USA; total of 120 microarrays). Labeling of RNA probes, hybridization, and washing were carried out according to the supplier’s instructions. Microarrays were analyzed by laser scanning (Gene Scanner, Hewlett-Packard, Palo Alto, CA, USA). The data for the stimuli TNF- α and TGF- β are accessible through Gene Expression Omnibus series accession number GSE13837; the data for the stimuli IL1- β and PDGF-D through Gene Expression Omnibus series accession number GSE58203. Since several studies have demonstrated that alternative Chip Definition Files (CDF) for gene annotation resolve the problem of choosing reliable and non-contradictory probe sets for each transcript, the CDF presented by Ferrari et al. was used in the present study . This allows to reduce the effects of cross-hybridization and other system-based biases [30–33]. Robust Multichip Average (RMA) was used with the default configuration for background adjustment and normalization .
For batch correction of microarray data, the non-parametric prior method of ComBat was used, which represents an Empirical Bayesian (EB) method . EB methods are well-suited for the analysis of microarrays since they are able to handle high-dimensional data from small sample sizes in a robust manner. For this purpose, EB methods borrow information from certain genes in order to obtain improved estimates for the expression of all genes. This capability of shrinking the variance across all genes has been shown in several publications [36–43]. Based on this advantage, Johnson et al. extended the EB methods with a location and scale (L/S) adjustment, which adjusts batches with small sample sizes to each other and avoids special normalization procedures [35, 44]. Based on the available R-package of Johnson, a modified method was developed for using RMA-normalized instead of dChip-normalized data . The Sample Information File was created, which, besides the creation date, contains the microarray name, time point (0, 1, 2, 4, and 12 h), and treatment (IL-1 β, TNF- α, TGF- β, and PDGF-D) of each sample, all serving as covariates for the ComBat method. The creation date was tagged as ‘batch effect’, necessary for the correction of possible batch effects by ComBat.
For the identification of differentially expressed genes (DEGs) in microarray experiments, the R-package LIMMA was used . LIMMA is commonly used in the analysis of microarray data, designed to analyze complex experiments involving simultaneous comparisons between many RNA targets. Regarding the identification of DEGs, the cardinal concept is to fit a linear model to the expression data for each gene. The analysis is applied to microarray expression data, which are represented in a matrix consisting of probe sets (genes; rows) and arrays (columns). LIMMA requires a design matrix representing different RNA targets, as well as a contrast matrix assigning the coefficients of the design matrix to the contrasts of interest (i.e., expression over time, disease status, and/or treatment). In our case, the contrasts of interests reflect the differential gene expression over time for each stimulus (IL-1 β, TNF- α, TGF- β, and PDGF-D). The lmFit function fits the gene-wise linear model to the microarray data. DEGs were obtained using the implemented top-table function and user-specific thresholds (i.e., 2-fold change; p-value ≤ 0.05) as recommended in . The Benjamini and Hochberg method was used to correct the obtained p-values for multiple testing .
To test the association between GO categories and the identified DEGs, the R-package GOstats was used. The implemented conditional hypergeometric (hg) test uses the relationships among GO categories to address the hierarchical structure of the GO database . For testing the associations between GO and the list of selected genes, a universe has to be defined containing all genes on the microarray. In addition, the cutoff for the adjusted p-value of the hg test has to be set to 0.05. The summary containing the enriched GO categories and their significance level is represented in the GOHyperGResult. The result of GOstats also contains the actual gene count for each of the significant GO categories.
Pathway Studio 9 (PS9), formerly known as PathwayAssist, was used to extract prior knowledge for the DEGs . PS9 is an analysis software for pathways, gene regulation networks, and protein interaction maps. It allows for an interpretation of microarray and proteomics data, classification of proteins, drawing of pathway diagrams, as well as export, import, and filtering of data. PS9 includes the proprietary ResNet Mammalian Database 9 built from 20,000,000 abstracts in PubMed, as well as over 880,000 full-text articles . All literature information regarding the analyzed genes was manually validated in the respective publications by 2 long-term experts in the field of (experimental) rheumatology (R.H.; R.W.K.), as described previously . This curation focused on the appearance and the temporal behavior of the following features: (i) constitutive vs. induced gene expression; (ii) co-expression vs. divergent expression of mediators, TFs, and target genes; (iii) expression of mediators/transcription factors vs. expression of target genes; (iv) regulation of target gene expression based on the expression of different transcription factors; (v) expression of individual genes vs. expression of their functional groups; and (vi) discrepancies to the literature. Subsequently, the extracted interactions were assessed with respect to biological coherence and relevance.
The sum of weighted gene expression of N genes and weighted input u(t) describes the dynamic change of expression x i of gene i. The terms u k (t) represent the external stimulation (IL-1 β, TNF- α, TGF- β, and PDGF-D; M = 4) in a step function (u k (t < 0) = 0 and u k (t > = 0) = 1). The interaction parameters a i, j and the input parameters b i,k model the regulatory interactions. The number of potential interactions sums up to N 2 + M · N, where N is the number of genes finally selected for modeling and M is the number of stimuli.
To denote potential (positive and negative) interactions in the network, the parameters a i, j and b i,k were estimated, with positive values representing activating connections and negative values representing repression. The model interaction parameters, which have to be determined by the NetGenerator V2.0 algorithm, characterize the structure of the GRN. The main component of the heuristic algorithm is an optimization procedure which minimizes the number of non-zero parameters (represented as edges in the network) required to achieve a good fit of the simulated model kinetics to the measured time series data. The NetGenerator V2.0 applies explicit structure optimization involving iterative construction of a sparse sub-model [12, 51–53].
One of the advantages of the NetGenerator V2.0 algorithm is the possibility to integrate prior knowledge. This knowledge is received from diverse resources such as the published literature. Since the extracted knowledge is independent of the time series data, it provides additional information for the inference process. To represent the prior knowledge, a separate interaction matrix is established assigning specific values to the interactions among genes. This information is encoded by the values shown in Additional file 1. To obtain the best results for the fitting of the inferred network model to the expression data, it is necessary that the provided knowledge is flexibly integrated (referred to as soft integration in the literature ). This avoids under-fitting or over-fitting as often observed upon hard knowledge integration . The best solution is a balance between low network complexity (i.e., the lowest possible number of edges) and the lowest possible deviation between simulated and measured expression values (i.e., average mean squared error - MSE). The most important parameter is the ‘allowedError’, which controls this balance. In this context, the ‘allowedError’ represents the maximally allowed error for any gene, stimulus, and/or time point. To determine the optimal value for the ‘allowedError’, the resulting models of the inference runs were analyzed. The optimized ‘allowedError’ resulting in a low deviation (MSE), a low number of network edges, and a high number of prior knowledge edges, was chosen for further analysis. The selected model was validated regarding the robustness of the inferred network against small changes in the time series data, reflecting measurement error due to technical or biological variance. Therefore, Gaussian distributed () random noise was added to the original data. This procedure was repeated 100 times, leading to a series of inferred models. These models were analyzed concerning the occurrence of the edges of the undisturbed network model. Edges with an absolute frequency > 51 were regarded as stable and consequently integrated into the final consensus model.
Analysis of differentially expressed genes
Number of DEGs and their union
GO analysis and the top 10 GO terms resulting from GO analysis
Response to cytokine stimulus
Cellular response to cytokine stimulus
Cytokine-mediated signaling pathway
Regulation of cell proliferation
Regulation of cartilage development
Cellular biosynthetic process
Cellular macromolecule biosynthetic process
Cellular response to chemical stimulus
Positive regulation of cellular process
Genes contained in the most significant GO term ‘cartilage development’
Odd-skipped related 2 (Drosophila)
Distal-less homeobox 2
Fibroblast growth factor 2 (basic)
Growth hormone receptor
GLI family zinc finger 2
GLI family zinc finger 3
Low density lipoprotein receptor-related protein 6
SMAD family member 3
msh homeobox 1
Carbohydrate (chondroitin 4) sulfotransferase 11
Paired related homeobox 2
Parathyroid hormone-like hormone
v-rel reticuloendotheliosis viral oncogene homolog A (avian)
Bone morphogenetic protein 4
Snail homolog 2 (Drosophila)
Snail homolog 1 (Drosophila)
SRY (sex determining region Y)-box 9
Transforming growth factor, beta 1
Wingless-type MMTV integration site family, member 9A
High mobility group AT-hook 2
Wingless-type MMTV integration site family, member 5B
Fibroblast growth factor 18
We have used GO analysis for an approach to the potential functional relevance of the 1914 differentially expressed genes, since it is one of the best available sources of information for TF-gene or gene-gene interactions and their biological importance. The choice of the 24 DEGs from the GO term was further supported by the fact that 11/24 genes also appeared in the independent GO term ‘regulation of cartilage development’. In addition, the 24 chosen genes were manually checked by the above 2 long-term experts in the field of (experimental) rheumatology (R.H.; R.W.K.) for consistency with the known literature as published previously . Also, the 24 chosen genes contain 13 transcription factors, 8 secreted factors, and 3 genes with other functions (now color-highlighted in Table 5), indicating a good balance between regulating factors and effector/target molecules likely to represent the (patho)physiological processes.
Prior knowledge for the 24 DEGs of the GO term ‘cartilage development’ was extracted using PS9 (see Additional file 3). The gene-to-gene interactions were subsequently formatted for input into the network inference tool NetGenerator V2.0 (see Additional file 1). All extracted knowledge was manually verified in order to avoid doubtful and misleading statements in the respective publications.
Subsequently, a stability analysis (also called ‘internal validation’) was performed to investigate the robustness of the model. The main reason for this investigation is to avoid an over-fitting of the inferred model to the measured data. Since minor data variability (i.e., noise) should not change the structure of the GRN model, the disturbed models should show a high structural similarity to the initial network. Therefore, noise was added to the expression values and further network inference was performed.
Three sources of additional literature knowledge were used for ‘external’ validation of the consensus model: (i) stimuli-to-gene interactions extracted by PS9 (2013/02/18) not considered for network model inference (see Additional file 5); (ii) gene-to-gene interactions integrated into PS9 in the period between network inference (2012/09/12) and external model validation (2013/02/18); (iii) interactions resulting from manual literature searches. This ‘external’ model validation resulted in the confirmation of 10 more edges (1 gene-to-gene and 9 stimuli-to-gene interactions; indicated by the respective 10 green edges in Figure 3). The 7/10 interactions derived from PS9 analysis included 6 interactions as named above in (i) and one interaction as named in (ii) (the gene-to-gene edge EDN1 to SNAI1 ).
Regarding the whole process of ‘internal’ and ‘external’ validation, the number of edges confirmed by literature knowledge increased from 17 to 25 from the initial to the consensus network (compare Figure 3 to Additional file 4).
To our knowledge, this is the first network modeling the complex, time-resolved concerted action of 4 different disease-relevant mediators (TNF- α, IL-1 β, TGF- β, and PDGF-D) on the gene expression in RA-SFBs.
On one hand, several known effects of each stimulus (e.g., induction of transcription factors such as SMAD3, SNAI2, and GLI2 by TGF- β[62–64] or RELA by TNF- α and IL-1 β) were represented in the present network. Such transcription factors regulate the expression of their target genes, thus controlling the subsequent cellular response to the (combination of) different stimuli. This includes the steering of important pathophysiological processes in RA, such as ECM formation and fibrosis induction (guided by SMAD3; ), proliferation, cell survival, and inflammation (driven by RelA; ), as well as differentiation or dedifferentiation (controlled by SNAI2  or GLI2 ). On the other hand, the present network predicted unexpected regulatory connections (e.g., induction of GLI2 by PDGF-D).
Regarding the network inference, NetGenerator V2.0 was successfully applied to model a medium-scale (24 genes and 4 stimuli), robust network of moderate complexity (57 edges). The full interaction matrix of the linear model for N genes used in the present study has N ∗ N elements (a i,j in the eq. in the chapter ‘Network Inference’ above) and thus requires at least N samples with N gene expression values each. In addition, for M stimuli (i.e., 4 in the present study) a number of M∗N elements (b j, k ) has to be added. Thus, in the present study with 24 genes and 4 stimuli, a total of at least 28 samples would be required to identify 672 edges. As a consequence, the number of inferred edges (57/672) is only 8.5% of that in the fully connected network. This selection of the most reliable edges is driven by both the criterion of sparseness and prior knowledge. There was an impressive fit of the simulated expression values to the measured values with a small average mean squared error, in particular with regard to the fact that the response of RA-SFB to 4 different stimuli was simultaneously modeled. In addition, a high percentage of literature-confirmed edges (15/57) could be used for inference of the consensus network. Strikingly, post-inference, ‘external’ network validation resulted in the confirmation of 10 additional network edges, thus yielding a total of 25/57 (44%) literature-supported network edges. The remaining 32 predicted edges (consisting of only one contradictory edge and 31 edges without literature information) still require verification by future wet-lab experiments.
However, the algorithm used for network inference in the present study also shows some limitations. A detailed mechanistic model for transcriptional GRN would also have to consider important steps such as processing, transport, translation, and degradation of mRNA, or else the parallel existence of numerous interacting molecules such as transcription factors, (phosphorylated) proteins, micro-RNA etc. In addition, effects of post-translational protein modifications, interactions with co-factors, and intracellular localization should be considered. However, if medium scale models - such as in the present study - are to be analyzed, mechanistic modeling reflecting all of the above processes is presently impossible due to the lack of detailed data and/or knowledge. Our current aim was to use the gene expression profiles of stimulated RA-SFBs to uncover relevant and unknown stimuli-to-gene or gene-to-gene relations, and thus to identify potential key regulators for RA.
In contrast to biological systems consisting of numerous individual genes and their multiple reciprocal interactions, network modeling based on limited gene expression data sets (120 microarrays derived from 5 time points, 6 biological replicates, and 4 stimuli) requires restricted complexity (i.e., a low number of genes; 24 in the present study) and a limited number of interactions (i.e., stimuli-to-gene or gene-to-gene interactions; 57 in the present study). The aim of restricted network complexity was also achieved by only assuming linear associations between the individual components of the network, allowing a more reliable modeling of data with considerable noise. The inference with complex, unknown non-linear elements would require a larger set of noise-free data (which is not realistic for in-vivo and in-vitro studies, and likely difficult to finance) to test for the best type of non-linearity and to identify its additional parameters (see section 3.3.2 in ). Although known non-linear relationships could be easily integrated into the network inference algorithm, this extension would still require prior knowledge about the type of relationship and a substantially larger data set (i.e., several hundred microarrays). By using cells from 6 different biological replicates (i.e., RA patients; see Table 2) in the present data set, a special attempt was made to consider the biological variance in diseased individuals.
The main novelty of the present approach for dynamic GRN inference is the possibility to consider the effects of multiple stimuli, which permits the advanced, simultaneous analysis of multi-experiment time series expression data, in our case together with the option of soft integration of prior knowledge. For example, NetGenerator was capable of modeling the response following the initial stimulation of RA-SFBs by TNF- α, IL-1 β, TGF- β, and PDGF-D, e.g., sequential activation of signaling molecules/transcription factors (see 13/24 potential transcription factors in Table 5) and/or protein secretion of (growth) factors (8/24 potential factors in Table 5), which may support enduring activation of RA-SFBs  and persistence of joint inflammation and destruction in RA . In addition, NetGenerator suggested that a variety of factors is induced in response to more than one stimulus. For instance, TGF- β, WNT9A (also known as WNT14), and PTHLH (also known as PTHRP; mediated via SMAD3) are induced by both TGF- β and PDGF-D, which is in good agreement with the literature [56, 58, 59].
The multi-stimuli approach also creates the basis to predict possible cross-talk between signaling pathways distal to the mediator-receptor level with potential relevance for the pathogenesis of RA. For instance, the network suggests parallel possibilities for the induction of the gene expression of RELA by IL-1 β and TFN- α, either directly or mediated by suppression of the transcriptional repressor OSR2 (an example for a positive feed-forward network motif). In addition, our model suggests an induction of FGF2 and FGF18 by all 4 stimuli via the transcription factors GLI2, PRRX2, MSX1, and/or OSR2. In this context, the predicted suppression of the transcriptional repressor OSR2 by all 4 stimuli analyzed in this study (as well as the influence of 2 or 3 stimuli on the genes TGF- β, RELA, GLI2, WNT5B, and WNT9A) may well indicate the existence and/or importance of respective cross-talk among these mediators in RA-SFBs. Thus, our predictive model provides new information about the sequential cascade regulation of gene expression, since the underlying indirect activation pathways are still inadequately characterized in the literature.
However, the network also revealed opposing regulatory effects, e.g., suppression of WNT5B by TNF- α and IL-1 β, but induction by TGF- β, thus reflecting the differential characteristics of the respective mediators. This may reflect the complex, mostly indirect induction of gene expression via subsequent signaling molecules by external pro-inflammatory cytokines/ growth factors such as IL-1 β, TNF- α, and TGF- β.
Also gene-regulatory (sub-) cycles were identified within the network. Two positive feedback loops were predicted by the inferred network (i.e., the loop between BMP4, GLI3, and PTHLH, and the loop between TGFB1, EDN1, and FGF2). These positive feedback loops bear pathogenetic potential for RA since deregulated expression of these relevant genes could contribute to pro-inflammatory or pro-destructive processes in RA [69–76]. For instance, TGF- β/SMAD3-induced PTHLH could suppress BMP4 expression, thereby suppressing BMP4-driven transcription of the negative transcription factor GLI3, and thereby enhancing its own transcription. This circuit could provide a functional basis for the enhanced expression of PTHLH and reduced amounts of BMP4 in RA synovial membrane and fluid, in analogy to previous reports [72, 77]. However, since no explicit evidence for the existence of such feedback loops was found in the respective literature, this issue remains a target of future studies.
Interestingly, the inferred GRN showed a relatively high number of stimuli-to-gene edges in the case of TGF- β and PDGF-D (9 each) and a lower number of such edges for TNF- α and IL-1 β (5 and 3, respectively). In addition, 6/9 of the TGF- β-to-gene edges were confirmed by literature, whereas only 1 PDGF-D-to-gene edge was validated by literature, underlining the novelty of pathogenetic PDGF-D effects in RA (; present model). The relative preference of stimuli-to-gene edges for TGF- β and PDGF-D is likely caused by the overrepresentation of DEGs for PDGF-D (see Table 3) and the exclusive use of the most significant GO term ‘cartilage development’ for network inference (see Table 4). Of the 8 non-confirmed, novel PDGF-D-to-gene edges, particularly WNT9A (stability score of 92/100; Figure 3) and possibly MSX1 (score of 86/100) would be the most attractive ‘key regulator’ targets for experimental validation.
Among the genes emphasized in the above discussion for the structural features of the inferred model, WNT9A, PTHLH, FGF2, and FGF18 are secreted proteins involved in tissue development, especially formation of cartilage and bone [69, 72, 75]. In addition, PTHLH has been shown to mediate anti-proliferative effects and to induce the production of matrix-degrading enzymes , which may participate in cartilage and bone destruction . Finally, FGF2 supports angiogenesis and inflammation in the synovial membrane . Therefore, these factors may trigger the activation of RA-SFBs and other articular cell types such as chondrocytes or osteoblasts [20, 71, 74, 75, 78].
NetGenerator V2.0 was successfully applied to model a dynamic, medium-scale (24 genes and 4 stimuli), robust GRN of moderate complexity (57 edges). In addition, the predicted GRN showed a high reliability, since 10 predicted edges were independently validated by literature findings after completion of the inference process. Also, the model reflects known network motifs that are crucial for dynamic cellular behavior, such as cross-talk among pathways, positive feed-forward motifs, and positive feed-back loops. Finally, the model provides new insight into the response of pathogenetically relevant RA-SFBs to multiple stimuli implicated in the pathogenesis of RA, especially the ‘novel’ potent growth factor PDGF-D. In particular, transcription factors such as SMAD3, SNAI2, and GLI2 (induced by TGF- β), or RELA (induced by TNF- α and IL-1 β), as well as the secreted factors WNT5B (suppressed by TNF- α and IL-1 β, but induced by TGF- β), and WNT9A (induced by both TGF- β and PDGF-D) are examples of complex network interactions resulting from the present study and may indicate very attractive ‘key regulator’ targets for experimental validation and/or therapy.
This work was supported by grants from the German Federal Ministry of Education and Research (BMBF FKZ 0315719A and FKZ 0315719B, ERASysBio PLUS, LINCONET; FKZ 0315736, Virtual Liver Network; FKZ 01EC1009A; ArthroMark) as well as the EU (IMI grant agreement no. 115142; BeTheCure).
- Babu MM, Luscombe NM, Aravind L, Gerstein M, Teichmann SA: Structure and evolution of transcriptional regulatory networks. Curr Opin Struct Biol. 2004, 14 (3): 283-291. [PMID: 15193307].View ArticlePubMedGoogle Scholar
- Davidson EH, Rast JP, Oliveri P, Ransick A, Calestani C, Yuh CH, Minokawa T, Amore G, Hinman V, Arenas-Mena C, Otim O, Brown CT, Livi CB, Lee PY, Revilla R, Rust AG, Schilstra MJ, Clarke PJC, Arnone MI, Rowen L, Cameron RA, McClay DR, Hood L, Bolouri H, Pan Zj: A genomic regulatory network for development. Science. 2002, 295 (5560): 1669-1678. [PMID: 11872831].View ArticlePubMedGoogle Scholar
- Bolstad A, Van Veen BD, Nowak R: Causal network inference via group sparse regularization. IEEE Trans Signal Process. 2011, 59 (6): 2628-2641. [PMID: 21918591].View ArticlePubMedPubMed CentralGoogle Scholar
- Meinshausen N, Bühlmann P: High-dimensional graphs and variable selection with the Lasso. Ann Stat. 2006, 34 (3): 1436-1462.View ArticleGoogle Scholar
- Morrissey ER, Juárez MA, Denby KJ, Burroughs NJ: On reverse engineering of gene interaction networks using time course data with repeated measurements. Bioinformatics. 2010, 26 (18): 2305-2312. [PMID: 20639410].View ArticlePubMedGoogle Scholar
- Opgen-Rhein R, Strimmer K: Learning causal networks from systems biology time course data: an effective model selection procedure for the vector autoregressive process. BMC, Bioinformatics. 2007, 8 Suppl 2: S3-[PMID: 17493252].View ArticlePubMedGoogle Scholar
- Kim SY, Imoto S, Miyano S: Inferring gene networks from time series microarray data using dynamic Bayesian networks. Brief Bioinform. 2003, 4 (3): 228-235. [PMID: 14582517].View ArticlePubMedGoogle Scholar
- Rau A, Jaffrézic F, Foulley JL, Doerge RW: An empirical Bayesian method for estimating biological networks from temporal microarray data. Stat Appl Genet Mol Biol. 2010, 9: Article 9-[PMID: 20196759].PubMedGoogle Scholar
- Nam D, Yoon SH, Kim JF: Ensemble learning of genetic networks from time-series expression data. Bioinformatics. 2007, 23 (23): 3225-3231. [PMID: 17977884].View ArticlePubMedGoogle Scholar
- Bansal M, di Bernardo D: Inference of gene networks from temporal gene expression profiles. IET Syst Biol. 2007, 1 (5): 306-312. [PMID: 17907680].View ArticlePubMedGoogle Scholar
- Li CW, Chen BS: Identifying functional mechanisms of gene and protein regulatory networks in response to a broader range of environmental stresses. Comp Funct Genom. 2010, 2010: 20-doi:10.1155/2010/408705.View ArticleGoogle Scholar
- Weber M, Henkel SG, Vlaic S, Guthke R, van Zoelen EJ, Driesch D: Inference of dynamical gene-regulatory networks based on time-resolved multi-stimuli multi-experiment data applying NetGenerator V2.0. BMC, Syst Biol. 2013, 7: 1-[PMID: 23280066].View ArticleGoogle Scholar
- Wang Y, Joshi T, Zhang XS, Xu D, Chen L: Inferring gene regulatory networks from multiple microarray datasets. Bioinformatics. 2006, 22 (19): 2413-2420. [PMID: 16864593].View ArticlePubMedGoogle Scholar
- Gupta R, Stincone A, Antczak P, Durant S, Bicknell R, Bikfalvi A, Falciani F: A computational framework for gene regulatory network inference that combines multiple methods and datasets. BMC, Syst Biol. 2011, 5: 52-[PMID: 21489290].View ArticleGoogle Scholar
- Kupfer P, Sebastian V, Huber R, Kinne RW, Guthke R: Different stimuli for inference of gene regulatory network in rheumatoid arthritis. Proceedings of the International Conference on Bioinformatics Models, Methods and Algorithms, Volume 1. Edited by: Gamboa H., Fred ALN, Solé-Casals J, Fernandes P. 2013, Barcelona: Scitepress - Science and Technology Publications, 282-287. [ISBN: 978-989-8565-35-8].Google Scholar
- Smolen JS, Steiner G: Therapeutic strategies for rheumatoid arthritis. Nat Rev Drug Discov. 2003, 2 (6): 473-488. [PMID: 12776222].View ArticlePubMedGoogle Scholar
- Kinne RW, Palombo-Kinne E, Emmrich F: Activation of synovial fibroblasts in rheumatoid arthritis. Ann Rheum Dis. 1995, 54 (6): 501-504. [PMID: 7632096].View ArticlePubMedPubMed CentralGoogle Scholar
- Firestein GS: Evolving concepts of rheumatoid arthritis. Nature. 2003, 423 (6937): 356-361. [PMID: 12748655].View ArticlePubMedGoogle Scholar
- Choy E: Understanding the dynamics: pathways involved in the pathogenesis of rheumatoid arthritis. Rheumatology. 2012, 51 Suppl 5: v3-11. [PMID: 22718924].View ArticlePubMedGoogle Scholar
- Karouzakis E, Neidhart M, Gay RE, Gay S: Molecular and cellular basis of rheumatoid joint destruction. Immunology Lett. 2006, 106: 8-13. [PMID: 16824621].View ArticleGoogle Scholar
- Pohlers D, Huber R, Ukena B, Kinne RW: Expression of platelet-derived growth factors C and D in the synovial membrane of patients with rheumatoid arthritis and osteoarthritis. Arthritis Rheum. 2006, 54 (3): 788-794. [PMID: 16508943].View ArticlePubMedGoogle Scholar
- Schett G, McInnes IB: Cytokines in the pathogenesis of rheumatoid arthritis. Nat Rev Immunol. 2007, 7 (6): 429-442. [PMID: 17525752].View ArticlePubMedGoogle Scholar
- Wilder RL, Lafyatis R, Roberts AB, Case JP, Kumkumian GK, Sano H, Sporn MB, Remmers EF: Transforming growth factor-beta in rheumatoid arthritis. Ann N Y Acad Sci. 1990, 593: 197-207. [PMID: 2165375].View ArticlePubMedGoogle Scholar
- Sakuma M, Hatsushika K, Koyama K, Katoh R, Ando T, Watanabe Y, Wako M, Kanzaki M, Takano S, Sugiyama H, Hamada Y, Ogawa H, Okumura K, Nakao A: TGF-beta type I receptor kinase inhibitor down-regulates rheumatoid synoviocytes and prevents the arthritis induced by type II collagen antibody. Int Immunol. 2007, 19 (2): 117-126. [PMID: 17135447].View ArticlePubMedGoogle Scholar
- Niedermeier M, Pap T, Korb A: Therapeutic opportunities in fibroblasts in inflammatory arthritis. Best Pract Res Clin Rheumatol. 2010, 24 (4): 527-540. [PMID: 20732650].View ArticlePubMedGoogle Scholar
- Rosengren S, Corr M, Boyle DL: Platelet-derived growth factor and transforming growth factor beta synergistically potentiate inflammatory mediator synthesis by fibroblast-like synoviocytes. Arthritis Res Ther. 2010, 12 (2): R65-[PMID: 20380722].View ArticlePubMedPubMed CentralGoogle Scholar
- Arnett FC, Edworthy SM, Bloch DA, Fries JF, Cooper NS, Healey LA, Kaplan SR, Liang MH, Luthra HS, McShane DJ: The American Rheumatism Association 1987 revised criteria for the classification of rheumatoid arthritis. Arthritis Rheum. 1988, 31 (3): 315-324. [PMID: 3358796].View ArticlePubMedGoogle Scholar
- Zimmermann T, Kunisch E, Pfeiffer R, Hirth A, Stahl HD, Sack U, Laube A, Liesaus E, Roth A, Palombo-Kinne E, Emmrich F, Kinne RW: Isolation and characterization of rheumatoid arthritis synovial fibroblasts from primary culture–primary culture cells markedly differ from fourth-passage cells. Arthritis Res. 2001, 3: 72-76. [PMID: 11178129].View ArticlePubMedGoogle Scholar
- Ferrari F, Bortoluzzi S, Coppe A, Sirota A, Safran M, Shmoish M, Ferrari S, Lancet D, Danieli GA, Bicciato S: Novel definition files for human GeneChips based on GeneAnnot. BMC, Bioinformatics. 2007, 8: 446-[PMID: 18005434].View ArticlePubMedGoogle Scholar
- Lu J, Lee JC, Salit ML, Cam MC: Transcript-based redefinition of grouped oligonucleotide probe sets using AceView: high-resolution annotation for microarrays. BMC Bioinformatics. 2007, 8: 108-[PMID: 17394657].View ArticlePubMedPubMed CentralGoogle Scholar
- Chalifa-Caspi V, Yanai I, Ophir R, Rosen N, Shmoish M, Benjamin-Rodrig H, Shklar M, Stein TI, Shmueli O, Safran M, Lancet D: GeneAnnot: comprehensive two-way linking between oligonucleotide array probesets and GeneCards genes. Bioinformatics. 2004, 20 (9): 1457-1458. [PMID: 17394657].View ArticlePubMedGoogle Scholar
- Draghici S, Khatri P, Eklund AC, Szallasi Z: Reliability and reproducibility issues in DNA microarray measurements. Trends Genet. 2006, 22 (2): 101-9. [PMID: 16380191].View ArticlePubMedGoogle Scholar
- Harbig J, Sprinkle R, Enkemann SA: A sequence-based identification of the genes detected by probesets on the Affymetrix U133 plus 2.0 array. Nucleic Acids Res. 2005, 33 (3): e31-[PMID: 16380191].View ArticlePubMedPubMed CentralGoogle Scholar
- Bolstad BM, Irizarry RA, Astrand M, Speed TP: A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics. 2003, 19 (2): 185-193. [PMID: 12538238].View ArticlePubMedGoogle Scholar
- Johnson WE, Li C, Rabinovic A: Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics. 2007, 8: 118-127. [PMID: 16632515].View ArticlePubMedGoogle Scholar
- Kendziorski CM, Newton MA, Lan H, Gould MN: On parametric empirical Bayes methods for comparing multiple groups using replicated gene expression profiles. Stat Med. 2003, 22 (24): 3899-3914. [PMID: 14673946].View ArticlePubMedGoogle Scholar
- Chen Y, Dougherty ER, Bittner ML: Ratio-based decisions and the quantitative analysis of cDNA microarray images. J Biomed Opt. 1997, 2 (4): 364-374. [PMID: 23014960].View ArticlePubMedGoogle Scholar
- Efron B, Tibshirani R, Storey JD, Tusher V: Empirical Bayes analysis of a microarray experiment. J Am Stat Assoc. 2001, 96 (456): 1151-1160. [http://dx.doi.org/10.2307/3085878].View ArticleGoogle Scholar
- Newton MA, Kendziorski CM, Richmond CS, Blattner FR, Tsui KW: On differential variability of expression ratios: improving statistical inference about gene expression changes from microarray data. J Comput Biol. 2001, 8: 37-52. [PMID: 11339905].View ArticlePubMedGoogle Scholar
- Tusher VG, Tibshirani R, Chu G: Significance analysis of microarrays applied to the ionizing radiation response. Proc Natl Acad Sci U S A. 2001, 98 (9): 5116-5121. [PMID: 11309499].View ArticlePubMedPubMed CentralGoogle Scholar
- Smyth GK: Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol. 2004, 3: Article3-[PMID: 16646809].PubMedGoogle Scholar
- Lönnstedt I, Rimini R, Nilsson P: Empirical bayes microarray ANOVA and grouping cell lines by equal expression levels. Stat Appl Genet Mol Biol. 2005, 4: Article7-[PMID: 16646860].PubMedGoogle Scholar
- Pan W: Incorporating biological information as a prior in an empirical bayes approach to analyzing microarray data. Stat Appl Genet Mol Biol. 2005, 4: Article12-[PMID: 16646829].PubMedGoogle Scholar
- Combat. [http://www.bu.edu/jlab/wp-assets/ComBat/Abstract.html].
- Kupfer P, Guthke R, Pohlers D, Huber R, Koczan D, Kinne RW: Batch correction of microarray data substantially improves the identification of genes differentially expressed in rheumatoid arthritis and osteoarthritis. BMC Med Genomics. 2012, 5: 23-[PMID: 22682473].View ArticlePubMedPubMed CentralGoogle Scholar
- Smyth GK: Limma: linear models for microarray data. Bioinformatics and Computational Biology Solutions using R and Bioconductor. Edited by: Gentleman R, Carey V, Dudoit S, Irizarry R, Huber W. 2005, New York: Springer, 397-420.View ArticleGoogle Scholar
- Shi L, Jones WD, Jensen RV, Harris SC, Perkins RG, Goodsaid FM, Guo L, Croner LJ, Boysen C, Fang H, Qian F, Amur S, Bao W, Barbacioru CC, Bertholet V, Chu TM, Collins PJ, Frueh FW, Fuscoe JC, Guo X, Han J, Herman D, Hong H, Kawasaki ES, Li QZ, Luo Y, Ma Y, Mei N, CaoXM, et al: The balance of reproducibility, sensitivity, and specificity of lists of differentially expressed genes in microarray studies. BMC, Bioinformatics. 2008, 9 (Suppl 9): S10-[PMID: 18793455 PMCID: 2537561].View ArticlePubMedGoogle Scholar
- Falcon S, Gentleman R: Using GOstats to test gene lists for GO term association. Bioinformatics. 2007, 23 (2): 257-8. [PMID: 17098774].View ArticlePubMedGoogle Scholar
- Nikitin A, Egorov S, Daraselia N, Mazo I: Pathway studio—the analysis and navigation of molecular networks. Bioinformatics. 2003, 19 (16): 2155-2157. [http://www.ncbi.nlm.nih.gov/pubmed/14594725]. [PMID: 14594725].View ArticlePubMedGoogle Scholar
- Wollbold J, Huber R, Pohlers D, Koczan D, Guthke R, Kinne RW, Gausmann U: Adapted Boolean network models for extracellular matrix formation. BMC Syst Biol. 2009, 3: 77.View ArticlePubMedPubMed CentralGoogle Scholar
- Guthke R, Möller U, Hoffmann M, Thies F, Töpfer S: Dynamic network reconstruction from gene expression data applied to immune response during bacterial infection. Bioinformatics. 2005, 21 (8): 1626-1634. [PMID: 15613398].View ArticlePubMedGoogle Scholar
- Toepfer S, Guthke R, Driesch D, Woetzel D, Pfaff M: Knowledge Discovery and Emergent Complexity in Bioinformatics, Volume 4366 of Lecture Notes in Computer Science. Edited by: Nowé A, Saeys Y, Westra R, Tuyls K. 2007, Springer Berlin Heidelberg, 119-130. [http://dx.doi.org/10.1007/978-3-540-71037-0_8].Google Scholar
- Hecker M, Lambeck S, Toepfer S, van Someren E, Guthke R: Gene regulatory network inference: data integration in dynamic models-a review. Biosystems. 2009, 96: 86-103. [PMID: 19150482].View ArticlePubMedGoogle Scholar
- Gustafsson M, Hörnquist M: Gene expression prediction by soft integration and the elastic net-best performance of the DREAM3 gene expression challenge. PLoS One. 2010, 5 (2): e9134-[PMID: 20169069].View ArticlePubMedPubMed CentralGoogle Scholar
- Clines GA, Mohammad KS, Bao Y, Stephens OW, Suva LJ, Shaughnessy JD Jr, Fox JW, Chirgwin JM, Guise TA: Dickkopf homolog 1 mediates endothelin-1-stimulated new bone formation. Mol Endocrinol. 2007, 21 (2): 486-98. [PMID: 17068196].View ArticlePubMedGoogle Scholar
- Bascom CC, Wolfshohl JR, Madisen L, Webb NR, Purchio AR, Moses HL, Coffey RJ Jr: Complex regulation of transforming growth factor beta 1, beta 2, and beta 3 mRNA expression in mouse fibroblasts and keratinocytes by transforming growth factors beta 1 and beta 2. Mol Cell Biol. 1989, 9 (12): 5508-5515. [PMID: 2586525].View ArticlePubMedPubMed CentralGoogle Scholar
- Kim SJ, Angel P, Lafyatis R, Hattori K, Kim KY, Sporn MB, Karin M, Roberts AB: Autoinduction of transforming growth factor beta 1 is mediated by the AP-1 complex. Mol Cell Biol. 1990, 10 (4): 1492-1497. [PMID: 2108318].View ArticlePubMedPubMed CentralGoogle Scholar
- Kiriyama T, Gillespie MT, Glatz JA, Fukumoto S, Moseley JM, Martin TJ: Transforming growth factor beta stimulation of parathyroid hormone-related protein (PTHrP): a paracrine regulator?. Mol Cell Endocrinol. 1993, 92: 55-62. [PMID: 8472867].View ArticlePubMedGoogle Scholar
- Spagnoli A, O’Rear L, Chandler RL, Granero-Molto F, Mortlock DP, Gorska AE, Weis JA, Longobardi L, Chytil A, Shimer K, Moses HL: TGF-beta signaling is essential for joint morphogenesis. J Cell Biol. 2007, 177 (6): 1105-1117. [PMID: 17576802].View ArticlePubMedPubMed CentralGoogle Scholar
- Ruebel KH, Leontovich AA, Tanizaki Y, Jin L, Stilling GA, Zhang S, Coonse K, Scheithauer BW, Lombardero M, Kovacs K, Lloyd RV: Effects of TGFbeta1 on gene expression in the HP75 human pituitary tumor cell line identified by gene expression profiling. Endocrine. 2008, 33: 62-76. [PMID: 18401765].View ArticlePubMedGoogle Scholar
- Javelaud D, Pierrat MJ, Mauviel A: Crosstalk between TGF-βand hedgehog signaling in cancer. FEBS Lett. 2012, 586 (14): 2016-2025. [PMID: 22609357].View ArticlePubMedGoogle Scholar
- Roberts AB, Russo A, Felici A, Flanders KC: Smad3: a key player in pathogenetic mechanisms dependent on TGF-beta. Ann N Y Acad Sci. 2003, 995: 1-10. [PMID: 12814934].View ArticlePubMedGoogle Scholar
- Sterling JA, Guelcher SA: Bone structural components regulating sites of tumor metastasis. Curr Osteoporos Rep. 2011, 9 (2): 89-95. [PMID: 21424744].View ArticlePubMedPubMed CentralGoogle Scholar
- Aomatsu K, Arao T, Sugioka K, Matsumoto K, Tamura D, Kudo K, Kaneda H, Fujita Y, Shimomura Y, Nishio K, Tanaka K: TGF-beta induces sustained upregulation of SNAI1 and SNAI2 through Smad and non-Smad pathways in a human corneal epithelial cell line. Invest Ophthalmol Vis Sci. 2011, 52 (5): 2437-2443. [PMID: 21169525].View ArticlePubMedGoogle Scholar
- Vallabhapurapu S, Karin M: Regulation and function of NF-kappaB, transcription factors in the immune system. Annu Rev Immunol. 2009, 27: 693-733. [PMID: 19302050].View ArticlePubMedGoogle Scholar
- Brown KD, Claudio E, Siebenlist U: The roles of the classical and alternative nuclear factor-kappaB pathways: potential implications for autoimmunity and rheumatoid arthritis. Arthritis Res Ther. 2008, 10 (4): 212-[PMID: 18771589].View ArticlePubMedPubMed CentralGoogle Scholar
- Nieto MA, Sargent MG, Wilkinson DG, Cooke J: Control of cell behavior during vertebrate development by Slug, a zinc finger gene. Science. 1994, 264 (5160): 835-839. [PMID: 7513443].View ArticlePubMedGoogle Scholar
- Mo R, Freer AM, Zinyk DL, Crackower MA, Michaud J, Heng HH, Chik KW, Shi XM, Tsui LC, Cheng SH, Joyner AL, Hui C: Specific and redundant functions of Gli2 and Gli3 zinc finger genes in skeletal patterning and development. Development. 1997, 124: 113-123. [PMID: 9006072].PubMedGoogle Scholar
- Karaplis AC, Luz A, Glowacki J, Bronson RT, Tybulewicz VL, Kronenberg HM, Mulligan RC: Lethal skeletal dysplasia from targeted disruption of the parathyroid hormone-related peptide gene. Genes Dev. 1994, 8 (3): 277-289. [PMID: 8314082].View ArticlePubMedGoogle Scholar
- Funk JL, Cordaro LA, Wei H, Benjamin JB, Yocum DE: Synovium as a source of increased amino-terminal parathyroid hormone-related protein expression in rheumatoid arthritis. A possible role for locally produced parathyroid hormone-related protein in the pathogenesis of rheumatoid arthritis. J Clin Invest. 1998, 101 (7): 1362-1371. [PMID: 9525978].View ArticlePubMedPubMed CentralGoogle Scholar
- Amizuka N, Henderson JE, White JH, Karaplis AC, Goltzman D, Sasaki T, Ozawa H: Recent studies on the biological action of parathyroid hormone (PTH)-related peptide (PTHrP) and PTH/PTHrP receptor in cartilage and bone. Histol Histopathol. 2000, 15 (3): 957-970. [PMID: 10963138].PubMedGoogle Scholar
- Hartmann C, Tabin CJ: Wnt-14 plays a pivotal role in inducing synovial joint formation in the developing appendicular skeleton. Cell. 2001, 104 (3): 341-351. [PMID: 11239392].View ArticlePubMedGoogle Scholar
- Maioli E, Fortino V, Torricelli C, Arezzini B, Gardi C: Effect of parathyroid hormone-related protein on fibroblast proliferation and collagen metabolism in human skin. Exp Dermatol. 2002, 11 (4): 302-310. [PMID: 12190938].View ArticlePubMedGoogle Scholar
- Guo X, Day TF, Jiang X, Garrett-Beal L, Topol L, Yang Y: Wnt/beta-catenin signaling is sufficient and necessary for synovial joint formation. Genes Dev. 2004, 18 (19): 2404-2417. [PMID: 15371327].View ArticlePubMedPubMed CentralGoogle Scholar
- Haque T, Nakada S, Hamdy RC: A review of FGF18: its expression, signaling pathways and possible functions during embryogenesis and post-natal development. Histol Histopathol. 2007, 22: 97-105. [PMID: 17128416].PubMedGoogle Scholar
- Presta M, Andrés G, Leali D, Ronca R, Dell’Era P: Inflammatory cells and chemokines sustain FGF2-induced angiogenesis. Eur Cytokine Netw. 2009, 20 (2): 39-50. [PMID: 19541589].PubMedGoogle Scholar
- Bramlage CP, Häupl T, Kaps C, Ungethüm U, Krenn V, Pruss A, Strutz F, Burmester GR, Müller GA: Decrease in expression of bone morphogenetic proteins 4 and 5 in synovial tissue of patients with osteoarthritis and rheumatoid arthritis. Arthritis Res Ther. 2006, 8 (3): R58-[PMID: 16542506].View ArticlePubMedPubMed CentralGoogle Scholar
- Huber LC, Distler O, Tarner I, Gay RE, Gay S, Pap T: Synovial fibroblasts: key players in rheumatoid arthritis. Rheumatology. 2006, 45 (6): 669-675. [PMID: 16567358].View ArticlePubMedGoogle Scholar
- The pre-publication history for this paper can be accessed here:http://www.biomedcentral.com/1755-8794/7/40/prepub
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited.