Novel factors in the pathogenesis of psoriasis and potential drug candidates are found with systems biology approach.

Psoriasis is a multifactorial inflammatory skin disease characterized by increased proliferation of keratinocytes, activation of immune cells and susceptibility to metabolic syndrome. Systems biology approach makes it possible to reveal novel important factors in the pathogenesis of the disease. Protein-protein, protein-DNA, merged (containing both protein-protein and protein-DNA interactions) and chemical-protein interaction networks were constructed consisting of differentially expressed genes (DEG) between lesional and non-lesional skin samples of psoriatic patients and/or the encoded proteins. DEGs were determined by microarray meta-analysis using MetaOMICS package. We used STRING for protein-protein, CisRED for protein-DNA and STITCH for chemical-protein interaction network construction. General network-, cluster- and motif-analysis were carried out in each network. Many DEG-coded proteins (CCNA2, FYN, PIK3R1, CTGF, F3) and transcription factors (AR, TFDP1, MEF2A, MECOM) were identified as central nodes, suggesting their potential role in psoriasis pathogenesis. CCNA2, TFDP1 and MECOM might play role in the hyperproliferation of keratinocytes, whereas FYN may be involved in the disturbed immunity in psoriasis. AR can be an important link between inflammation and insulin resistance, while MEF2A has role in insulin signaling. A controller sub-network was constructed from interlinked positive feedback loops that with the capability to maintain psoriatic lesional phenotype. Analysis of chemical-protein interaction networks detected 34 drugs with previously confirmed disease-modifying effects, 23 drugs with some experimental evidences, and 21 drugs with case reports suggesting their positive or negative effects. In addition, 99 unpublished drug candidates were also found, that might serve future treatments for psoriasis.

Introduction

Psoriasis is a multifactorial inflammatory skin disease. A recent systematic review reported a prevalence from 0% (Taiwan) to 2.1% (Italy) in children and from 0.91% (United States) to 8.5% (Norway) in adults.[1] Genetic predisposition and environmental factors are both important in disease etiology. Several genome-wide association studies have been carried out and until now 36 susceptibility loci have been identified.[2] Environmental triggers are also reported such as drugs, smoking, mental stress, skin injury, Streptococcal infection, hormonal changes etc.[3] Psoriasis is an immune-mediated disease. Important immune cells and cytokines have been identified in disease pathogenesis such as IL6, IL17A, TNF etc.[4] Autoimmune basis for chronic inflammation is supposed, although no consistent antigen has been found. Patients with psoriasis have higher risk for metabolic syndrome, and risk increases with disease severity. Both diseases have immunological basis with common cytokines and genetic risk loci like CDKAL1.[5] Keratinocyte hyperproliferation is present in lesional phenotype and is responsible for scale formation. Keratinocyte differentiation markers like keratin 1 and keratin 10 are downregulated and parakeratosis (keratinocytes with nuclei in the stratum granulosum) is also present.[3]

Psoriasis is one of the most studied skin diseases. By now more than 34000 hits are available in PubMed for the keyword “psoriasis” and the number is increasing. No spontaneous psoriasis-like skin disease is known in animals. Induced mouse models are available which are similar, but not the same as psoriasis in human.[6] Therefore drug discovery is difficult in such models what makes in silico analysis more essential. “Omics” data gives the opportunity to examine the disease with systems biology approach.

Stationary changes in gene expression are responsible for fixing phenotypes such as lesional skin areas in psoriasis. Several microarray studies have been carried out to characterize gene expression in healthy and psoriatic skin samples (Table 1). Microarray meta-analysis gives the opportunity to evade biological, regional, and study design-caused variation between studies.[7] Network analysis is a novel and highly developing area of systems biology. Considering gene expression data it is possible to explain alterations in intracellular processes with the analysis of protein-protein and protein-DNA (or gene regulatory) interaction networks. These networks consist of proteins and/or regulated genes as nodes and undirected or directed edges between them. Network centralities like degree or stress are suitable for ranking nodes. Total edge number belong to one node equals its degree in undirected networks. Nodes have in- and out-degrees based on edge directions in directed networks. Degree distribution follows a scale-free power law distribution in biological networks. This fact indicates that highly connected vertices have a large chance of occurring. Nodes with highest degree are called hubs and are essential in network stability.[8] Stress centrality indicates the number of shortest paths (from all shortest paths between any two nodes in the network) passing through the given node thus the capability of a protein for holding together communicating nodes.[9] Interconnecting nodes make up network motifs. Several, such as feed-forward or bifan motif are significantly enriched in biological networks compared to random networks. These elements have important role in network dynamics.[10]

Table 1

Study information and QC measure summary.

	Study	MIAME	GEO ID	Platform/Chip	NL	L	IQC	EQC	CQCg	CQCp	AQCg	AQCp	Rank
1	Gudjonsson et al.[110]	Available	GSE13355	GPL570/Affymetrix HU133 Plus 2.0	54	53	4.18	4	307.65	307.65	95.2	292.19	2.17
2	Yao et al.[111]	Available	GSE14905	GPL570/Affymetrix HU133 Plus 2.0	27	32	5.58	4	307.65	307.65	81.32	185.34	2.67
3	Zaba et al.[112]	Available	GSE11903	GPL571/Affymetrix HU133A 2.0	15	12	7.34	3	307.65	307.65	79.24	260.95	2.75
4	Suarez-Farinas et al.[113]	Available	GSE30999	GPL570/Affymetrix HU133 Plus 2.0	79	80	0.86*	4	307.65	307.65	33	193.93	3.67
5	Reischl et al.[114]	Available	GSE6710	GPL96/Affymetrix HU133A	12	12	2.7	4	307.65	271.23	40.3	118.68	3.92
6	Johnson-Huang et al.[115]	Available	GSE30768	GPL571/Affymetrix HU133A 2.0	1	4	Excluded by Array Quality Metrics package

MIAME information was available for all study Studies were downloaded from Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/). All studies were carried out on Affymetrix platforms. Lesional and Non-Lesional sample count is shown. Stars in table indicate non-statistical significance of QC measures. Study no 6 was already excluded by sample filtering by arrayQualityMetrics. Other studies had high quality and no outlier study was present. IQC: Internal Quality index, EQC: External Quality index, CQCg and CQCp: Consistency Quality Control indexes, AQCg and AQCp: Accuracy Quality Control indexes, NL: non-lesional sample count, L: lesional sample count.

