BACKGROUND
the globe figure of people with diabetics is increasing rapidly  <cit> .
the diabetes epidemic worldwide is due to an interaction between environment and genetic risk factors  <cit> .
the modern environment causes diabetes in many ways, such as stress, increased availability of unhealthy food, and decreased physical activities  <cit> .
our body system is a robustness system to keep our blood glucose within normal ranges with various perturbations.
however, in genetically susceptible individuals, long term unfavorable environmental factors will affect epigenetics, thereafter gene expressions, and eventually lead to diabetes.
t2dm is chronic with nature history lasting for more than twenty years, which has been divided into five stages: latent stage, transition stage, impaired glucose tolerance stage , impaired fasting glucose stage , and overt stage  <cit> .
igt and ift stages are called prediabetes.
during the first  <dig> stages, the sub-health status is still able to return to normals.
once reached stage  <dig>  overt stage, t2dm is diagnosed.
the systems of diabetes are also robust: even with food restriction, increased physical activity, and multidrug therapies, diseases are usually impossible to return back to normals  <cit> .
in order to detailed study diabetes, several animal models have been developed.
goto-kakizaki  rat, a spontaneous non insulin dependent diabetes model with a heterogeneous background, is recognized as one of the best model for human t2dm.
the colony was first produced in japan by selective repeated inbreeding nondiabetic wistar-kyoto  rats with minor glucose intolerance  <cit> .
the diabetic state became spontaneous and stable after  <dig> generations.
the characteristics of gk subcolonies are slightly different.
however the important hallmarks are the same, including inherent decreased beta cell mass, moderate hyperglycemia, insulin resistance, and a non-obese phenotype  <cit> .
at embryonic day  <dig>  beta cell mass of gk rats is only 50% of that in normal wky controls.
gk fetuses show decreased insulin levels and decreased beta cell mass.
before  <dig> weeks of age, gk babies show normal blood glucose, but decreased insulin levels.
basal hyperglycemia has been detected at 3- <dig> weeks.
gk rats show unstable blood glucose levels between 6- <dig> weeks and hyperglycemia became consistent in gk rats older than  <dig> weeks of age.
although it exhibits similar metabolic disorders to the human diabetes, gk is non obese without hyperlipidemia at the beginning.
thus it only represents a subset of human t2dm.
t2dm is a systemic metabolic disease.
the two major characters are insulin resistance and beta cells fail to compensate.
liver plays a key role in not only energy metabolism but also insulin resistance, thus liver gene expression changes play a role in the progression of diabetes.
now most scientists agree that the risk of developing t2dm is low with only single gene mutation  <cit> .
environmental factors act on predisposing individuals, changing their dna modification and mrna expression to certain levels until the system is not able to return to normals.
microarray technology makes it easy and accurate to measure significantly changed gene expressions  <cit> .
however, to understand the real meaningful hints from the information ocean and to elucidate the connections between changed biological molecules and diseases seem quite challenge.
it has been recognized that a complex disease cannot be fully understood by merely analyzing individual genes or biomolecules.
it is interactions or networks of those components that are ultimately responsible for malfunctions of the system.
therefore, instead of picking up single interesting gene, we are using network screening to analyze the active networks or pathways based on the high throughput data, a promising approach to investigate associations between biological molecules and phenotypes.
a knowledge-based network is constructed first by extracting as many relationships identified by experimental studies as possible and then superimposing them to microarray data.
recently, we proposes a method  <cit>  to estimate the consistency of a given network with the measured data: i) the network is quantified into a log-likelihood from the measured data, based on the gaussian network, and ii) the probability of the likelihood corresponding to the measured data, named the graph consistency probability , is estimated based on the generalized extreme value distribution.
in this paper, we survey the active regulatory networks in gk and wky rats liver in a comprehensive manner by network screening.
the microarray data measured previously for five liver samples of both groups at each of  <dig> time points  <cit>  are analyzed by the standard statistical techniques and the network screening.
the analyses reveal the expression signatures different between gk and wky rats and the network signatures that are composed of the networks well consistent between the network structure and the graph structure.
as a result, we present the candidates of active regulatory networks, which including new and reasonable networks, as well as the networks previously reported as to be essential to diabetes.
furthermore, we discuss merits and pitfalls of the present approach for surveying the active regulatory networks for a special disease.
materials and methods
• network screening
overview
the candidates of active regulatory networks are detected by network screening in the following manner.
first, the regulatory network sets are generated by combining the binary relationships between transcriptional factors  and their regulating genes, which are compiled in transfac database  <cit> , and the functional gene sets defined in the molecular signatures database   <cit> .
then, we calculate the graph consistency probability   <cit> , which expresses the consistency of a given network structure with the monitored expression data of the constituent genes in this study, for each of the network structures obtained at the first step.
in addition, in each reference network, the enrichment probability of the genes with the significant differences between gk and wky rats  is further tested.
for this purpose, the expression signature is determined using the student’s t-test .
the number of genes included in the expression signature is tested for each network, based on the hyper-geometric probability.
thus, we refine the candidates of active regulatory networks, in terms of both the network structure by gcp and the extent of gene expression by enrichment probability.
the significance of both thresholds is set to  <dig> .
the details of the reference network and the gcp are described, below.
reference network set construction
in the present study, the gcp is estimated for the ensemble of reference networks, to extract the candidate activated networks in gk and wky rats.
the reference networks are constructed using the binary relationships between transcriptional factors and their regulating genes and the classification scheme for gene function.
as for the reference networks, the orthologous genes in rat corresponding all genes in the human binary relationships from transfac database  <cit>  are first investigated, and then the binary relationships in rat that are composed of the orthologous genes to human are constructed.
based on the binary relationships, transcriptional networks are constructed, according to the functional gene sets previously defined in the molecular signatures database   <cit> .
in each gene set, the regulated genes in the binary relationships are searched, and if at least one gene is found in the gene set, then the corresponding binary relationships are regarded as a regulatory network characterized by the gene set.
the set of constructed networks is used as the reference network for network screening.
in present study, the reference network is composed of  <dig>  regulatory networks that are constructed from  <dig>  transcriptional factor-regulated gene relationships.
graph consistency probability
network analysis is based on the procedure for estimating the consistency of a network structure  with the measured data for the constituent variables in the graph  <cit> .
first, the joint density function for a given network  is recursively factorized into conditional density functions according to the parent-descent relationship in the graph  <cit> .
suppose a causal graph is a directed acyclic graph , g, where vi is a vertex  and ej is an edge  in the graph.
the dag can be factorized into subgraphs according to the parent-descent relationships  <cit> .
then, the joint density function f, corresponding to vi for the graph g, can be factorized into the conditional density functions according to the graph, as follows:
where pa{xi} is the set of variables corresponding to the parents of vi in the graph.
second, the causal graph meets the measured data based on the gaussian graphical model   <cit> .
on the assumption that the probability variable xi is subjected to a multiple normal distribution, each conditional function in equation  is obtained by linear regression for the measured data of the constituent nodes  measured at m points, i.e.,
where xik is the measured value of xi, at the k-th point, and ni is the number of variables corresponding to the parents of vi.
thus, the joint density function in equation  is expressed by the regression for the measured data in equation .
finally, the logarithm of the likelihood of the equation  is calculated for the measured data as
thus, the gn allows us to quantify a given network into the corresponding numerical value from the measured data, according to the network form.
note that the calculation of likelihood itself requires no assumptions on the relationships between variables.
indeed, the likelihood can be calculated in the case of non-linear regressions, such as spline regression.
finally, the probability of the log-likelihood for the network structure  is estimated by the distribution of log-likelihoods for many networks generated under the condition that the generated networks shared the same numbers of nodes and edges as those of the given network.
in previous paper, we assume that the generated networks follow the extreme value distribution  <cit> .
in this paper, we generate nr networks under the same condition, and the gcp is simply defined as
where nr is total number of generated networks, and ns is the number of networks with larger log-likelihoods than log-likelihood of tested network.
in the present study, nr is set to  <dig> .
the significance gcp of the given network is set at  <dig>  in this analysis.
enrichment probability
the network signature is additionally evaluated by the number of constituent genes included in the expression signature.
the enrichment probability of the genes in the expression signature for each network is estimated based on the hyper-geometric probability.
when the network is composed of k genes, and l genes are detected in the expression signature, then the probability is obtained by
where m and n are total number of genes in the expression signature, and total number of genes in the reference networks, respectively.
• microarray data
microarray dataset is cited from the national center for biotechnology information  gene expression omnibus  database .
the data are composed of  <dig>  probes measured by using affymetrix microarray suite  <dig>  , which are reduced into  <dig>  genes, for  <dig> samples of male goto-kakizaki  spontaneously diabetic rats and wky rats at each of  <dig> time points .
hyperglycemia begins to show at  <dig> weeks of age and stabilize after  <dig> weeks in gk, thus we divided data into three functional groups: 4w, 8-12w, and 16-20w.
RESULTS
• activated pathways revealed by network screening and their functions
we estimate active regulatory networks among the reference regulatory network set that is generated by the combination of the binary regulatory relationships in transfac database and the functional gene sets defined in the molecular signatures database .
in addition, in each reference network, the enrichment probability of the genes with the significant differences between gk and wky rats is further tested.
finally, we identify a total of  <dig> and  <dig> differentially activating transcriptional regulatory networks in gk and wky rats, respectively.
table  <dig> presents detailed significant networks information separated by ages and strains.
there are fewer pathways activating at 4w and 16-20w in gk rats which are at the beginning and the steady state of diabetes.
while during 8-12w, more pathways are significantly activated, which indicates a dynamic process involving dysfunctions and compensations in the development of diabetes, as showed outside blood glucose fluctuations.
there are more active pathways in the 4w and 8-12w than those in the 16-20w in wky, which may be due to body growth and development.
it is worth pointing out that many activating pathways in wky are diminished in gk rats at 4w, suggesting that those pathways in the liver important to keep glucose metabolism homeostasis are dysfunction at very early stages of diseases.
gk
wky
nguyen_kerato_up
p21_p53_middle_dn
uvb_nhek1_c2
vegfpathway
vegf_huvec_30min_up
yagi_aml_prog_assoc
glycerophospholipid_metabolism
golub_all_vs_aml_up
hohenkirk_monocyte_dend_up
hsc_lateprogenitors_adult
integrinpathway
integrin_mediated_cell_adhesion_kegg
lindstedt_dend_8h_vs_48h_dn
longevitypathway
mef2dpathway
p35alzheimerspathway
rcc_nl_up
brentani_protein_modification
cell_death
hcc_survival_good_vs_poor_up
hsc_lateprogenitors_shared
icf_up
ni2_lung_dn
park_raralpha_mod
schuringa_stat5a_up
brca_brca1_neg
lei_hoxc8_dn
testis_expressed_genes
tsadac_rkoexp_up
apart from the view of differentially activated networks along the time points, the networks in the gk and wky strains can be classified into  <dig> functional categories in table  <dig>  which are metabolism, immune, transcription, and signal transduction.
note that some activated pathways share their functions.
in that case, they are listed under several functional groups as long as the condition met.
then, we combine the activated networks belonging to the same functional category, if any constituent genes of transcriptional factor  and its regulated gene share each other in the networks.
thus tf-gene expression networks for each functional category are created , where the appearance of sub-networks depending on time points is distinguished by colored nodes and edges.
interestingly, significantly activated networks in gk and wky strains are very different even in the same functional category.
we will describe the details of the activated networks in  <dig> functional categories, below.
gk
wky
atria_up
glycerophospholipid_metabolism
golub_all_vs_aml_up
hohenkirk_monocyte_dend_up
hsc_lateprogenitors_adult
longevitypathway
vegf_huvec_30min_up
yagi_aml_prog_assoc
lindstedt_dend_8h_vs_48h_dn
lei_hoxc8_dn
testis_expressed_genes
atria_up
golub_all_vs_aml_up
hohenkirk_monocyte_dend_up
hsc_lateprogenitors_adult
mef2dpathway
hcc_survival_good_vs_poor_up
hsc_lateprogenitors_shared
schuringa_stat5a_up
nuclear_receptors
cell_death
ni2_lung_dn
park_raralpha_mod
nuclear_receptors
integrin_mediated_cell_adhesion_kegg
mef2dpathway
p35alzheimerspathway
rcc_nl_up
vhl_normal_up
aston_oligodendroglia_myelination_ subset
brca_brca1_neg
lei_hoxc8_dn
testis_expressed_genes
tsadac_rkoexp_up
vegfpathway
uvb_nhek1_c2
alkpathway
brentani_protein_modification
cell_death
ni2_lung_dn
park_raralpha_mod
tgfbpathway
• metabolism
metabolic tf regulatory network in wky rats reveals increased expression of several genes are important to keep metabolic homeostasis, e.g.
bone gamma-carboxyglutamic acid-containing protein , hepatocyte nuclear factor  <dig> alpha  and lipoprotein lipase  .
in addition to its role in bone-building, bglap, also known as osteocalcin, acts as a hormone on metabolic regulation.
bglap stimulates pancreatic beta cells releasing more insulin and increases insulin sensitivity via enhancing adipocytes adiponectin secretion  <cit> .
hnf4a plays a key role in liver development.
mutations in this gene have been associated with maturity-onset non-insulin-dependent diabetes of the young   <cit> .
our analysis indicates that reduced hnf4a expression may also favor t2dm development in gk rats.
lpl is an enzyme that hydrolyzes triglyceride in lipoproteins such as very low-density lipoproteins  and reforms high-density lipoproteins .
lipoprotein lipase deficiency leads to elevated levels of triglycerides in the bloodstream.
increment of lpl activity leads to decreased triglycerides levels, elevated hdl levels, a significant fall in fasting glucose and glycohemoglobin, and delayed complication occurrence  <cit> .
interestingly, like hnf4a, lpl is also suggested to be a diabetes susceptibility gene by human studies  <cit> .
metabolic networks in gk rats are more complicated than those in wky rats .
besides the reduced expression of three genes described in the previous paragraph in diabetic gk rats, some pathways identified by network screening further contribute to metabolism disorders.
several signal transducers and activators of transcription  are found in gk tf regulatory metabolic network.
diabetic gk rats are reported to have higher levels of cytokines  <cit> .
cytokines induce activation of janus kinase -stat pathway leading to expression of various suppressors of cytokine signaling  .
checking original microarray data we found that expression of socs <dig> and stat <dig> but not socs <dig> is decreased in gk rats.
decreased expression of socs <dig> leads to enlarged internal organs, which consists with the description in the original paper that liver weight as a percentage of total body weight is significantly larger in gk  <cit> .
insulin directly stimulates socs <dig> and stat <dig> expression, and the decreased socs <dig> and stat <dig> levels are due to insulin deficiency or resistance.
beta cell mass after birth is only half in gk compared to wky rats.
the higher plasma insulin levels in gk measured via millipore ri-13k rat insulin ria kit may be due to cross reaction with elevated proinsulin.
at later stage, insulin resistance also occurs.
igf- <dig>  has a function similar to insulin, and it can also improve blood sugar profiles in type  <dig> diabetics.
igf- <dig> deficiency mice were very insulin insensitive, while administration of igf- <dig> shows the insulin resistance improvement  <cit> .
igf- <dig> levels are increased at 4w, but significantly decreased, thereafter.
while igf <dig> receptor  is exclusively down-regulated, decreased igf1r signaling pathway may partially explain the insulin resistance after  <dig> weeks of age in gk rats.
we also observed some compensative pathways activation in gk to fight against insulin resistance.
for instance, insulin receptor substrate  <dig>  is up-regulated and socs <dig> is down-regulated at 8-12w.
cytokine-induced socs- <dig> interacts with the phosphorylated insulin receptor and promotes ubiquitination  and degradation of ir-irs complex, thereby preventing insulin signaling pathways  <cit> .
decreased socs- <dig> is correlated to insulin sensitivity.
however, compensations fail to stop development of diabetes.
• immune
many proinflammatory pathways are activating in gk compared to wky rats .
from the tf-regulatory gene expression networks in gk rats, two hubs which play important role in immune damages are displayed.
cytochrome b- <dig>  beta polypeptide  is a gene encoding gp <dig> protein, a phagocyte nadph oxidase.
the protein is also known as p91-phox and nox <dig>  reactive oxygen species  produced by nox <dig> are able to kill phagocytized bacteria.
because of its highly reactive nature, cybb has been considered harmful mediators of inflammation  <cit> .
nf-kb and interferon-gamma further increase cybb expression.
prolonged highly cybb expression enhanced production of reactive oxygen species, which are critical sources mediating neurovascular damage.
significantly overexpressed cybb in gk stain is a critical contributor to the microvascular complications associated with diabetes.
activating transcription factor  <dig>  is a stress-inducible gene and encodes atf <dig> transcription factors.
atf <dig> expression has been reported up-regulated in insulitis and type  <dig> or type  <dig> diabetics.
induction of atf <dig> is mediated by proinflammatory factors, such as nitric oxide and nf-κb.
importantly, the induction of atf <dig> leads cell apoptosis, while signals without atf <dig> up-regulation do not cause cell damage  <cit> .
increased gene expression of atf <dig> in gk rats are related to increased immune response and apoptosis.
besides these two hubs, about  <dig> immune related genes are changed in gk strain.
some are up-regulated, such as high affinity immunoglobulin gamma fc receptor i .
some are down-regulated, such as cell surface associated , which protects the body from infection by binding to pathogens.
in sum, inflammation is significantly increased in diabetic gk rats.
• transcription
pathways analysis reveals that wky transcriptional network is a balanced and well-controlled system.
several pathways  are involved in cell replication, good survival and self renewal.
others, including p21-p53_middle_dn, ubv_nhek1_c <dig>  and tgfbpathway, emphasize anticancer and cell cycle checkpoints regulation .
in gk rats, two out of  <dig> pathways are related to apoptosis .
caspase  <dig> , which has been shown to induce cell apoptosis, is overexpressed.
transforming growth factor alpha , which stimulates neural cell proliferation, is inhibited.
interestingly, diabetes activates several genes involving in neurodegenerative disorders.
alzheimer's disease shares many commons with t2dm, so that some scientists proposed to call alzheimer’s disease "type  <dig> diabetes" or "diabetes of the brain."
calpain small subunit  <dig> , a highly-conserved cysteine protease, which have been implicated in neurodegenerative processes after oxidative stress stimulation, is more active in gk.
casein kinase i isoform alpha , also called ck1α, is associated with phosphorylate tau and amyloid formation  <cit> .
reduction in ck1α expression induces tau phosphorylation inhibition.
the expression of ck1α gene is much higher in gk.
• signal transduction
the key difference in signal transduction category is activation of hypoxia and coagulation related pathways in gk rats .
coagulation factor xiii a chain  is the last zymogen activating in the blood coagulation cascade, which stabilize clots  <cit> .
in gk rats, f13a <dig> gene expression levels are significantly elevated which enhance thrombosis.
macrophages expressing high affinity immunoglobulin gamma fc receptor i  also display coagulation function via binding platelets and initiate thrombosis.
<cit> .
tissue plasminogen activator  breakdowns blood clots.
gk rats present significantly higher plat expression levels, which may explain hemolysis and thrombosis co-existing in diabetics.
dr. auwerx reported in diabetics, plat and plasminogen activator  inhibitor are both activated  <cit> .
the elevated levels of pa-inhibitor activity abolish plat activity inducing a reduced fibrinolytic capacity.
rcc_nl_up and vhl_normal_up are two networks involved in hypoxia.
the von hippel-lindau tumor suppressor -hypoxia-inducible factor  pathways are key players in tumor hypoxia survival.
many genes involved in such pathways include interferon regulation factor  <dig>  gtpase hras , and vhl, are negatively expressed in gk strain.
t2dm shows increased incidence and delayed recovery from hypoxia.
reduced hypoxia network activity potentially plays a pivotal role in this phenomenon.
• dynamic changes of regulatory networks
in order to understand the dynamical changes of regulatory networks in the development of diabetes, we drew the active networks at each time segments .
among the genes in the networks, some can be seen in more than one time segment, which are considered to be more important than others, and are distinguished by the colored nodes and edges according to their appearance in which time segments.
furthermore, the information on the expression degree is also important in comparison with gk and wky, and is indicated in the node form in each network.
at the beginning of hyperglycemia, tf regulatory network in 4w gk displays  <dig> genes differently expressed between gk and wky .
those genes can be divided into  <dig> functional groups: immune, metabolism, proliferation and apotosis.
f13a <dig>  cybb, fcgr1a, hck, ctss are involved in inflammation and their expression levels are exclusively increased in gk at  <dig> weeks of age.
previously we have talked about overexpression of cybb and fcgr1a inducing inflammation.
although f13a <dig> is related to thrombosis, it is also been recognized as an inflammation-related gene.
tyrosine-protein kinase  is an enzyme predominantly expressed in hemopoietic cell types.
overexpression of hck contributes to inflammation by promoting neutrophil migration and degranulation as well as couple the fc receptor to the activation of the respiratory burst  <cit> .
cathepsin s  encodes a lysosomal protease that participates in macrophage activation by the degradation of antigens to peptides for presentation  <cit> .
metabolism group includes higher expression of ugdh, abcb <dig>  and soat <dig> genes in gk.
udp-glucose 6-dehydrogenase  converting udp-glucose to udp-glucuronate is significantly increased in dm.
the enhanced expression of ugdh is due to excess glucose load.
multidrug resistance protein  <dig> is a protein that encoded by abcb <dig> gene, which transports phospholipids from hepatocytes into bile.
overexpression is associated with progressive familial intrahepatic cholestasis type  <dig>  sterol o-acyltransferase  <dig> , also known as acyl-coenzyme a: cholesterol acyltransferase, forms cholesterol esters from cholesterol located in the endoplasmic reticulum.
abcb <dig> and soat <dig> are reported coexpressed in gallbladder tissue and participate in bile metabolism  <cit> .
overexpression of soat <dig> functions to atherosclerosis and accumulates cholesterol in the gallbladder mucosa.
recent studies show that bile metabolism is in close contact with occurrence of t2dm.
disturbed bile metabolism has been reported in animal and human diabetes.
bile acid–binding resin prevents and treats diabetes.
diabetes remission after bariatric surgeries is also suggested to be related to changed bile acid metabolism.
analyzing genes in proliferation and apoptosis groups reveal decreased replication.
the proliferation functional group includes reduced expression of hpse, pbk and pold.
heparanase  plays an important role in metastasis and angiogenesis.
lymphokine-activated killer t-cell-originated protein kinase  encodes a mitotic kinase related to mitogen-activated protein kinase kinase  family.
dna polymerase delta catalytic subunit  is a dna polymerase involves in dna repair synthesis after damage.
the apoptosis group including casp <dig>  overexpressing casp <dig> at 4w gk strain causes increased apoptosis.
the tf regulatory network contributing to initial hyperglycemia at 4w continues to be active in 8-12w diabetic gk strain.
in this middle term diabetes, networks making the chaos directing to the diabetes and networks compensate are both active .
a good example is the increased expression level of cathepsin d .
animal and human data suggest that ctsd selectively degrades macrophage inflammatory proteins and is possibly used by tumor to escape antitumoral immune response.
higher expression of ctsd may be secondary to the increased inflammation in diabetics.
however, ctsd will enhance receptor-mediated insulin degradation in vivo, thus inducing insulin resistance  <cit> .
the insulting and compensating battle slowly progress diabetes to next stage.
at stable hyperglycemia stage, fewer networks are activated compared to middle term stage.
however, the insult factors expressed in this stage make diabetes a robust system and unable to return to normals.
• important networks keep normal or diabetes robustness
hyperglycemia is consistent in 16-20w gk rat.
thus we believe that genes expressed at this stage in wky and gk rats are important to keep a steady normal or disease phase.
the first compelling result is the importance of nuclear receptors to maintain the non-diabetic robustness after analyzing tf-regulatory network in the 16-20w wky .
nuclear receptors directly bind to dna, thereby controlling essential biology functions, such as development, homeostasis, and metabolism.
hnf4a, nr3c <dig>  esr <dig>  ar, pparg, nr1d <dig> all belong to nuclear receptor family.
hnf4a belongs to nuclear receptor subfamily  <dig>  nr3c <dig>  esr <dig> and ar are members of subfamily  <dig>  while pparg and nr1d <dig> are included in subfamily  <dig>  they work in concert to defense the disturbance outside.
disease states such as diabetes may be induced by the opposite activities of these receptors.
hepatocyte nuclear factor  <dig> alpha  has been described previously in metabolism section.
it directly regulates genes involved in glucose transport and glycolysis.
estrogen receptor alpha  and androgen receptor  are activated by the sex hormone estrogen and androgen, respectively.
numerous data suggest that estrogen improves glucose metabolism and plasma lipids in t2dm  <cit> .
ar deficiency plays key roles in the development of insulin and leptin resistance, which explains increased diabetes incidence in elder male  <cit> .
the glucocorticoid receptor, also known as nr3c <dig>  is expressed in almost every cell controlling the development, metabolism, especially immune response.
nr3c <dig> decreases inflammation.
peroxisome proliferator–activated receptor-γ  regulates fatty acid storage and glucose metabolism, thus improve insulin sensitivity without increased insulin secretion.
many insulin sensitizing drugs are pparg agonists  <cit> .
n subfamily  <dig>  group d, member  <dig>  also known as rev-erba activates histone deacetylation, thereby regulating gene expression.
publications indicate that snps in these nuclear receptors associate with obesity and/or diabetes.
our data suggest that decreased expression of hnf4a, nr3c <dig>  esr <dig>  ar, pparg and nr1d <dig> overexpression contribute to t2dm.
in gk rats, some compensational pathways still exist, for example a no synthesis pathway is up-regulated.
three genes nitric oxide synthase  <dig> , argininosuccinate synthetase , and nadh: quinone oxidoreductase  related to this pathway are overexpressed.
it is well-known that no decreases blood pressure and promotes vascular actions of insulin.
nos <dig> catalyzes arginine, oxygen and nadph to no and citrulline.
ass <dig> and nq01contribute to this metabolism cycle.
many cytokines increase no regeneration several folds.
increased no synthesis pathway indicates an inflammation environment in the liver in gk rats.
because reduced cell no action has been reported in diabetes, the beneficial effects of increased no production is uncertain.
data analyze reveal increased insulin resistance, hypertension and apoptosis are important to push diabetes to next stage .
protein kinase c alpha  is mostly expressed in hepatocytes promoting glycogenolysis and gluconeogenesis.
activation of prkca mediates serine/threonine phosphorylation of the insulin receptor resulting in decreased active form of insulin receptor, inducing insulin resistance  <cit> .
angiotensin i converting enzyme  <dig>  is an exopeptidase that catalyses angiotensin peptides and has opposite effects on ras axis.
thus decreased expression levels of ace <dig> accelerate the pathologic process such as hypertension, inflammation, fibrosis and inflammation.
gap junction alpha- <dig> also known as connexin- <dig>  is a component of gap junctions providing a route for cell to cell communication via diffusion materials.
decreased gja <dig> expression particularly in hyperglycemia accelerates apoptosis.
• advantages of network screening over single gene based method
when comparing our results to the original study conducted by dr. almon  <cit> , network screening is clearly superior to the single gene-based analysis.
one good example is to explain how liver insulin resistance  develops.
ir is the major character of t2dm and also present in gk rats after  <dig> weeks of age.
in the original study, authors notice higher expression of p <dig>  thus suspecting interaction of p <dig> with irs leading to ir.
however, we believe that the developing ir is a dynamic process involving many steps.
the first step could be significantly decreased igf-1r expression after  <dig> weeks inducing ir in gk.
after that, higher expression of ctsd accelerates ir.
compensational pathways also occur, which includes irs <dig> overexpression at 8-12w in gk.
however as pkc overexpression plus decreased expression of many nuclear factors such as pparg at 16-20w, ir deteriorates and diabetes becomes un-returnable.
our method is based on the networks and is very different from the gene-based method of identifying the differential expression.
discussion and
CONCLUSIONS
t2dm is a complex disease, which is usually not caused by individual gene changes, thereby requiring systems biology methods to understand their mechanisms.
in this work, we have performed comprehensive active regulatory network survey by network screening to the published gk vs. wky liver microarray data  <cit> .
available resources from msigdb and transfac are combined together to identify the significant pathways responsive to the status of diabetes or normals.
after combining the networks according to features or time points, we built functional or time series tf regulatory network graphs.
analyzing the graphs reveals:  <dig>  more pathways are active during inter-middle stage diabetes;  <dig>  inflammation, hypoxia, increased apoptosis, decreased proliferation, and altered metabolism are characteristics in gk strain, and displayed as early as 4w.
<dig>  diabetes progression accompanies insults and compensations.
<dig>  nuclear receptors work in concert to maintain normal glycemic robustness system.
network-based analysis based on high throughput data is a challenging issue, which is expected to help us understand complex disease such as diabetes and further elucidate the essential mechanisms of living organisms which would escape conventional single gene-based analysis.
in this paper, instead of picking up differently expressed genes from high-throughput data, we use known functional pathways to screen datasets and evaluate significantly activated pathways.
then genes with no annotated linkages to tf are overlooked and the available gene regulatory relationships are integrated to form a comprehensive tf regulatory network, which cannot be achieved by single gene based method.
the network shows a whole picture of activated tf regulated functional gene sets under certain conditions and is much easier to bring the biological insights to us.
to our knowledge, two conclusions have not been reported before.
the first one comes out from tf regulatory network at 4w gk.
it is well-known that the major cause of diabetes in gk rats is insulin secreting beta cell dysfuction.
beta cell mass in gk is only half of that in wky after birth.
to be surprised, we find that at very early age liver already exhibits serious gene expression alternations involving in bile metabolism dysfunction, inflammation, increased apoptosis and decreased proliferation, which greatly contribute to diabetes development.
another interesting finding is that the  <dig> nuclear receptors working in concert to maintain robustness of normal blood glucose.
although the relationships of those nuclear receptors with diabetes have been investigated individually before, it is the first time to report how they work together as a fine tune.
restoring their network regulation may have important therapeutic potentials.
this is the first time to use network screening to explain the role of liver in development of diabetes and the underline mechanism.
the results provide many important rational information and insights into guiding experiments design.
it is worth pointing out that the molecular relationships change dynamically, depending on the conditions in a living cell, which suggests implicitly that all of the relationships in the knowledge-based network do not always exist.
note that some methods are proposed for identifying the active networks from measured data  <cit> .
our method evaluates the networks from only one set of data measured under one condition to estimate the absolute consistency between network structure and the data, while the other methods generally need the two sets of data to estimate their relative difference by some criteria such as mutual information.
we combined various resources together to identify the significant regulatory networks related to the development stages of diabetes.
the matching between networks and gene expression profiling was identified by the evaluation of network screening.
the active regulatory networks are the potential disease signatures from the comparison of gk and wky rats.
the dynamics of regulatory networks indicate the dysfunctional progression from the network perspective.
in conclusion, network screening is a superior approach to analyze complex disease such as diabetes.
the conclusions drawn from this method are more complete and systemic, which gives biologist better guidance for further experiment design.
actually, we are now extending this approach for screening general biomolecular networks  <cit>  with both directed and undirected edges, and in future possibly for studying the problem of networkomics  which covers all stable forms of biomolecular networks  <cit>  not only at different biological conditions but also at different spatiotemporal situations.
abbreviations
t2dm: type  <dig> diabetes mellitus; gk: goto-kakizaki; wky: wistar-kyoto; igt: impaired glucose tolerance stage; ift: impaired fasting glucose stage; gcp: graph consistency probability; tf: transcriptional factor; msigdb: molecular signatures database; fdr: false discovery rate; dag: directed acyclic graph; gn: gaussian network; geo: gene expression omnibus; mody: maturity-onset non-insulin-dependent diabetes of the young
competing interests
the authors declare that they have no competing interests.
authors' contributions
hz, kh and lc conceived the research.
hz, ss and zpl performed the study.
gp and jw gave valuable suggestions and improvements.
lc and hz supervised the project.
hz and zpl drafted a version of the manuscript.
all authors wrote and approved the manuscript.

 the globe figure of people with diabetics is increasing rapidly. the diabetes epidemic worldwide is due to an interaction between environment and genetic risk factors. the modern environment causes diabetes in many ways, such as stress, increased availability of unhealthy food, and decreased physical activities. our body system is a robustness system to keep our blood glucose within normal ranges with various perturbations. however, in genetically susceptible individuals, long term unfavorable environmental factors will affect epigenetics, thereafter gene expressions, and eventually lead to diabetes. t2dm is chronic with nature history lasting for more than twenty years, which has been divided into five stages: latent stage, transition stage, impaired glucose tolerance stage, impaired fasting glucose stage, and overt stage. igt and ift stages are called prediabetes. during the first stages, the sub-health status is still able to return to normals. once reached stage overt stage, t2dm is diagnosed. the systems of diabetes are also robust: even with food restriction, increased physical activity, and multidrug therapies, diseases are usually impossible to return back to normals. in order to detailed study diabetes, several animal models have been developed. goto-kakizaki rat, a spontaneous non insulin dependent diabetes model with a heterogeneous background, is recognized as one of the best model for human t2dm. the colony was first produced in japan by selective repeated inbreeding nondiabetic wistar-kyoto rats with minor glucose intolerance. the diabetic state became spontaneous and stable after generations. the characteristics of gk subcolonies are slightly different. however the important hallmarks are the same, including inherent decreased beta cell mass, moderate hyperglycemia, insulin resistance, and a non-obese phenotype. at embryonic day beta cell mass of gk rats is only 50% of that in normal wky controls. gk fetuses show decreased insulin levels and decreased beta cell mass. before weeks of age, gk babies show normal blood glucose, but decreased insulin levels. basal hyperglycemia has been detected at 3- weeks. gk rats show unstable blood glucose levels between 6- weeks and hyperglycemia became consistent in gk rats older than weeks of age. although it exhibits similar metabolic disorders to the human diabetes, gk is non obese without hyperlipidemia at the beginning. thus it only represents a subset of human t2dm. t2dm is a systemic metabolic disease. the two major characters are insulin resistance and beta cells fail to compensate. liver plays a key role in not only energy metabolism but also insulin resistance, thus liver gene expression changes play a role in the progression of diabetes. now most scientists agree that the risk of developing t2dm is low with only single gene mutation. environmental factors act on predisposing individuals, changing their dna modification and mrna expression to certain levels until the system is not able to return to normals. microarray technology makes it easy and accurate to measure significantly changed gene expressions. however, to understand the real meaningful hints from the information ocean and to elucidate the connections between changed biological molecules and diseases seem quite challenge. it has been recognized that a complex disease cannot be fully understood by merely analyzing individual genes or biomolecules. it is interactions or networks of those components that are ultimately responsible for malfunctions of the system. therefore, instead of picking up single interesting gene, we are using network screening to analyze the active networks or pathways based on the high throughput data, a promising approach to investigate associations between biological molecules and phenotypes. a knowledge-based network is constructed first by extracting as many relationships identified by experimental studies as possible and then superimposing them to microarray data. recently, we proposes a method to estimate the consistency of a given network with the measured data: i) the network is quantified into a log-likelihood from the measured data, based on the gaussian network, and ii) the probability of the likelihood corresponding to the measured data, named the graph consistency probability, is estimated based on the generalized extreme value distribution. in this paper, we survey the active regulatory networks in gk and wky rats liver in a comprehensive manner by network screening. the microarray data measured previously for five liver samples of both groups at each of time points are analyzed by the standard statistical techniques and the network screening. the analyses reveal the expression signatures different between gk and wky rats and the network signatures that are composed of the networks well consistent between the network structure and the graph structure. as a result, we present the candidates of active regulatory networks, which including new and reasonable networks, as well as the networks previously reported as to be essential to diabetes. furthermore, we discuss merits and pitfalls of the present approach for surveying the active regulatory networks for a special disease. materials and methods • network screening overview the candidates of active regulatory networks are detected by network screening in the following manner. first, the regulatory network sets are generated by combining the binary relationships between transcriptional factors and their regulating genes, which are compiled in transfac database, and the functional gene sets defined in the molecular signatures database. then, we calculate the graph consistency probability, which expresses the consistency of a given network structure with the monitored expression data of the constituent genes in this study, for each of the network structures obtained at the first step. in addition, in each reference network, the enrichment probability of the genes with the significant differences between gk and wky rats is further tested. for this purpose, the expression signature is determined using the student’s t-test. the number of genes included in the expression signature is tested for each network, based on the hyper-geometric probability. thus, we refine the candidates of active regulatory networks, in terms of both the network structure by gcp and the extent of gene expression by enrichment probability. the significance of both thresholds is set to. the details of the reference network and the gcp are described, below. reference network set construction in the present study, the gcp is estimated for the ensemble of reference networks, to extract the candidate activated networks in gk and wky rats. the reference networks are constructed using the binary relationships between transcriptional factors and their regulating genes and the classification scheme for gene function. as for the reference networks, the orthologous genes in rat corresponding all genes in the human binary relationships from transfac database are first investigated, and then the binary relationships in rat that are composed of the orthologous genes to human are constructed. based on the binary relationships, transcriptional networks are constructed, according to the functional gene sets previously defined in the molecular signatures database. in each gene set, the regulated genes in the binary relationships are searched, and if at least one gene is found in the gene set, then the corresponding binary relationships are regarded as a regulatory network characterized by the gene set. the set of constructed networks is used as the reference network for network screening. in present study, the reference network is composed of regulatory networks that are constructed from transcriptional factor-regulated gene relationships. graph consistency probability network analysis is based on the procedure for estimating the consistency of a network structure with the measured data for the constituent variables in the graph. first, the joint density function for a given network is recursively factorized into conditional density functions according to the parent-descent relationship in the graph. suppose a causal graph is a directed acyclic graph, g, where vi is a vertex and ej is an edge in the graph. the dag can be factorized into subgraphs according to the parent-descent relationships. then, the joint density function f, corresponding to vi for the graph g, can be factorized into the conditional density functions according to the graph, as follows: where pa{xi} is the set of variables corresponding to the parents of vi in the graph. second, the causal graph meets the measured data based on the gaussian graphical model. on the assumption that the probability variable xi is subjected to a multiple normal distribution, each conditional function in equation is obtained by linear regression for the measured data of the constituent nodes measured at m points, i.e., where xik is the measured value of xi, at the k-th point, and ni is the number of variables corresponding to the parents of vi. thus, the joint density function in equation is expressed by the regression for the measured data in equation. finally, the logarithm of the likelihood of the equation is calculated for the measured data as thus, the gn allows us to quantify a given network into the corresponding numerical value from the measured data, according to the network form. note that the calculation of likelihood itself requires no assumptions on the relationships between variables. indeed, the likelihood can be calculated in the case of non-linear regressions, such as spline regression. finally, the probability of the log-likelihood for the network structure is estimated by the distribution of log-likelihoods for many networks generated under the condition that the generated networks shared the same numbers of nodes and edges as those of the given network. in previous paper, we assume that the generated networks follow the extreme value distribution. in this paper, we generate nr networks under the same condition, and the gcp is simply defined as where nr is total number of generated networks, and ns is the number of networks with larger log-likelihoods than log-likelihood of tested network. in the present study, nr is set to. the significance gcp of the given network is set at in this analysis. enrichment probability the network signature is additionally evaluated by the number of constituent genes included in the expression signature. the enrichment probability of the genes in the expression signature for each network is estimated based on the hyper-geometric probability. when the network is composed of k genes, and l genes are detected in the expression signature, then the probability is obtained by where m and n are total number of genes in the expression signature, and total number of genes in the reference networks, respectively. • microarray data microarray dataset is cited from the national center for biotechnology information gene expression omnibus database. the data are composed of probes measured by using affymetrix microarray suite, which are reduced into genes, for samples of male goto-kakizaki spontaneously diabetic rats and wky rats at each of time points. hyperglycemia begins to show at weeks of age and stabilize after weeks in gk, thus we divided data into three functional groups: 4w, 8-12w, and 16-20w.  • activated pathways revealed by network screening and their functions we estimate active regulatory networks among the reference regulatory network set that is generated by the combination of the binary regulatory relationships in transfac database and the functional gene sets defined in the molecular signatures database. in addition, in each reference network, the enrichment probability of the genes with the significant differences between gk and wky rats is further tested. finally, we identify a total of and differentially activating transcriptional regulatory networks in gk and wky rats, respectively. table presents detailed significant networks information separated by ages and strains. there are fewer pathways activating at 4w and 16-20w in gk rats which are at the beginning and the steady state of diabetes. while during 8-12w, more pathways are significantly activated, which indicates a dynamic process involving dysfunctions and compensations in the development of diabetes, as showed outside blood glucose fluctuations. there are more active pathways in the 4w and 8-12w than those in the 16-20w in wky, which may be due to body growth and development. it is worth pointing out that many activating pathways in wky are diminished in gk rats at 4w, suggesting that those pathways in the liver important to keep glucose metabolism homeostasis are dysfunction at very early stages of diseases. gk wky nguyen_kerato_up p21_p53_middle_dn uvb_nhek1_c2 vegfpathway vegf_huvec_30min_up yagi_aml_prog_assoc glycerophospholipid_metabolism golub_all_vs_aml_up hohenkirk_monocyte_dend_up hsc_lateprogenitors_adult integrinpathway integrin_mediated_cell_adhesion_kegg lindstedt_dend_8h_vs_48h_dn longevitypathway mef2dpathway p35alzheimerspathway rcc_nl_up brentani_protein_modification cell_death hcc_survival_good_vs_poor_up hsc_lateprogenitors_shared icf_up ni2_lung_dn park_raralpha_mod schuringa_stat5a_up brca_brca1_neg lei_hoxc8_dn testis_expressed_genes tsadac_rkoexp_up apart from the view of differentially activated networks along the time points, the networks in the gk and wky strains can be classified into functional categories in table which are metabolism, immune, transcription, and signal transduction. note that some activated pathways share their functions. in that case, they are listed under several functional groups as long as the condition met. then, we combine the activated networks belonging to the same functional category, if any constituent genes of transcriptional factor and its regulated gene share each other in the networks. thus tf-gene expression networks for each functional category are created, where the appearance of sub-networks depending on time points is distinguished by colored nodes and edges. interestingly, significantly activated networks in gk and wky strains are very different even in the same functional category. we will describe the details of the activated networks in functional categories, below. gk wky atria_up glycerophospholipid_metabolism golub_all_vs_aml_up hohenkirk_monocyte_dend_up hsc_lateprogenitors_adult longevitypathway vegf_huvec_30min_up yagi_aml_prog_assoc lindstedt_dend_8h_vs_48h_dn lei_hoxc8_dn testis_expressed_genes atria_up golub_all_vs_aml_up hohenkirk_monocyte_dend_up hsc_lateprogenitors_adult mef2dpathway hcc_survival_good_vs_poor_up hsc_lateprogenitors_shared schuringa_stat5a_up nuclear_receptors cell_death ni2_lung_dn park_raralpha_mod nuclear_receptors integrin_mediated_cell_adhesion_kegg mef2dpathway p35alzheimerspathway rcc_nl_up vhl_normal_up aston_oligodendroglia_myelination_ subset brca_brca1_neg lei_hoxc8_dn testis_expressed_genes tsadac_rkoexp_up vegfpathway uvb_nhek1_c2 alkpathway brentani_protein_modification cell_death ni2_lung_dn park_raralpha_mod tgfbpathway • metabolism metabolic tf regulatory network in wky rats reveals increased expression of several genes are important to keep metabolic homeostasis, e.g. bone gamma-carboxyglutamic acid-containing protein, hepatocyte nuclear factor alpha and lipoprotein lipase. in addition to its role in bone-building, bglap, also known as osteocalcin, acts as a hormone on metabolic regulation. bglap stimulates pancreatic beta cells releasing more insulin and increases insulin sensitivity via enhancing adipocytes adiponectin secretion. hnf4a plays a key role in liver development. mutations in this gene have been associated with maturity-onset non-insulin-dependent diabetes of the young. our analysis indicates that reduced hnf4a expression may also favor t2dm development in gk rats. lpl is an enzyme that hydrolyzes triglyceride in lipoproteins such as very low-density lipoproteins and reforms high-density lipoproteins. lipoprotein lipase deficiency leads to elevated levels of triglycerides in the bloodstream. increment of lpl activity leads to decreased triglycerides levels, elevated hdl levels, a significant fall in fasting glucose and glycohemoglobin, and delayed complication occurrence. interestingly, like hnf4a, lpl is also suggested to be a diabetes susceptibility gene by human studies. metabolic networks in gk rats are more complicated than those in wky rats. besides the reduced expression of three genes described in the previous paragraph in diabetic gk rats, some pathways identified by network screening further contribute to metabolism disorders. several signal transducers and activators of transcription are found in gk tf regulatory metabolic network. diabetic gk rats are reported to have higher levels of cytokines. cytokines induce activation of janus kinase -stat pathway leading to expression of various suppressors of cytokine signaling. checking original microarray data we found that expression of socs and stat but not socs is decreased in gk rats. decreased expression of socs leads to enlarged internal organs, which consists with the description in the original paper that liver weight as a percentage of total body weight is significantly larger in gk. insulin directly stimulates socs and stat expression, and the decreased socs and stat levels are due to insulin deficiency or resistance. beta cell mass after birth is only half in gk compared to wky rats. the higher plasma insulin levels in gk measured via millipore ri-13k rat insulin ria kit may be due to cross reaction with elevated proinsulin. at later stage, insulin resistance also occurs. igf- has a function similar to insulin, and it can also improve blood sugar profiles in type diabetics. igf- deficiency mice were very insulin insensitive, while administration of igf- shows the insulin resistance improvement. igf- levels are increased at 4w, but significantly decreased, thereafter. while igf receptor is exclusively down-regulated, decreased igf1r signaling pathway may partially explain the insulin resistance after weeks of age in gk rats. we also observed some compensative pathways activation in gk to fight against insulin resistance. for instance, insulin receptor substrate is up-regulated and socs is down-regulated at 8-12w. cytokine-induced socs- interacts with the phosphorylated insulin receptor and promotes ubiquitination and degradation of ir-irs complex, thereby preventing insulin signaling pathways. decreased socs- is correlated to insulin sensitivity. however, compensations fail to stop development of diabetes. • immune many proinflammatory pathways are activating in gk compared to wky rats. from the tf-regulatory gene expression networks in gk rats, two hubs which play important role in immune damages are displayed. cytochrome b- beta polypeptide is a gene encoding gp protein, a phagocyte nadph oxidase. the protein is also known as p91-phox and nox reactive oxygen species produced by nox are able to kill phagocytized bacteria. because of its highly reactive nature, cybb has been considered harmful mediators of inflammation. nf-kb and interferon-gamma further increase cybb expression. prolonged highly cybb expression enhanced production of reactive oxygen species, which are critical sources mediating neurovascular damage. significantly overexpressed cybb in gk stain is a critical contributor to the microvascular complications associated with diabetes. activating transcription factor is a stress-inducible gene and encodes atf transcription factors. atf expression has been reported up-regulated in insulitis and type or type diabetics. induction of atf is mediated by proinflammatory factors, such as nitric oxide and nf-κb. importantly, the induction of atf leads cell apoptosis, while signals without atf up-regulation do not cause cell damage. increased gene expression of atf in gk rats are related to increased immune response and apoptosis. besides these two hubs, about immune related genes are changed in gk strain. some are up-regulated, such as high affinity immunoglobulin gamma fc receptor i. some are down-regulated, such as cell surface associated, which protects the body from infection by binding to pathogens. in sum, inflammation is significantly increased in diabetic gk rats. • transcription pathways analysis reveals that wky transcriptional network is a balanced and well-controlled system. several pathways are involved in cell replication, good survival and self renewal. others, including p21-p53_middle_dn, ubv_nhek1_c and tgfbpathway, emphasize anticancer and cell cycle checkpoints regulation. in gk rats, two out of pathways are related to apoptosis. caspase, which has been shown to induce cell apoptosis, is overexpressed. transforming growth factor alpha, which stimulates neural cell proliferation, is inhibited. interestingly, diabetes activates several genes involving in neurodegenerative disorders. alzheimer's disease shares many commons with t2dm, so that some scientists proposed to call alzheimer’s disease "type diabetes" or "diabetes of the brain." calpain small subunit, a highly-conserved cysteine protease, which have been implicated in neurodegenerative processes after oxidative stress stimulation, is more active in gk. casein kinase i isoform alpha, also called ck1α, is associated with phosphorylate tau and amyloid formation. reduction in ck1α expression induces tau phosphorylation inhibition. the expression of ck1α gene is much higher in gk. • signal transduction the key difference in signal transduction category is activation of hypoxia and coagulation related pathways in gk rats. coagulation factor xiii a chain is the last zymogen activating in the blood coagulation cascade, which stabilize clots. in gk rats, f13a gene expression levels are significantly elevated which enhance thrombosis. macrophages expressing high affinity immunoglobulin gamma fc receptor i also display coagulation function via binding platelets and initiate thrombosis.. tissue plasminogen activator breakdowns blood clots. gk rats present significantly higher plat expression levels, which may explain hemolysis and thrombosis co-existing in diabetics. dr. auwerx reported in diabetics, plat and plasminogen activator inhibitor are both activated. the elevated levels of pa-inhibitor activity abolish plat activity inducing a reduced fibrinolytic capacity. rcc_nl_up and vhl_normal_up are two networks involved in hypoxia. the von hippel-lindau tumor suppressor -hypoxia-inducible factor pathways are key players in tumor hypoxia survival. many genes involved in such pathways include interferon regulation factor gtpase hras, and vhl, are negatively expressed in gk strain. t2dm shows increased incidence and delayed recovery from hypoxia. reduced hypoxia network activity potentially plays a pivotal role in this phenomenon. • dynamic changes of regulatory networks in order to understand the dynamical changes of regulatory networks in the development of diabetes, we drew the active networks at each time segments. among the genes in the networks, some can be seen in more than one time segment, which are considered to be more important than others, and are distinguished by the colored nodes and edges according to their appearance in which time segments. furthermore, the information on the expression degree is also important in comparison with gk and wky, and is indicated in the node form in each network. at the beginning of hyperglycemia, tf regulatory network in 4w gk displays genes differently expressed between gk and wky. those genes can be divided into functional groups: immune, metabolism, proliferation and apotosis. f13a cybb, fcgr1a, hck, ctss are involved in inflammation and their expression levels are exclusively increased in gk at weeks of age. previously we have talked about overexpression of cybb and fcgr1a inducing inflammation. although f13a is related to thrombosis, it is also been recognized as an inflammation-related gene. tyrosine-protein kinase is an enzyme predominantly expressed in hemopoietic cell types. overexpression of hck contributes to inflammation by promoting neutrophil migration and degranulation as well as couple the fc receptor to the activation of the respiratory burst. cathepsin s encodes a lysosomal protease that participates in macrophage activation by the degradation of antigens to peptides for presentation. metabolism group includes higher expression of ugdh, abcb and soat genes in gk. udp-glucose 6-dehydrogenase converting udp-glucose to udp-glucuronate is significantly increased in dm. the enhanced expression of ugdh is due to excess glucose load. multidrug resistance protein is a protein that encoded by abcb gene, which transports phospholipids from hepatocytes into bile. overexpression is associated with progressive familial intrahepatic cholestasis type sterol o-acyltransferase, also known as acyl-coenzyme a: cholesterol acyltransferase, forms cholesterol esters from cholesterol located in the endoplasmic reticulum. abcb and soat are reported coexpressed in gallbladder tissue and participate in bile metabolism. overexpression of soat functions to atherosclerosis and accumulates cholesterol in the gallbladder mucosa. recent studies show that bile metabolism is in close contact with occurrence of t2dm. disturbed bile metabolism has been reported in animal and human diabetes. bile acid–binding resin prevents and treats diabetes. diabetes remission after bariatric surgeries is also suggested to be related to changed bile acid metabolism. analyzing genes in proliferation and apoptosis groups reveal decreased replication. the proliferation functional group includes reduced expression of hpse, pbk and pold. heparanase plays an important role in metastasis and angiogenesis. lymphokine-activated killer t-cell-originated protein kinase encodes a mitotic kinase related to mitogen-activated protein kinase kinase family. dna polymerase delta catalytic subunit is a dna polymerase involves in dna repair synthesis after damage. the apoptosis group including casp overexpressing casp at 4w gk strain causes increased apoptosis. the tf regulatory network contributing to initial hyperglycemia at 4w continues to be active in 8-12w diabetic gk strain. in this middle term diabetes, networks making the chaos directing to the diabetes and networks compensate are both active. a good example is the increased expression level of cathepsin d. animal and human data suggest that ctsd selectively degrades macrophage inflammatory proteins and is possibly used by tumor to escape antitumoral immune response. higher expression of ctsd may be secondary to the increased inflammation in diabetics. however, ctsd will enhance receptor-mediated insulin degradation in vivo, thus inducing insulin resistance. the insulting and compensating battle slowly progress diabetes to next stage. at stable hyperglycemia stage, fewer networks are activated compared to middle term stage. however, the insult factors expressed in this stage make diabetes a robust system and unable to return to normals. • important networks keep normal or diabetes robustness hyperglycemia is consistent in 16-20w gk rat. thus we believe that genes expressed at this stage in wky and gk rats are important to keep a steady normal or disease phase. the first compelling result is the importance of nuclear receptors to maintain the non-diabetic robustness after analyzing tf-regulatory network in the 16-20w wky. nuclear receptors directly bind to dna, thereby controlling essential biology functions, such as development, homeostasis, and metabolism. hnf4a, nr3c esr ar, pparg, nr1d all belong to nuclear receptor family. hnf4a belongs to nuclear receptor subfamily nr3c esr and ar are members of subfamily while pparg and nr1d are included in subfamily they work in concert to defense the disturbance outside. disease states such as diabetes may be induced by the opposite activities of these receptors. hepatocyte nuclear factor alpha has been described previously in metabolism section. it directly regulates genes involved in glucose transport and glycolysis. estrogen receptor alpha and androgen receptor are activated by the sex hormone estrogen and androgen, respectively. numerous data suggest that estrogen improves glucose metabolism and plasma lipids in t2dm. ar deficiency plays key roles in the development of insulin and leptin resistance, which explains increased diabetes incidence in elder male. the glucocorticoid receptor, also known as nr3c is expressed in almost every cell controlling the development, metabolism, especially immune response. nr3c decreases inflammation. peroxisome proliferator–activated receptor-γ regulates fatty acid storage and glucose metabolism, thus improve insulin sensitivity without increased insulin secretion. many insulin sensitizing drugs are pparg agonists. n subfamily group d, member also known as rev-erba activates histone deacetylation, thereby regulating gene expression. publications indicate that snps in these nuclear receptors associate with obesity and/or diabetes. our data suggest that decreased expression of hnf4a, nr3c esr ar, pparg and nr1d overexpression contribute to t2dm. in gk rats, some compensational pathways still exist, for example a no synthesis pathway is up-regulated. three genes nitric oxide synthase, argininosuccinate synthetase, and nadh: quinone oxidoreductase related to this pathway are overexpressed. it is well-known that no decreases blood pressure and promotes vascular actions of insulin. nos catalyzes arginine, oxygen and nadph to no and citrulline. ass and nq01contribute to this metabolism cycle. many cytokines increase no regeneration several folds. increased no synthesis pathway indicates an inflammation environment in the liver in gk rats. because reduced cell no action has been reported in diabetes, the beneficial effects of increased no production is uncertain. data analyze reveal increased insulin resistance, hypertension and apoptosis are important to push diabetes to next stage. protein kinase c alpha is mostly expressed in hepatocytes promoting glycogenolysis and gluconeogenesis. activation of prkca mediates serine/threonine phosphorylation of the insulin receptor resulting in decreased active form of insulin receptor, inducing insulin resistance. angiotensin i converting enzyme is an exopeptidase that catalyses angiotensin peptides and has opposite effects on ras axis. thus decreased expression levels of ace accelerate the pathologic process such as hypertension, inflammation, fibrosis and inflammation. gap junction alpha- also known as connexin- is a component of gap junctions providing a route for cell to cell communication via diffusion materials. decreased gja expression particularly in hyperglycemia accelerates apoptosis. • advantages of network screening over single gene based method when comparing our results to the original study conducted by dr. almon, network screening is clearly superior to the single gene-based analysis. one good example is to explain how liver insulin resistance develops. ir is the major character of t2dm and also present in gk rats after weeks of age. in the original study, authors notice higher expression of p thus suspecting interaction of p with irs leading to ir. however, we believe that the developing ir is a dynamic process involving many steps. the first step could be significantly decreased igf-1r expression after weeks inducing ir in gk. after that, higher expression of ctsd accelerates ir. compensational pathways also occur, which includes irs overexpression at 8-12w in gk. however as pkc overexpression plus decreased expression of many nuclear factors such as pparg at 16-20w, ir deteriorates and diabetes becomes un-returnable. our method is based on the networks and is very different from the gene-based method of identifying the differential expression. discussion and  t2dm is a complex disease, which is usually not caused by individual gene changes, thereby requiring systems biology methods to understand their mechanisms. in this work, we have performed comprehensive active regulatory network survey by network screening to the published gk vs. wky liver microarray data. available resources from msigdb and transfac are combined together to identify the significant pathways responsive to the status of diabetes or normals. after combining the networks according to features or time points, we built functional or time series tf regulatory network graphs. analyzing the graphs reveals: more pathways are active during inter-middle stage diabetes; inflammation, hypoxia, increased apoptosis, decreased proliferation, and altered metabolism are characteristics in gk strain, and displayed as early as 4w. diabetes progression accompanies insults and compensations. nuclear receptors work in concert to maintain normal glycemic robustness system. network-based analysis based on high throughput data is a challenging issue, which is expected to help us understand complex disease such as diabetes and further elucidate the essential mechanisms of living organisms which would escape conventional single gene-based analysis. in this paper, instead of picking up differently expressed genes from high-throughput data, we use known functional pathways to screen datasets and evaluate significantly activated pathways. then genes with no annotated linkages to tf are overlooked and the available gene regulatory relationships are integrated to form a comprehensive tf regulatory network, which cannot be achieved by single gene based method. the network shows a whole picture of activated tf regulated functional gene sets under certain conditions and is much easier to bring the biological insights to us. to our knowledge, two conclusions have not been reported before. the first one comes out from tf regulatory network at 4w gk. it is well-known that the major cause of diabetes in gk rats is insulin secreting beta cell dysfuction. beta cell mass in gk is only half of that in wky after birth. to be surprised, we find that at very early age liver already exhibits serious gene expression alternations involving in bile metabolism dysfunction, inflammation, increased apoptosis and decreased proliferation, which greatly contribute to diabetes development. another interesting finding is that the nuclear receptors working in concert to maintain robustness of normal blood glucose. although the relationships of those nuclear receptors with diabetes have been investigated individually before, it is the first time to report how they work together as a fine tune. restoring their network regulation may have important therapeutic potentials. this is the first time to use network screening to explain the role of liver in development of diabetes and the underline mechanism. the results provide many important rational information and insights into guiding experiments design. it is worth pointing out that the molecular relationships change dynamically, depending on the conditions in a living cell, which suggests implicitly that all of the relationships in the knowledge-based network do not always exist. note that some methods are proposed for identifying the active networks from measured data. our method evaluates the networks from only one set of data measured under one condition to estimate the absolute consistency between network structure and the data, while the other methods generally need the two sets of data to estimate their relative difference by some criteria such as mutual information. we combined various resources together to identify the significant regulatory networks related to the development stages of diabetes. the matching between networks and gene expression profiling was identified by the evaluation of network screening. the active regulatory networks are the potential disease signatures from the comparison of gk and wky rats. the dynamics of regulatory networks indicate the dysfunctional progression from the network perspective. in conclusion, network screening is a superior approach to analyze complex disease such as diabetes. the conclusions drawn from this method are more complete and systemic, which gives biologist better guidance for further experiment design. actually, we are now extending this approach for screening general biomolecular networks with both directed and undirected edges, and in future possibly for studying the problem of networkomics which covers all stable forms of biomolecular networks not only at different biological conditions but also at different spatiotemporal situations. abbreviations t2dm: type diabetes mellitus; gk: goto-kakizaki; wky: wistar-kyoto; igt: impaired glucose tolerance stage; ift: impaired fasting glucose stage; gcp: graph consistency probability; tf: transcriptional factor; msigdb: molecular signatures database; fdr: false discovery rate; dag: directed acyclic graph; gn: gaussian network; geo: gene expression omnibus; mody: maturity-onset non-insulin-dependent diabetes of the young competing interests the authors declare that they have no competing interests. authors' contributions hz, kh and lc conceived the research. hz, ss and zpl performed the study. gp and jw gave valuable suggestions and improvements. lc and hz supervised the project. hz and zpl drafted a version of the manuscript. all authors wrote and approved the manuscript. 
