the effects of the phosphorylation of serine 273 on the protein structure of ppar were studied by comparing the structures of phosphorylated serine 273 and the wild type , with a particular focus on the dbd of the receptor complex , using md simulations .
although a detailed study on dbd dynamical properties under perturbing stimuli and its relation to dr1-binding specificity is not the subject of this paper , a simple comparison of time - averaged molecular motions can indicate regions of altered stability .
indeed , serine 273 phosphorylation affects the stability of rxr , mostly in the zinc - finger helix which is in contact with dna 's major groove ( figure 2a ) .
the general change in the structure of dbd - rxr is relatively small : 0.5   of averaged atomic displacement ( root - mean s.d . ) ( figure 2b ) .
nevertheless , the whole domain became less stable , as shown in differences in the root mean square fluctuation profiles ( figure 2c ) , with the most remarkable fluctuation occurring in the helical region of one of the zinc fingers .
this indicates that phosphorylation introduces changes in the dna - binding domain , shifting the protein population towards a state which could contribute to an increase in disease - related gene expression .
the lbd of ppar is not affected on a large scale ; despite some differences , both profiles are following a similar trend ( figure 2d ) . to study the effects of modified ppar gene regulation , we reconstructed an interaction network consisting of 235 degs , 152 of which were upregulated in the phosphorylated state ( wild type ) and 48 of which were upregulated in the mutated state ( s273a mutant ) .
subsequently , we searched a sub - network that makes an important contribution to the stability of the global network .
we then selected the top 10 statistically significant motif types consisting of 2 , 3 , 4 or 5 nodes using a z - value - based ranking , all having a p - value < 0.01 ( only 10 of 1000 randomly generated networks contained motifs at a higher frequency ) .
we examined the expression profile of the genes in each motif to ensure that the expression profile was consistent in all members of the identified motif : among the top ranked 2- , 3- , 4- or 5-node motifs that exhibited multi - stable behavior , we identified a number of 3- and 4-node motifs that were consistent with expression values ( figure 3 ) .
we thus identified 39 genes involved in 55 switches forming a single cluster , which we call the  core ' network ( figure 4 ) . computation of the stability of the core network cluster reveals two attractors , in which all genes are either in an on or an off state , which exactly match the expression pattern .
the two most commonly found motifs are positive forward loop motifs with off / off / off changed to the on / on / on states during the course of the disease .
consistent with this , the in silico perturbation of some nodes of the core cluster trigger a transition from the off to the on attractor , but the opposite transition ( from on to off ) is not triggered by the perturbation of any node in the core . in all , 34 nodes in the cluster are capable of triggering the transition from off to on when they are perturbed . in contrast , five nodes ( col1a1 ( collagen type i-1 ) , col1a2 ( collagen type i-2 ) , kruppel - like factor 5 ( klf5 ) , perforin-1 and runt - related transcription factor 1 ) constitute a set of genes that potentially could be in the on state , taking part in different processes without leading to the activation of the cluster as a consequence .
we can see in figure 5 that perturbation of klf5 , involved in a motif with early growth response 1 ( egr1 ) , does not trigger the transition of the cluster from off to on ( see supplementary figure 1 for all perturbation experiments ) .
furthermore , we compared properties of genes in the core network with remaining genes in the global network ( supplementary table 1 ) .
specificity of inter - gene interactions may be reflected in the network modularity . using the newman 
girvan algorithm , we identified 11 clusters of which clusters 1 and 2 are mainly occupied by core genes .
to identify inter - cluster connectors , we calculated the participation coefficient that on average was significantly different when comparing nodes in the core and global networks .
the median participation in the core and global networks was 0.45 and 0 , respectively ; the distributions in the two groups differed significantly ( mann  whitney 
the betweenness centrality of a gene is a centrality measure that is proportional to the number of shortest paths between genes in the network that go through the gene in question .
high betweenness centrality corresponds to high level of inter - node communication , and is therefore an appropriate measure for highlighting which genes link different molecular processes and pathways .
median betweenness centrality in the core and global networks are 243 and 0 , respectively ; the distributions in the two groups differ significantly ( mann  whitney  wilcoxon w=1430 , n1=44 , n2=191 , p - value=3.65e-14 ) .
we identified a group of genes , which act as potential hubs and exist in the core network ( having a betweenness centrality score 0.10 in the core network , as well as 2900 in the global network , being the top four ranking genes in both cases ) .
the genes with the highest betweenness centrality were hif1a ( hypoxia - inducible factor 1-subunit ) , egr1 , stat1 ( signal transducers and activators of transcription 1 ) and cxcl12 ( chemokine ( c  x  c motif ) ligand 12 ) and they are all overexpressed in the case of phosphorylated ppar ( the putative disease state ) .
we examined these genes in more detail with regard to their association with adipogenesis , for which ppar is the master regulator under any condition , to confirm that our core network is consistent with experimental studies .
we found that this was indeed the case : the role of the transcription factor hif1a in adipocyte differentiation has been described previously , egr1 functions as a transcriptional regulator the expression of which is rapidly induced during the differentiation of murine 3t3-l1 adipocytes , stat1 also acts as a transcription activator , which is rapidly activated in the 3t3-l1 adipocyte cell culture model , and cxcl12 is a chemokine , which demonstrates a significant increase in expression in differentiating 3t3-l1 preadipocytes .
the metabolic syndrome is a common disease and is associated with increased incidence of type ii diabetes mellitus and cardiovascular mortality .
its treatment is designed to target the different components of the syndrome ( such as high triglyceride , cholesterol and glucose levels and high blood pressure ) .
thus , the discovery of ppar as a master regulator of adipogenesis , the decreased activity of which was also shown to be linked to insulin resistance , was very important . it was assumed that increasing ppar activity would potentially represent a unifying strategy to target the metabolic syndrome as an entity .
unfortunately , unwanted side effects such as weight gain , heart failure and even osteoporosis , as well as discrepant results between genetic and pharmacological modulations of ppar , outbalance the observed clinical benefits , which mostly rely on improved blood sugar control and lipid metabolism .
ppar still represents a desirable drug target , but research over the last decade showed that its biological function with regard to adipogenesis is a balanced and thrifty response .
aspects such as co - factor recruitment and function and/or post - translational modifications have to be taken into consideration , before more specific / distinct new approaches to modify ppar action can lead to new compounds .
ppar expression is also influenced by gata - binding protein 2 ( gata2 ) , a transcription factor specifically expressed in the stromal vascular fraction of the adipose tissue , which is also where preadipocytes are located .
the transition from preadipocyte to adipocyte is blocked , at least in part , by gata2 's suppression of ppar ; insulin has downstream effects on the kinase v - akt murine thymoma viral oncogene homolog 1 ( akt ) which phosphorylates gata2 , impairing gata2 's translocation to the nucleus , leading to a reduced suppression of ppar in the nucleus , thus promoting adipogenesis . in the insulin - resistant state
, the phosphorylation of gata2 by akt is disturbed and consequently adipogenesis regulation is dysregulated .
ppar-driven adipogenesis is reduced , whereas the remaining preadipocytes start expressing pro - inflammatory cytokines , an observation contributing to the understanding of the puzzling phenomenon of improved metabolic outcome of adipogenesis - promoting ppar agonists .
gata2 was not among the degs in our study , but this is not unexpected as gata2 's activity is regulated through phosphorylation .
metalloproteinase inhibitor 3 ( timp3 ) is a member of our core network and was originally linked to sorsby 's fundus dystrophy , which is an autosomal dominant macular degeneration disorder .
this gene has been linked to a number of other diseases including , more recently , the metabolic syndrome , in particular with insulin resistance , hepatic inflammation , dyslipidemia and atherosclerosis in diabetic subjects . a possible mechanism of how ppar and timp3 are connected may be that timp3 expression influences c / ebp ( ccaat / enhancer - binding protein beta ) expression , which in turn induces c / ebp ( ccaat / enhancer - binding protein alpha ) and ppar during adipocyte differentiation .
post - translational modification can be seen as perturbations of the environment in which a protein exists .
the effect of such a post - translational modification on a protein is a good situation to examine the hypothesis of the population shift of molecules , which proposes that the external stimuli stabilize and therefore populate pre - existing , although previously non - overrepresented molecule states . in biological terms
, a population of protein conformers may be shifted towards malfunctioning , diseased - related conformation(s ) , that is , the sustained perturbation of protein structure is a potential trigger for changes of molecular networks associated with specific disease phenotypes .
choi et al . studied the role of ppar in the context of post - translational modification , that is , the phosphorylation of serine 273 in the lbd of ppar ; a number of genes were identified the differential expression of which is caused by a change in the conformation of the pparrxr dimer when serine 273 is phosphorylated .
this modification was shown to influence gene expression of disease ( obesity , insulin resistance)-relevant genes independently from the general receptor transcriptional activity .
to study the effects of perturbation on a structural level , we performed a md simulation to examine the impact of phosphorylation on the pparrxr dimer , demonstrating that phosphorylated ppar destabilizes the pparrxr dimer , which may be a mechanism explaining the altered transcriptional regulation of ppar target genes .
interestingly , phosphorylation at another residue ( ser112 ) modulates ppar activity with positive effects on glucose and lipid metabolism in vivo possibly by decreasing its affinity for endogenous ligands .
furthermore , we analyzed perturbations of biological networks associated with changes in ppar , by examining the degs comparing wild - type phosphorylated ppar with mutant non - phosphorylated ppar. a network was constructed connecting these degs using the biomedical literature - based resnet mammalian database .
therefore , it can be expected that there will be little bias introduced by orthologs from exotic species , that is , species distant to mouse in the phylogenetic tree .
we found an enrichment of bi - stable switches in wild - type phosphorylated state containing genes that significantly change their expression levels in an adipogenesis - dependent manner .
of particular interest were the genes with the highest betweenness centrality : hif1a , egr1 , stat1 and cxcl12 .
these genes are all overexpressed in the putative disease state and all have been previously described to have a role in adipogenesis .
betweenness centrality is a measure of the influence that a node has over the spread of the information through the network , and our finding that the hub genes with the highest levels of betweenness centrality consist of three transcriptional regulators and one chemokine is consistent with them functioning as information propagators in the network .
the bi - stable motifs we identified assemble an interconnected core network of genes , which have two steady - stable states corresponding to gene functionality states .
this stable core network gives robustness to the cellular network by keeping the system in the non - phosphorylated putative healthy state ; however , once the system is pulled out of this state , it stabilizes the network in a phosphorylated putative disease - related state .
the therapeutic ligand present in the intact structure of the pparrxrdna complex ( pdb i d : 3dzy ) was replaced with docosahexanoic acid , which is a natural ligand of ppar. md simulations were performed with gromacs and served as a tool to generate perturbed states of proteins .
implicit solvent generalized born simulations with the amber99sb - ildn force field were carried out for wild - type and mutant proteins : energy minimizations were performed and the production runs were 50  ns in length . for our analysis , we extracted degs from the results of gene expression analysis experiments by choi et al .
, in which ppar-null mouse embryonic fibroblasts were transfected with wild - type ppar ( phosphorylated ) or the s273a ppar mutant ( not phosphorylated ) .
the cutoffs for selecting the degs were uncorrected p - values 0.008 and corrected false discovery rate p - values 0.15 .
subsequently , the resnet mammalian database from ariadne genomics ( http://www.ariadnegenomics.com/ ) was used to construct an interaction network of 235 degs with directed interactions ( the  global ' network ) .
this database includes biological relationships and associations , which have been extracted from the biomedical literature using ariadne 's medscan technology .
medscan processes sentences from pubmed abstracts and produces a set of regularized logical structures representing the meaning of each sentence .
the resnet mammalian database stores information harvested from the entire pubmed , including more than 715  000 relations for 106  139 proteins , 1220 small molecules , 2175 cellular processes and 3930 diseases .
motif detection was performed using the fanmod algorithm for the global network and limited to motifs consisting of 2 , 3 , 4 or 5 nodes each .
each resultant topology was analyzed in comparison with 1000 separately randomized versions of the initial network . using both the original and the randomized versions , z - scores and p - values
the z - score of a motif is the original frequency minus the random frequency divided by s.d .
the p - value of a motif is the number of random networks in which it occurred more often than in the original network , divided by the total number of random networks . of the bi - stable switch motifs we identified
, we also examined the expression profile of the genes in each motif to ensure that the expression profile was consistent in all motif members .
we constructed a network ( the  core ' network ) consisting of 39 degs that occur in bi - stable motifs , which are significantly overrepresented in the global network and the genes of which are in the stable state , matching the experimental expression values . to compute the attractors of the core network
the program converts the network into a continuous dynamical system based on ordinary differential equations . in the absence of detailed kinetic parameters ,
the program interpolates a sigmoid curve between the states completely on and completely off for each node .
squad first calculates the steady states found in a discrete dynamical system ( boolean model ) and then uses these states as a guide to localize the steady states in the continuous model .
each perturbation is a single pulse that changes the state of the node from 0 to 1 in the off attractor or from 1 to 0 in the on attractor as initial states of the system .
the degs of the core and global networks were divided into groups by topological modularization using the newman  girvan algorithm .
participation ( p ) of a node in intra - modular communication is calculated as follows :   where nm = number of modules , ks = number of connections with module s , ki = total degree of module k. analysis of betweenness centrality was calculated with igraph library in r comparing the core and global networks .
the participation coefficient and betweenness centrality results for the core and global network groups were compared with the mann 
we also examined these genes in more detail with regard to their association with adipogenesis , for which ppar is the master regulator under any condition , to confirm that our core network is consistent with experimental studies .
the therapeutic ligand present in the intact structure of the pparrxrdna complex ( pdb i d : 3dzy ) was replaced with docosahexanoic acid , which is a natural ligand of ppar. md simulations were performed with gromacs and served as a tool to generate perturbed states of proteins .
implicit solvent generalized born simulations with the amber99sb - ildn force field were carried out for wild - type and mutant proteins : energy minimizations were performed and the production runs were 50  ns in length .
for our analysis , we extracted degs from the results of gene expression analysis experiments by choi et al .
, in which ppar-null mouse embryonic fibroblasts were transfected with wild - type ppar ( phosphorylated ) or the s273a ppar mutant ( not phosphorylated ) .
the cutoffs for selecting the degs were uncorrected p - values 0.008 and corrected false discovery rate p - values 0.15 .
subsequently , the resnet mammalian database from ariadne genomics ( http://www.ariadnegenomics.com/ ) was used to construct an interaction network of 235 degs with directed interactions ( the  global ' network ) .
this database includes biological relationships and associations , which have been extracted from the biomedical literature using ariadne 's medscan technology .
medscan processes sentences from pubmed abstracts and produces a set of regularized logical structures representing the meaning of each sentence . the resnet mammalian database stores information harvested from the entire pubmed , including more than 715  000 relations for 106  139 proteins , 1220 small molecules , 2175 cellular processes and 3930 diseases .
motif detection was performed using the fanmod algorithm for the global network and limited to motifs consisting of 2 , 3 , 4 or 5 nodes each .
each resultant topology was analyzed in comparison with 1000 separately randomized versions of the initial network . using both the original and the randomized versions , z - scores and p - values
the z - score of a motif is the original frequency minus the random frequency divided by s.d .
the p - value of a motif is the number of random networks in which it occurred more often than in the original network , divided by the total number of random networks . of the bi - stable switch motifs we identified
, we also examined the expression profile of the genes in each motif to ensure that the expression profile was consistent in all motif members .
we constructed a network ( the  core ' network ) consisting of 39 degs that occur in bi - stable motifs , which are significantly overrepresented in the global network and the genes of which are in the stable state , matching the experimental expression values . to compute the attractors of the core network
the program converts the network into a continuous dynamical system based on ordinary differential equations . in the absence of detailed kinetic parameters ,
the program interpolates a sigmoid curve between the states completely on and completely off for each node .
squad first calculates the steady states found in a discrete dynamical system ( boolean model ) and then uses these states as a guide to localize the steady states in the continuous model .
each perturbation is a single pulse that changes the state of the node from 0 to 1 in the off attractor or from 1 to 0 in the on attractor as initial states of the system .
the degs of the core and global networks were divided into groups by topological modularization using the newman  girvan algorithm .
participation ( p ) of a node in intra - modular communication is calculated as follows :   where nm = number of modules , ks = number of connections with module s , ki = total degree of module k. analysis of betweenness centrality was calculated with igraph library in r comparing the core and global networks .
the participation coefficient and betweenness centrality results for the core and global network groups were compared with the mann 
we also examined these genes in more detail with regard to their association with adipogenesis , for which ppar is the master regulator under any condition , to confirm that our core network is consistent with experimental studies .
 peroxisome proliferator - activated receptor gamma ( ppar ) is a key regulator of adipocyte differentiation and has an important role in metabolic syndrome . 
 phosphorylation of the receptor 's ligand - binding domain at serine 273 has been shown to change the expression of a large number of genes implicated in obesity . 
 the difference in gene expression seen when comparing wild - type phosphorylated with mutant non - phosphorylated ppar may have important consequences for the cellular molecular network , the state of which can be shifted from the healthy to a stable diseased state . 
 we found that a group of differentially expressed genes are involved in bi - stable switches and form a core network , the state of which changes with disease progression . 
 these findings support the idea that bi - stable switches may be a mechanism for locking the core gene network into a diseased state and for efficiently propagating perturbations to more distant regions of the network . a structural analysis of the pparrxr dimer complex supports the hypothesis of a major structural change between the two states , and 
 this may represent an important mechanism leading to the differential expression observed in the core network . 
Results
Discussion
Materials and Methods
Molecular modeling
Network analysis
Supplementary Material