We hypothesized that it could be possible to find novel elements of psoriasis pathogenesis with detailed analysis of precisely constructed networks. Network motif enrichment caused by changes in gene expression could have important role in disease development and sustainment. It could be also possible to detect potential drug candidates by analyzing chemical–protein networks. Thus our goal was to construct reliable but yet detailed protein-protein, protein-DNA, merged (containing both protein-protein and protein-DNA interactions) and chemical-protein interaction networks consisting of differentially expressed genes (DEG) between lesional and non-lesional skin samples and/or the coded proteins. Detailed analysis of these networks could help us to reveal novel players in disease pathomechanism and to identify network motifs and sub-networks with the ability to sustain lesional phenotype.

Methods

Microarray Meta-analysis

Six microarray studies examining lesional and non-lesional skin biopsy samples of psoriatic patients were found in Gene Expression Omnibus (GEO) (Table 1). “Minimum Information About a Microarray Experiment” (MIAME) was available for each study. Only non-lesional and lesional samples from affected individuals were used for analysis, samples from healthy people were excluded. Raw.CEL files were downloaded and quality of each sample was assessed with the R package arrayQualityMetrics.[11] This package defines sample quality with 5 different methods and generates plots for outlier detection. A sample was excluded if it was obviously an outlier in at least 1 measure or had borderline values in at least 2 measures (analysis results are in Dataset S1 compressed file; outliers and argument of exclusion is listed in Table S1). Raw data normalization of remaining samples was carried out with the R package Easy Microarray data Analysis (EMA).[12] GCRMA normalization method was used and probe sets with expression level below 3.5 were discarded. Probe set with the highest interquartile range (IQR) was chosen for common HUGO Gene Nomenclature Committee (HGNC) gene identifiers. Original findings were confirmed with published statistics. For this EMA was used after GCRMA normalization. More DEGs were found in some cases, which might be caused by the pre-filtering process with arrayQualityMetrics (Table S2). The R package MetaQC was used for filtering out low quality studies.[13] The fifty most prevalent gene set were chosen with the software Gene Set Enrichment Analysis (GSEA) and used for external quality control (EQC) score calculation.[14] GSEA was carried out for each study with the following settings: 1000 permutations; minimum set size was 5 and the gene set database was c2.all.4.0.symbols. The resultant study-level p values of a gene set were combined with Fisher's combined probability test. The fifty gene sets with the lowest meta-analysis p value were chosen as input for EQC score calculation. C2.all.4.0.symbols gene set database was chosen as input for consistency quality control (CQCp) value calculation. GSEA input expression matrices contained gene IDs that were present in all studies after EMA filtering. MetaDE package was used to determine DEGs in lesional samples compared to non-lesional ones.[15] DEG p value in individual studies was calculated by two sample T test with unequal variances. Fisher's combined probability test was chosen for meta-analysis statistical method.[16] Fold change of gene expression was given by the ratio between geometrical means of gene expression in lesional and non-lesional samples.[17] Genes with false discovery rate (FDR) less than 0.001 and with fold change higher than 1.5 or less than −1.5 were accepted as DEGs.

Construction of protein-protein, protein-DNA and chemical-protein interaction networks

STRING database 9.0 was used as resource for protein-protein interactions (PPI).[18] Both directed and undirected networks were created by selecting all interactions between DEG – coded proteins in downloaded raw data. Interaction confidence score cutoff was 900 (“highest confidence” group) in case of undirected and 800 (containing a part of “high confidence” and all “highest confidence” interactions) in case of directed interactions. Only directed interactions with “activation” or “ptmod” actions were used. Chemical-protein interactions between potential drugs, intra- and extracellular compounds and DEG-coded proteins were collected from STITCH database 3.1.[19] The way of interaction confidence score calculation is the same in this database as in STRING thus interactions with the described confidence score cutoff values were selected for network construction. Protein-DNA interaction (PDI) network consisting of DEGs and DEG-coded transcription factors (TF) was created using cis-Regulatory Element Database (CisRED).[20] Regulatory element motifs with were collected from DEG promoter regions. Motifs were coupled with TFs or TF complexes using TRANSFAC and JASPAR databases.[21], [22] Motifs without respective TFs were excluded. Merged DEG-derived network containing PPI and PDI interactions and a network containing only DEG-coded TFs were also generated. Complete PPI, PDI, merged, TF-TF and chemical-protein interaction networks were created for controls using all available interactions in databases with the same statistical threshold as in DEG-derived network construction.

General network analysis, identification of central nodes and motif detection

General network analysis and node centrality value calculation were carried out with NetworkAnalyzer Cytoscape plugin.[23] Isolated nodes and node groups (without connection with the main PPI network) were deleted from graph in order to evade false results. Curve fitting on node degree and stress value distributions was done with MATLAB Curve Fitting Tool (MATLAB R2012b, The Mathworks Inc., Natick, MA). Curve of power law distribution was assessed with Trust-Region algorithm. Goodness of fitting was assessed by R-square and corrected R-square values which prove power law distribution of these node centralities (Table 2). As power law distribution is asymmetric with a long tail, nodes with centralities above average cannot be assessed using arithmetic mean. A variable with a power-law distribution has a probability of taking a value following the function , where is constant. First moment (mean value) of a power-law distributed quantity equals:

Table 2

Results of node centrality analysis.

Network	Centrality	Curve	Cutoff	R-square	Adjusted R-square
PPI Undirected	Degree	0.8555◯-1.649	27.21	0.9957	0.9956
	Stress	47.1◯-0.8034	427072.25	0.9795	0.9793
PPI Directed	In-Degree	0.5925◯-1.808	5.100152	0.9969	0.9968
	Out-Degree	0.5462◯-1.759	23.493461	0.9983	0.9983
	Stress	15.34◯-0.961	8504.103	0.8753	0.8748