in the aim of identifying significant transcriptional regulatory networks in the liver contributing to diabetes, we have performed comprehensive active regulatory network survey by network screening in weeks, 8- w, and 18- w goto-kakizaki rat liver microarray data. we identify active regulatory networks in gk rat by network screening in the following procedure. first, the regulatory networks are compiled by using the known binary relationships between the transcriptional factors and their regulated genes and the biological classification scheme, and second, the consistency of each regulatory network with the microarray data measured in gk rat is estimated to detect the active networks under the corresponding conditions. our findings also implicate that network screening is able to help us understand complex disease such as diabetes, and demonstrate the power of network systems biology approach to elucidate the essential mechanisms which would escape conventional single gene-based analysis. the comprehensive survey of the consistency between the networks and the measured data by the network screening approach in the case of non-insulin dependent diabetes in the gk rat reveals: more pathways are active during inter-middle stage diabetes; inflammation, hypoxia, increased apoptosis, decreased proliferation, and altered metabolism are characteristics and display as early as 4weeks in gk strain; diabetes progression accompanies insults and compensations; nuclear receptors work in concert to maintain normal glycemic robustness system. although gene expression levels are able to be measured via microarray since it is difficult to evaluate the contributions of one altered gene expression to a specific disease. notably this is the first comprehensive network screening study of non-insulin dependent diabetes in the gk rat based on high throughput data of the liver.
 type diabetes mellitus is a complex systemic disease, with significant disorders of metabolism. the liver, a central energy metabolic organ, plays a critical role in the development of diabetes. although gene expression levels are able to be measured via microarray since it is difficult to evaluate the contributions of one altered gene expression to a specific disease. one of the reasons is that a whole network picture responsible for a specific phase of diabetes is missing, while a single gene has to be put into a network picture to evaluate its importance. in the aim of identifying significant transcriptional regulatory networks in the liver contributing to diabetes, we have performed comprehensive active regulatory network survey by network screening in weeks, 8- w, and 18- w goto-kakizaki rat liver microarray data.  we identify active regulatory networks in gk rat by network screening in the following procedure. first, the regulatory networks are compiled by using the known binary relationships between the transcriptional factors and their regulated genes and the biological classification scheme, and second, the consistency of each regulatory network with the microarray data measured in gk rat is estimated to detect the active networks under the corresponding conditions. the comprehensive survey of the consistency between the networks and the measured data by the network screening approach in the case of non-insulin dependent diabetes in the gk rat reveals: more pathways are active during inter-middle stage diabetes; inflammation, hypoxia, increased apoptosis, decreased proliferation, and altered metabolism are characteristics and display as early as 4weeks in gk strain; diabetes progression accompanies insults and compensations; nuclear receptors work in concert to maintain normal glycemic robustness system.  notably this is the first comprehensive network screening study of non-insulin dependent diabetes in the gk rat based on high throughput data of the liver. several important pathways have been revealed playing critical roles in the diabetes progression. our findings also implicate that network screening is able to help us understand complex disease such as diabetes, and demonstrate the power of network systems biology approach to elucidate the essential mechanisms which would escape conventional single gene-based analysis. 9- september the 4th international conference on computational systems biology suzhou, p. r. china 