PDI	Out-Degree	13280◯-1.367	287.20865	0.9252	0.9002
CPI Undirected	Degree	0.8314◯-3.168	14.761	1	1
	Stress	2.41e1-2.432	6.63	0.9811	0.9811
CPI Directed	Out-Degree	0.7859◯-2.132	8.5757576	1	1

Distribution of node centrality values were assessed by curve fitting. Curve equations, goodness of fit (R-square and adjusted R-square) and the resultant cutoff values are shown. CPI: chemical – protein interaction network.

Second moment (variance) of a power-law distributed quantity equals:

The sum of first and second moment (mean value and variance) was used as cutoff for centralities with distribution exponent . Expression of variance becomes infinite, when , thus only first moment (mean value) was used as cutoff for centralities with distribution exponent. [24] Expression of mean value becomes infinite, if . In this case weighted mean was used to assess cutoff with the following formula:

As bidirectional connections are available in undirected PPI network, stress centrality is independent from edge directions thus both degree and stress had to be above cutoff for central protein selection. As directed networks contain unidirectional interactions, low stress values (i. e. low number of shortest paths cross through the node) can be caused by the dominance of incoming (in-degree) or outgoing (out-degree) interactions. Important nodes with high in-degree or out-degree can still have low stress centrality thus either out-degree or in-degree or stress had to be above cutoff in directed PPI network. As TFs have mainly outgoing interactions, out-degree was used for TF prioritization. Similarly to PPI networks degree and stress had to be above cutoff in undirected chemical - protein interaction network. Drugs with more targets in DEG-derived PPI-networks may have bigger disease modifying effect thus out-degree had to be above cutoff in directed chemical – protein interaction network for drug prioritization (Table 2).

NetMODE software was used for network motif statistical analysis. Frequency of 3 or 4 node motifs in DEG-derived and complete control networks were compared with 1000 random graphs. Local constant switching mode was used for edge switching method during random network generation. NetMODE p value indicates the number of random networks in which a motif occurred more often than in the input network, divided by total number of random networks. was used as cutoff.[25] Respective sub-networks of enriched motifs were identified with NetMatch Cytoscape plugin.[26] jActiveModules and ClusterONE were used for network module and protein complex detection. ClusterONE analysis was carried out with minimum cluster size of 3 with unweighted edges and default advanced parameters. jActiveModules considers gene expression for module search. Input gene expression values have to be between 0 and 1 so normalized expression values got with EMA were scaled between these numbers.[27], [28] Functional description of node groups was done with BinGO (“Biological function” GO terms were selected, FDR<0.001 was used for term enrichment).[29]

Results

Detection of DEGs with microarray meta-analysis

In order to get reliable data about gene expression in lesional psoriatic skin samples microarray meta-analysis was carried out. The study by Johnson-Huang et al. was already excluded after sample quality analysis with arrayQualityMetrics package, because at least two samples from one phenotype group are needed for MetaQC analysis and only one non-lesional sample remained after sample filtering. The overall quality of each study was assessed by MetaQC.[13] The software calculated six quality control (QC) measures then created principal component analysis (PCA) biplot and standardized mean rank summary (SMR) score to help in the identification of problematic studies. It was described by authors, that if a study is on the opposite side of arrows in the PCA biplot and has large SMR scores, it's strongly suggested to be excluded from meta-analysis. In contrary, if a study is on the same side of arrows in the PCA biplot and has small SMR scores, it should be included. All five studies were defined as usable based on quality values (Table 1, Figure 1). DEGs were identified by MetaDE.[15] 2307 upregulated and 3056 downregulated genes were found in lesional skin samples compared to non-lesional ones (Table S3). The relatively high number of DEGs can be the result of filtering out low quality samples, which could increase variance and using lower fold change cutoff values than in original studies. DEGs were used for network construction.

Figure 1

PCA biplot numbers on PCA biplot represents studies in Table 1.

Study number placed opposite to quality measure axes are of low quality and should be excluded. No outlier study was detected.

General Network analysis

Undirected and directed PPI networks with DEG – coded proteins, directed PDI networks with DEG – coded TFs and regulated DEGs and merged directed networks containing both PPIs and PDIs were created. A TF-TF network consisting of DEG-coded TFs was also generated. The Cytoscape plugin NetworkAnalyzer calculated main network properties for both DEG-derived and control complete networks (Table 3). DEG – derived networks had higher diameter (i. e. the length of the longest shortest path in the network) and average shortest path length than control full networks. This may be caused by the inverse correlation of node degree and fold change.[30] Nodes with lower fold change has higher degree. Genes with fold change under cutoff are filtered out from DEG derived networks (between red lines on Figure 2). The remaining nodes has smaller average degree, therefore connectivity of the network is lower resulting in higher diameter and average shortest path length value.

Table 3

Results of general network analysis.

Network	Nodes	Edges	Diameter	Average shortest path
PPI Undirected	1614 (9412)	5156 (55039)	14 (12)	4.79 (4.45)
PPI Directed	464 (4040)	815 (13377)	14 (12)	5.26 (4.35)
PDI	2840 (15839)	6398 (123210)	10 (7)	3.69 (3.029)

DEG derived and control networks has similar attributes, but average shortest path length and network diameter is lower in DEG derived networks, which can be explained by lower connectivity (Figure 2). Values for control networks are in brackets.

Figure 2

Degree-Fold Change relationship.

Nodes with higher degree has lower fold change of gene expression in all network types. Genes between red lines have higher average degree and are filtered out from network analysis. Remaining nodes in DEG-derived networks have lower average degree and connectivity.

Determination of hubs in DEG-derived networks

Most important nodes of DEG-derived networks were determined using degree and/or stress centralities (Table 2, full list of nodes and centralities is in Table S4). Numerous already published psoriasis-associated protein-coding genes were found (Table 4). CCNA2, FYN and PIK3R1 proteins are present in top rated hubs in undirected PPI network and are yet unpublished in association with the disease. CCNA2 have role in mitosis regulation.[31] FYN is important in interferon gamma (IFN gamma) signaling, while PIK3R1 is important in insulin-stimulated glucose uptake.[32], [33] FYN could be found in jActiveModules cluster with the 2nd highest score while PIK3R1 were found in cluster with the 3rd highest score (Figure S1, S2). Taking account BinGO results these clusters are responsible for signaling and for immune regulation as well (Table S5). A highly connected chemokine-chemokine receptor cluster was also found with ClusterONE analysis (Figure S3). Central nodes in directed and undirected PPI networks showed overlap (Table 4). CTGF is in top ranked proteins and yet not associated with psoriasis. CTGF is responsible for fibrosis downstream of TGFβ signaling. Downregulation of CTGF by psoriasis-associated cytokines INFγ and TNFα is already published.[34]

Table 4

Top rated nodes in DEG-derived networks.

PPI Undirected		PPI Directed		PDI
Name	Fold change	Name	Fold change	Name	Fold change
IL8	67.31113193	IL8	67.31113193	TFDP1	4.612130627
CCNB1	11.13277565	BIRC5	9.309154577	MECOM	1.705869235
BIRC5	9.309154577	MMP1	7.446458555	AR	−1.649992095
STAT1	9.038900879	SOD2	7.198087989	NF1	−1.707954442
CCNA2	8.737535122	IL1B	4.293906976	MEF2A	−1.738635445
CXCR4	5.109553129	STAT3	3.965626652
IL1B	4.293906976	MMP9	3.661047085
MAPK14	4.152927326	SOCS3	3.315643007
STAT3	3.965626652	HMOX1	3.207443671
MMP9	3.661047085	CCL2	2.896844503
LCK	3.609090653	BAX	1.9009731
AURKB	2.493884913	ICAM1	1.722246429
MAPK1	1.820524831	CD69	1.721780507
MYC	1.690987073	MYC	1.690987073
NFKB1	1.636019496	CD86	1.676295675
PCNA	1.623673041	CD28	1.640633244
CDKN1A	1.583889601	NFKB1	1.636019496
HDAC1	1.57828429	EGFR	−1.607280925
CYP1A1	−1.595883159	CTNNB1	−1.648110677
EGFR	−1.607280925	FN1	−1.75413351
CREBBP	−1.626480892	EDN1	−1.836157927
CTNNB1	−1.648110677	SP1	−1.923552267
FN1	−1.75413351	CTGF	−2.037178621
FYN	−1.849385591	NFATC1	−2.187942784
SP1	−1.923552267	IRS1	−2.277490062
SMAD4	−1.95145712	INS−IGF2	−2.33005624
INS-IGF2	−2.33005624	CCND1	−2.341844947
CCND1	−2.341844947	FOS	−2.362430819
FOS	−2.362430819	PPARG	−2.556455049
PPARG	−2.556455049	BCL2	−2.632996792
BCL2	−2.632996792	F3	−3.835078706
PIK3R1	−2.955639724	LEP	−6.266827433

Central proteins with centrality value(s) above cutoff are listed. Fold change between gene expression in lesional and non-lesional samples are also shown. Proteins with bold characters are yet non-published in terms of psoriasis.

PDI network contained DEG-coded TFs and regulated DEGs as nodes and directed edges pointing from the TFs to the regulated genes. TFs were ranked using out-degree centrality. Androgen receptor (AR) and TFDP1 were the highest ranked nodes. AR is a TF, regulating genes that have immunological functions and role in carbohydrate metabolism.[35], [36] TFDP1 controls cell cycle progression and is yet not associated with psoriasis.[37] BinGO analysis of TFDP1-regulated genes prove its central role in cell cycle activation (Table S5). MECOM and MEF2A are TFs above centrality cutoff and yet not associated with psoriasis. MECOM have role in cell proliferation and is associated with chronic myeloid leukemia.[38] MEF2A is responsible for the insulin dependent glucose transporter GLUT4 expression and is downregulated in insulin deficient diabetes mellitus.[39]

Motif analysis in DEG-derived networks

Motifs consisting of 3 or 4 nodes were analyzed in directed DEG-derived and control networks as well (Table 5, Figure 3). Analysis found motifs which were enriched in directed DEG-derived but were absent in control networks or vice versa. Some were already generally described in biological systems like convergent (no. 36), divergent (no. 6) and bifan (no. 204) motifs, but yet non-examined ones were detected like motif no. 924 in directed PPI networks, no. 332 in TF-TF networks and no. 6356 in merged networks etc. Cause of missing convergent, divergent and bifan motifs in DEG derived directed PPI or PDI networks compared to control was not investigated as uncertainty is present about the role of these network motifs in biological systems.[10] Identifying nodes making up motif no. 924 resulted in the high occurrence of central proteins found before. These proteins were associated with the immune system and carbohydrate metabolism. Motif 332 is enriched in the TF network of lesional skin. This motif is based on the TFDP1–AR reciprocal regulation. Importance of these TFs is already mentioned.

Table 5

Summary of network motif analysis.

	PPI directed		PDI		PDI +PPI
Motif no.	Psoriasis	Full	Psoriasis	Full	Psoriasis	Full
6 (divergent)	0.705	0.031	0.168	0.974	0.908	0.952
36 (convergent)	0.997	0.972	0.826	0.023	0.083	0.045
38 (feed-forward)	0	0	0.073	0.978	0.941	0.998
98 (feedback)	0.329	0.242	0.518	0.233	0.064	0.046
204 (bifan)	0.255	0	0.483	0.082	0.872	0.041
332	0.958	0.162	0.042 (TF network)	0.838 (TF network)	0.41	0.067
924	0.007	0.292	N/A	0.305	0.794	0.17
6356	0.025	0.02	N/A	0.916	0.001	0.512

Numbers are p values of motif enrichment compared to 1000 random networks. Values with bold characters are below 0.05 and thus significant. Significant enrichment was only found in TF-TF networks in case of motif no. 332. Network motif pictures are in Figure 3.

Figure 3

Network motifs with 3 or 4 nodes.

Analysis results of the respective motif can be found in Table 5.

An interesting result of motif analysis is the enrichment of feedback loops containing 3 nodes in merged networks compared to separate ones and the enrichment of motif no. 6356 in DEG-derived merged network compared to control. Motif no. 6356 consist of a positive feedback loop and all nodes of the loop are controlled by another separated node like IL1B or AR.

Controller sub-network construction

Both lesional and non-lesional skin areas can be found on patients at the same time. We wanted to highlight nodes which may be important in the “all or none” switch in lesional skin areas and sustain this phenotype for a long time. It has been argued that hubs in intracellular regulatory networks are enriched with either positive or negative regulatory links and cause much more positive feedback loops than negative ones.[40] It is also proven that positive feedback loops have fundamental role in maintaining autoimmune and autoinflammatory disease states.[41] Enrichment of motif no. 6356 consisting of a positive feedback loop with all nodes controlled by a separated one also suggests central role of positive feedback loops in lesional skin which may be activated by important central proteins like AR or IL1B. This is published that in biological systems interlinked slow and fast positive feedback loops allow systems to convert graded inputs (like several environmental and genetic factors in a psoriatic individual) into decisive all or none outputs (like lesional skin phenotype).[42] Transcriptional regulation needs time so we hypothesized that slow positive feedback loops may consist of at least one gene regulatory interaction. Fast loops may consist of only PPIs. Transcriptional changes of nodes in these loops may be able to sustain the “switched on” state.

In order to find most important slow and fast feedback loops containing 2, 3 or 4 nodes, a merged PPI and PDI network was constructed from proteins with centralities above cutoff value. All feedback loops were identified with NetMatch. A positive feedback loop was selected if and only if expression of all nodes changes in the direction of sustaining or suppressing the activity of the loop and “activation” or “inhibition” properties of all edges were proven by publications. Expression of all nodes was downregulated in two loops needed for carbohydrate metabolism: the INS-IGF2-EDN1-LEP-INS-IGF2 and the LEP-PPARG-INS-IGF2-LEP loop. The IL1B-NFKB1-CCL2-IL1B loop contained only upregulated nodes and has role in inflammation (Figure 4). The remaining loops contained inflammation and metabolism-related nodes as well. These may be key components in the metabolic-inflammatory interplay in the pathomechanism of psoriasis. “Slow” positive feedback loops containing gene regulatory interactions and “fast” loops containing only PPIs were also found. All positive feedback loops had common nodes, thus a merged network was generated containing interlinked slow and fast positive feedback loops (Figure 4). Transcriptional changes of all nodes and influence of all edges supported the sustainment of lesional phenotype in this sub-network. Boolean analysis of the resultant controller network was also performed. Nodes with downregulated expression got value of 0 and nodes with upregulated expression got value of 1. Future state of nodes was set based on interactions (Table 6). The output boolean values were the same as the input state values which prove the role of the controller network in the sustainment of present (lesional) phenotype. Chemical - protein interaction analysis further prove the importance of controller network.

Figure 4

Positive feedback loops and the merged controller sub-network in lesional psoriatic skin.

Individual positive feedback loops with 2, 3 or 4 nodes are shown. Node color is blue if the gene expression is decreased and red if increased. Merged controller sub-network is shown on the top. Node color is proportional with fold change. red line: gene regulatory interaction; blue line: protein-protein interaction; arrow-headed line: activation; bar-headed line: inhibition

Table 6

Boolean analysis of controller network.

	Input state	Relation	Future state(*)
NFATC*  =  FOS	0	0	0
FOS*  =  EDN1	0	0	0
EDN1*  =  NFATC1 and INS-IGF2 and LEP	0	0 and 0 and 0	0
INS-IGF2*  =  PPARG and LEP	0	0 and 0	0
LEP*  =  EDN1 and INS-IGF2	0	0 and 0	0
PPARG*  =  INS-IGF2 and LEP and AR	0	0 and 0 and 0	0
AR*  =  not (IL8 and NFKB1)	0	not (1 and 1)	0
STAT3*  =  not AR	1	not 0	1
IRF1*  =  STAT3	1	1	1
IL8*  =  not PPARG; STAT3 and IRF1 and NFKB1	1	not 0; 1 and 1 and 1	1
IL1B*  =  CCL2	1	1	1
NFKB1*  =  not AR; IL1B	1	not 0; 1	1
CCL2*  =  NFKB1 and IL1B	1	1 and 1	1

Logical relations can be seen in the first and third column. Input and future state of network is stationary.

Analysis of chemical-protein interaction networks

Undirected and directed chemical-protein interaction networks were constructed using STITCH database, which contains interactions between proteins and chemical compounds (internal non-protein substances, drugs and environmental substances).[19] Drugs or potential drugs were filtered out from chemicals and ranked by degree and stress centrality in case of undirected and out degree centrality in case of directed networks (Table S4).

Top ranked drugs were grouped into Anatomical Therapeutic Chemical (ATC) classes (Table 7).[43] KEGG DRUG was used for classification.[44] Results show a big overlap between undirected and directed network analysis. Best rated drugs consisted of retinoic acid, cholecalciferol, costicosteroids, methotrexate, sirolimus and tacrolimus, which can be already found in psoriasis guidelines and large clinical trials have proved their effectiveness.[45]

Table 7

Published Drugs.

ATC Class	Drugs
STUDIES AVAILABLE
Retinoids for topical use in acne	retinoic acid
Corticosteroids	dexamethasone, hydrocortisone, corticosterone, prednisolone
H2 receptor antagonists	cimetidine
Immunosupressants	sirolimus, tacrolimus
Antiinflammatory and antirheumatic drugs	indomethacin
Blood glucose lowering drugs excl. insulines	metformin, troglitazone, rosiglitazone, pioglitazone
Intestinal anti-inflammatory agents	sulfasalazine
Vitamins	cholecalciferol, folic acid
Antimycobacterials	rifampicin
Mineral supplements	selenium
Antifungals for topical use	salicylic acid
Antineoplastic agents	5-fluorouracil, methotrexate, paclitaxel, cycloheximide
Cardiac stimulants excl. cardiac glycosides	epinephrine-bitartrate, norepinephrine
Lipid-modifying agents, plain	simvastatin, atorvastatin-calcium
Calcium channel blockers	nifedipine
Psychoanaleptics	caffeine
Thyroid therapy	Liothyronine
Drugs for obstructive airway diseases	theophylline
N/A	berberine, curcumin, triptolide
EXPERIMENTAL EVIDENCE
Topical products for joint and muscular pain	capsaicin
Respiratory system	N-acetyl-L-cysteine
Antineoplastic agents	Velcade, celecoxib
Hormone antagonists and related agents	tamoxifen
Cardiac stimulants excl. cardiac glycosides	isoproterenol
Liver therapy	glycyrrhizinic acid
Antiinfectives and antiseptics, excl. combinations with corticosteroids	arsenic
Beta blocking agents	propranolol
Lipid-modifying agents, plain	clofibrate, bezafibrate, fluvastatin, pravastatin
Blood glucose lowering drugs excl. insulines	ciglitazone
N/A	N-ethylmaleimide, baicalein, apigenin, SB 202190, monensin, rolipram, eflornithine, calphostin C, trichostatin A, rottlerin
CASE REPORTS
Antivirals for systemic use	ritonavir
Antiinflammatory and antirheumatic drugs	diclofenac, ibuprofen, aspirin
Antigout preparations	colchicine
Antiprotozoals	chloroquine
Ophtalmologicals	atropine
Antineoplastic agents	cytarabine-hydrochloride, doxorubicin, cysplatin, imatinib, docetaxel, gefitinib
Cardiac stimulants excl. cardiac glycosides	phenylephrine
Antiadrenergic agents, centrally acting	clonidine
Agents acting on the renin-angiotensin system	captopril, losartan
Anaesthetics	lidocaine
Psycholeptics	olanzapine
Psychoanaleptics	fluoxetine
Other nervous system drugs	nicotine

Psoriasis studies are available for numerous potential drugs with high centralities. “Blood glucose lowering drugs” are promising drug candidates. The biguanide metformin is associated with reduced psoriasis risk in a population based case control study.[46] Many studies are available about “Thiazolidinedione” group. A recent meta-analysis showed significant decrease in Psoriasis Area and Severity Index (PASI) scores compared to placebo in case of pioglitazone and non-significant improvement in PASI 50/70 in case of rosiglitazone.[47] Troglitazone normalized histological features in psoriasis models and the lesional phenotype in a small clinical trial.[48] The “HMG CoA reductase inhibitor” drug simvastatin was effective in a pilot study, although atorvastatin in the same class showed only a non-significant improvement in a different study.[49], [50] Salicylic acid has antifungal effects and it's used as adjuvant because of its keratolytic effect in the treatment of psoriasis.[51] The “Antineoplastic agent” methotrexate is a well-known medication for psoriasis but several additional drugs in the same class were found in our analysis. Studies are available about 5-fluorouracil for the treatment of dystrophic psoriatic fingernails, but it showed only non-significant improvement.[52] Micellar paclitaxel significantly improved psoriasis in a prospective phase II study.[53] A study reported significant effectiveness of topical caffeine.[54] “Calcium channel blocker” nifedipine is found to be inductor of the disease in a case control study.[55] A study in 2005 reported significant PASI score reduction of 49.9% by topical theophylline ointment.[56] Mahonia aquafolium extract - consisting of berberine among others - is not classified into ATC classes, but three clinical trials already indicated improvement of psoriasis with this substance.[57] Multiple studies prove efficacy of the terpenoid triptolide in the treatment of psoriasis.[58] A recent study investigated effect of rifampicin on psoriasis and reported a 50.03% mean PASI reduction.[59] Study about the treatment of psoriasis with curcumin was carried out but reported only low response rate.[60]

In an in vitro experiment the “Lipid modifying agent” clofibrate, but not bezafibrate reversed UVB-light-mediated expression of psoriasis – related inflammatory cytokines (interleukin-6, interleukin-8).[61] Fluvastatin and pravastatin have the potential to inhibit Th17 cell chemotaxis thus lowering immune cell infiltration of psoriatic skin.[62] Anti-proliferative effect of novel COX2 inhibitors on HaCaT keratinocytes was proven in an in vitro experiment and possible therapeutic use in psoriasis was supposed. However no such experiment was carried out with celecoxib which was the only COX2 inhibitor in best rated drugs.[63] N-acetyl-cysteine attenuated TNF alpha – induced cytokine production in primary human keratinocytes, which suggests its anti-psoriatic potential.[64] The “Thiazolidinedione” ciglitazone was never used as a medication, but inhibited keratinocyte proliferation in a dose dependent fashion.[48] Histone – deacetylase inhibitor trichostatin A blocked the conversion of regulatory T cells to IL17 expressing T cells suggesting its beneficial role in treating psoriasis.[65] Tse et al. suppose that antiproliferative effect of arsenic compounds could have positive effects on psoriatic skin.[66] The phosphodiesterase inhibitor rolipram has the ability to block enterotoxin B-mediated induction of skin homing receptor on T lymphocytes and may have the potential to inhibit lymphocytic infiltration of lesional skin.[67] The natural polyphenolic compound rottlerin is a potent inhibitor of NFκB and may have disease modulating effects.[68]

Case reports are available about psoriasis induction by clonidine, “agents acting on the renin-angiotensin system” like captopril or losartan; the “protein kinase inhibitor” and “antineoplastic agent” imatinib; diclofenac, olanzapine, fluoxetine and chloroquine. Also case reports are available about the beneficial effects of ritonavir; “antineoplastic agents” like cytarabine, doxorubicin, and cysplatin; gefitinib, colchicine, lidocaine and nicotine.[69]–[83]

The 32 effective drugs of “Studies available” group in Table 7 were filtered out from STITCH data and target proteins were analyzed. All target proteins got an in-degree value reflecting the number of effective drugs acting on it. The group of proteins forming the controller sub-network was compared with the remaining target proteins. The controller sub-network protein group got significantly higher median value (10 vs. 1) using Mann-Whitney Rank Sum Test than the other one, which prove the importance of the controller sub-network in psoriatic lesions. (Figure 5) (p<0.001; in-degree has power law distribution, thus T-test could not be used) Higher median value could be caused by higher original degree centralities of controller network proteins in PPI networks, but only weak relation have been found between original degree centrality and the number of effective drugs acting on a protein, which cannot explain the big difference between the median of two groups (corrected R square value in regression analysis: 0.304)

Figure 5

Effect of anti-psoriatic drugs on controller network.

Higher number of effective anti-psoriatic drugs act on controller nodes than on other proteins. Totally the targets of 32 effective anti-psoriatic drugs were analyzed (median 10 vs. 1) *p<0.001

In summary, studies are available for 34 drugs found by our analysis, experimental evidence is available for 24 drugs, case reports suggest beneficial or disease-inductor effect of 21 drugs and 98 unpublished drug candidates for the treatment of psoriasis were also found (Table 7–8).

Table 8

Drug candidates unassociated with psoriasis.

ATC Class	Drug
Retinoids for topical use in acne	retinol
Blood glucose lowering drugs excl. insulines	glyburide
Vitamin K and other hemostatics	menadione
Antineoplastic agents	aldophosphamide, MLS003389283, etoposide, dasatinib, decitabine
Sex hormones and modulators of the genital system	(4–14c)pregn-4-ene-3,20-dione, mifepristone, testosterone-propionate, androstanolone, diethylstilbestrol, raloxifene
Hormone antagonists and related agents	flutamide, fulvestrant
Cardiac stimulants excl. cardiac glycosides	bucladesine
Cardiac glycosides	G-Strophantin
Drugs for obstructive airway diseases	salbutamol
Antiadrenergic agents, centrally acting	reserpine
Antiadrenergic agents, peripherally acting	prazosin
Lipid modifying agents, plain	lovastatin, pitavastatin, fenofibrate
Calcium channel blockers	verapamil
Diuretics	furosemide, spironolactone
Liver therapy	silibinin
Platelet aggregation inhibitors excl. heparin	dipyridamole, cilostazol, amiloride-hydrochloride
Agents acting on the renin-angiotensin system	telmisartan, valsartan
Anaesthetics	ketamine, propofol, cocaine, isoflurane
Analgesics	morphine
Psycholeptics	haloperidol, clozapine, diazepam
Psychoanaleptics	desipramine, amitriptyline, metamphetamine
Antiepileptics	phenobarbital, valproic acid
Antidotes	naloxone
Other nervous system drugs	carbacholin
N/A	cytochalasin D, aminoguanidine, Neurogard, paraquat, Y27632, oxidopamine, nitroarginine, AC1LA4H9, SL327, emodin, 2,3,7,8-tetrachlorodibenzo-dioxin, 3-(2-aminoethyl)-5-[(4-ethoxyphenyl)methylidene]-1,3-thiazolidine-2,4-dione, CHEMBL248238, geldanamycin, anisomycin, 8-bromocyclic GMP, tempol, MK-801, 1-(5-isoquinolinesulfonyl)-2-methylpiperazine, ionomycin, herbimycin, pyrrolidine dithiocarbamate, nordihydroguaiaretic acid, gamma-imino-ATP, forskolin, GMP-Pnp, roscovitine, flavopiridol, N-formyl-Met-Leu-Phe, ns-398, sodium butyrate, AC1L1I8V, tyrphostin B42, kainic acid, pirinixic acid, IBMX, bisindolmaleimide I, proline-dithiocarbamate, KBio2_002303, Zillal, thapsigargin, calcimycin, clenbuterol, indole-3-carbinol, 1,9-pyrazoloanthrone, herbimycin, kaempferol, daidzein, lithium-chloride, naringenin

Discussion

Previous meta-analysis of psoriasis microarray studies was carried out by Tian et al. 1120 DEGs were found using 5 studies and 1832 DEGs using 3 studies.[84] We used the same 5 studies, but samples with inadequate quality were excluded from each study using arrayQualityMetrics package. The high number of DEGs (5363) in our study may be surprising, but it can be caused by the lower gene expression fold change cutoff (1.5 and −1.5 instead of 2 and −2). The pre - filtering process of samples can decrease variance and can also increase the number of DEGs. Further analysis of DEGs was carried out with Ingenuity Pathway Analysis (IPA) by Tian et al. IPA uses published references, carry out gene set enrichment analysis and TF detection. We used fundamentally different analysis. We generated PPI networks based on the largest PPI database (STRING) available which not only contain experimentally proven interactions but highly reliable interactions based on prediction algorithms or data mining. PDI network was also generated using not only literally proven interactions but interactions based on high fidelity prediction algorithms. Using lower DEG fold-change cutoff and detailed analysis based on node centrality statistics made it possible to identify proteins yet not associated with the disease but may have remarkable impact on pathogenesis. A chemical – protein interaction network based on STITCH database was also created and disease – modifying drug prediction was also possible with this method.

Keratinocyte hyperproliferation and Psoriasis

Keratinocyte hyperproliferation and inhibition of apoptosis are well-known phenomena in psoriasis. Several proteins have been associated with these mechanisms like BCL2, BAX, NFATC1, PPARδ, EGF, mTOR, NF-κB etc.[85]–[88] Most of them were in central proteins detected by DEG-derived network analysis. Candidate DEG-coded proteins for hyperproliferation like CCNA2, TFDP1 and MECOM were also found. CCNA2 encodes Cyclin A2, that controls S phase and G2/M transition. Not only cell cycle progression is abnormal in lesional skin, but actin cytoskeleton organization as well.[89] A recent study reported that CCNA2 protein has role in cytoskeletal rearrangements and cell migration as well.[31] Cyclin A2 may take part in hyperproliferation and in aberrant actin cytoskeleton organization in psoriatic skin keratinocytes. TFDP1 encodes DP1 protein which is a dimerization partner of E2F transcription factor. The E2F/DP1 heterodimers regulate cell cycle via DNA replication control and apoptosis. DP1 has E2F-independent function as well: DP1 can stabilize Wnt-on and Wnt-off states in Wnt/β-catenin signaling and determine differential cell fates.[37] TFDP1-regulated genes belong to cell cycle progression as shown by BinGO analysis (Table S5). TFDP1 also has a reciprocal gene expression regulation with AR. This interaction was responsible for motif no. 332 enrichment in psoriasis PDI network compared to complete PDI network. This interaction may connect the hyperproliferation machinery to the merged controller sub-network.

Immunological-metabolic interplay in psoriasis

Psoriasis is an immune-mediated disease. Some proteins which are published as important factors in pathogenesis were absent from DEGs in our microarray-meta analysis, such as TNF alpha, which is an important target in psoriasis therapy. This could be explained by the fact, that increased TNF alpha in psoriatic plaques can be caused mainly by post-transcriptional mechanisms.[90]

Many proteins published in association with the immunopathogenesis of psoriasis were highly ranked hubs in PPI networks: IL1, IL8, TGFB1, SP1, STAT1, STAT3, NFKB1, IRF1 etc.[87], [91]–[97] A highly interconnected cluster mainly consisting of upregulated chemokines and chemokine receptors was also found by PPI analysis (Figure S3). The downregulation of src kinase FYN seems to be a counteracting compensatory mechanism as this protein is important in IFN gamma action, in TNF alpha induced COX2 expression and in adipose tissue - mediated inflammation leading to insulin resistance. These processes are important in the pathomechanism of psoriasis.[32], [98], [99] These data suggest that the FYN inhibitor KBio2_002303 may have beneficial effects in the treatment of psoriasis. An important node in controller sub-network is IL8. Although its role in psoriasis pathogenesis is published, no trial has been done with IL8 inhibitors.[100] This is true for CCL2 and IRF1 as well. Our study confirms their basic role in sustainment of lesional phenotype. Both can be found in highly ranked hubs and CCL2 is also essential in controller sub – network by activating two positive feedback loops related to inflammation.

Psoriasis and metabolic syndrome comorbidity is a well-known phenomenon. There is a complicated interaction between the two diseases mediated by inflammatory cytokines among others.[101] Numerous DEG-coded proteins associated with both diseases could be found in central proteins like PPARG, INS-IGF2, LEP etc.[102]–[104] Others, like PIK3R1, AR and MEF2A may have role in the development of metabolic syndrome in psoriasis. PI3KR1 is important in the development of insulin resistance, it propagates inflammatory response in obese mice and may be an important link between the obesity-inflammation interplay in psoriasis.[33] AR has important effect on insulin signaling and thus insulin resistance. It is published that AR knockout mice exhibit insulin resistance.[35] To our knowledge AR has not yet been associated with psoriasis. However it was found in 1981, that lower serum testosterone level therefore decreased AR activation can be detected in psoriatic patients.[105] AR and PPARG connect inflammation- and metabolism-related hubs in controller network thus modulation of these proteins can be beneficial in psoriatic patients, which was also proven by our drug target analysis (Figure 5). MEF2A is important for GLUT4 expression on insulin-responsive cells. Expression of MEF2A is downregulated in lesional skin samples which suggests a possible mechanism for insulin resistance in psoriasis.

Many drugs, which are already widely used as treatment for psoriasis could be found in highly ranked nodes of chemical-protein interaction networks such as methotrexate, retinoic acid, corticosteroids, sirolimus and tacrolimus. According to STITCH data all of them act through at least one of the hubs in controller sub-network. Top ranked ATC drug classes target members of controller sub-network as well. Blood glucose-lowering drugs act through PPARG and INS-IGF2 activation, which can be the basis of the positive effects of fibrate and HMG-CoA inhibitors in psoriasis as well.[47] Cardiac stimulants such as adrenergic agents also have high impact on lesional skin's PPI and PDI network, mainly by modulating hubs in controller sub-network. “Sex hormones and modulators of the genital system” ATC drug class act on AR. The “antineoplastic drug” methotrexate mainly acts through the accumulation of adenosine, but other antineoplastic agents may have their effect on keratinocyte hyperproliferation.[106] Studies or case reports already suggest efficacy of some antineoplastic drugs but several new possible agents were found in our analysis.[53], [107], [108] Mental stress is a known trigger for psoriasis and connection between the neuroendocrine system and skin immune system has been reported. [3], [109] This is not surprising that numerous drugs acting on the CNS are enriched in highly ranked drugs. A lot of other drugs which are either classified in ATC classes or just drug candidates are found like kainic acid, cocaine, the HDAC inhibitor sodium butyrate, the PKC inhibitor bisindolylmaleimide I etc. (Table 7)

In summary this is the first time PPI, PDI and chemical-protein interaction networks of psoriatic skin samples has been examined with detailed network analysis. Network-building DEGs were identified with fine-quality microarray meta-analysis of 187 non-lesional and 189 lesional samples. Several proteins were found which are yet not associated with psoriasis but may have high impact on the pathogenesis of the disease. Basic disease controller sub-network was also constructed consisting of central nodes coded by DEGs. Numerous anti-psoriatic drugs and drug candidates were also found acting mainly on these nodes.

FYN protein in the jActiveModules cluster with 2 Nodes with blue-shaded color are downregulated and nodes with red-shaded color are upregulated. Color intensity is proportional with fold change.

(PNG)

Click here for additional data file.

PIK3R1 protein in the jActiveModules cluster with 3 Nodes with blue-shaded color are downregulated and nodes with red-shaded color are upregulated. Color intensity is proportional with fold change.

(PNG)

Click here for additional data file.

Chemokine-chemokine receptor cluster found by ClusterONE.

(PNG)

Click here for additional data file.

Outlier samples in arrayQualityMetrics analysis and explanation of exclusion. Numbers from 1 to 5 indicate the number of method used by the software to assess quality.

(XLSX)

Click here for additional data file.

Results of repeated original statistics. The same statistics, cutoff and filtering was used as it is shown in table after arrayQualityMetrics sample filtering and GCRMA normalization.

(XLSX)

Click here for additional data file.

Study-level p values of T test for differential gene expression and meta-analysis FDR values. DEGs are highlighted with orange color.

(XLSX)

Click here for additional data file.

Node centralities in each network. Different networks can be found on different worksheets. Columns indicate centrality values calculated by NetworkAnalyzer. Central nodes are highlighted with orange color.

(XLSX)

Click here for additional data file.

Results of BinGO analysis. Significant GO terms are highlighted with orange color

(XLSX)

Click here for additional data file.

Results of arrayQualityMetrics analysis. Only html data can be found in directories, pdf files were deleted due to size restrictions.

(ZIP)

Click here for additional data file.