Pharmacodynamic/Pharmacogenomic Modeling of Insulin Resistance Genes in Rat Muscle After Methylprednisolone Treatment: Exploring Regulatory Signaling Cascades

Corticosteroids (CS) effects on insulin resistance related genes in rat skeletal muscle were studied. In our acute study, adrenalectomized (ADX) rats were given single doses of 50 mg/kg methylprednisolone (MPL) intravenously. In our chronic study, ADX rats were implanted with Alzet mini-pumps giving zero-order release rates of 0.3 mg/kg/h MPL and sacrificed at various times up to 7 days. Total RNA was extracted from gastrocnemius muscles and hybridized to Affymetrix GeneChips. Data mining and literature searches identified 6 insulin resistance related genes which exhibited complex regulatory pathways. Insulin receptor substrate-1 (IRS-1), uncoupling protein 3 (UCP3), pyruvate dehydrogenase kinase isoenzyme 4 (PDK4), fatty acid translocase (FAT) and glycerol-3-phosphate acyltransferase (GPAT) dynamic profiles were modeled with mutual effects by calculated nuclear drug-receptor complex (DR(N)) and transcription factors. The oscillatory feature of endothelin-1 (ET-1) expression was depicted by a negative feedback loop. These integrated models provide testable quantitative hypotheses for these regulatory cascades.


Introduction
A major function of the hypothalmic/pituitary/adrenal axis (HPA-axis) involves regulating the production, storage, use, and distribution of substrates for systemic energy metabolism. Glucocorticoids, the effecter hormones from the adrenal gland, were named for their glucose-regulating properties and have extensive impact on carbohydrate, lipid, and protein metabolism. Glucocorticoids are the major regulators of systemic gluconeogenesis, a multi-tissue phenomena involving substrate fl ow from muscle and fat coupled with enzymatic synthesis of glucose from these substrates in liver and kidney. Glucocorticoid output from the adrenal cortex, which occurs with a circadian periodicity, serves to maintain adequate blood glucose in order to assure that suffi cient glucose is available to fulfi ll the needs of vital organs, particularly the brain. This is accomplished by adjusting both supply and demand. On the supply side, glucocorticoids increase the synthesis and activity of enzymes in the gluconeogenic pathway in liver and kidney, as well as ancillary enzymes/pathways such as expression of aminotransferases and components of the urea cycle. They also control the release of gluconeogenic substrates from muscle and fat. On the demand side, they decrease glucose utilization in peripheral tissues which in part encompasses their anti-infl ammatory and immunosuppressive actions.
Synthetic glucocorticoids, corticosteroids (CS), are used therapeutically for a wide variety of conditions that require immune and/or infl ammatory modulation. Because corticosteroids pharmacologically magnify the physiological actions of glucocorticoids, therapeutic use of this class of drugs is accompanied by a wide range of adverse effects that include muscle wasting and hyperglycemia. Skeletal muscle is the major organ responsible for glucose disposal, accounting for about 80% of total glucose uptake under insulin stimulation. Corticosteroids exert their glucose sparing effect on skeletal muscle mainly by inhibiting glucose utilization and switching muscle energy metabolism from glucose to lipids. Inhibition of insulin-directed glucose disposal in skeletal muscle in conjunction with enhanced gluconeogenesis results in elevated plasma glucose levels (diabetes). In addition, prolonged exposure to corticosteroids and the attendant net degradation of muscle protein causes atrophy of the musculature.
Previously, we used a rat model and employed Affymetrix gene arrays to profi le the time course of changes in gene expression in skeletal muscle of male rats in response to a single bolus dose of the synthetic corticosteroid, methylprednisolone (MPL). ADX animals were used to eliminate circadian oscillation in corticosteroid responsive genes, which facilitated data mining. This time course involved 16 time points over a 72 hour period plus untreated controls. Because these experiments were initiated using adrenalectomized animals, the drug in effect acted as a stimulus that perturbed the homeostatic balance of the system, and the experiment monitored the deviation of the system and its return to the original state. The objective was to identify those changes in transcription induced by corticosteroids, which if perpetuated by repeated dosing, would manifest as chronic adverse effects. Using a fi ltering approach to eliminate unregulated probe sets, we identifi ed 653 probe sets out of 8799 with a high probability of being regulated (Almon et al. 2004). Of these, we identifi ed 29 probe sets (representing 22 gene transcripts) directly related to insulin resistance (Almon et al. 2005b).
Although very useful, a single time series only provides a one dimensional view of the dynamics of the system in response to the stimulus. A pharmacological time series is different from most time series studies (for example those assessing developmental changes) in that it can be challenged using a different dosing regimen. A second dosing regimen is valuable in two regards. First, it can serve to corroborate results of the fi rst dosing regimen. If changes in expression levels for a particular transcript are observed in multiple time points and in multiple animals from two different doses, these changes are unlikely to result from chip artifacts. Second, the results can be used to group genes with common mechanisms of regulation. If two or more genes have a common mechanism of regulation, then their response profi les should follow the same pattern regardless of the dosing regimen. More recently, we used microarrays to broadly characterize the response of skeletal muscle to a second dosing regimen which entailed chronic infusion of MPL in a group of ADX male rats. For these experiments MPL was infused at 0.3 mg/kg/h for 168 hours using Alzet osmotic pumps. Four animals were sacrifi ced at each of ten different times during the infusion period. Here the drug essentially was used as an unbalancing stimulus and the experiment evaluated the capacity of the system to rebalance in the continuous presence of the drug.
Microarrays provide a method of high throughput data collection that is necessary for constructing comprehensive information on the transcriptional basis of such complex systemic polygenic phenomena as insulin resistance. When microarrays are used in rich in vivo time series studies they yield temporal patterns of changes in gene expression that illustrate the cascade of molecular events in time that cause complex responses. With such data, the challenge then becomes one of constructing rational, quantitative, mechanism-based frameworks that describe the relationships between the elements of the cascades. Kinetic/dynamic modeling provides an approach to developing quantitative, testable, mechanism-based hypotheses concerning the relationship between drug kinetics and elements of cascades that continue long after the drug has dissipated. Such models can accommodate hierarchal cascades where one process generates molecules that become effectors or mediators for other processes. They can also accommodate convergence of cascades that commonly occur in the control of the expression of genes where binding sites for multiple nuclear factors participate in the regulation of the level of expression of a particular mRNA. Proteins (transcription factors) expressed from this mRNA in turn become endogenous mediators that may themselves become effector molecules for other processes.
In the present report we use the two response profi les (acute and chronic) to develop kinetic/ dynamic models for the regulation of the expression of six genes that contribute to corticosteroid induced insulin resistance in skeletal muscle. The fi rst is insulin receptor substrate-1 (IRS-1). Here we model the expression profi le as being driven by both MPL and Interleukin 6 (IL-6). The second and third are two nuclear encoded mitochondrial genes with similar expression profi les, uncoupling protein 3 (UPC3) and pyruvate dehydrogenase kinase isoenzyme 4 (PDK4). These two genes are modeled as being driven by both MPL and the muscle specifi c transcription factor MyoD. The fourth and fifth are two genes related to the increased use of lipids for energy, fatty acid translocase (FAT) and glycerol-3-phosphate acyltransferase (GPAT). These two genes are modeled as being under the combined infl uence of both MPL and sterol regulatory element binding protein-1c (SREBP-1c). The last is endothelin-1 (ET-1) which infl uences vascular resistance in skeletal muscle, thereby infl uencing both blood fl ow and glucose disposal. Although these are four separate models, they represent a beginning at developing more comprehensive models for the complex phenomena of glucocorticoid-induced insulin resistance in skeletal muscle. Such models gain more importance when the data linking glucocorticoids to the insulin resistance associated with metabolic syndrome are considered.

Experimental
Gastrocnemius muscles were obtained from animal studies performed in our laboratory (Sun et al. 1999;Ramakrishnan et al. 2002). In the acute study, male ADX Wistar rats weighing 225-250 g (Harlan Spague-Dawley Inc., IN) were subject to right external jugular vein cannulation under light ether anesthesia one day prior to the study. A single intravenous bolus of 50 mg/kg MPL (Solu-Medrol, Pharmacia-Upjohn Company, MI) was given to 47 animals via right jugular vein cannula. Rats were sacrifi ced by aortic exsanguination at 0. 25, 0.5, 0.75, 1, 2, 4, 5, 5.5, 6, 7, 8, 12, 18, 30, 48 and 72 h after dosing. The sampling time points were selected based on previous studies of receptor dynamics and enzyme induction in muscle and liver (Sun et al. 1998;Sun et al. 1999;Ramakrishnan et al. 2002). Four vehicle-treated rats were designated as controls and were considered as sacrifi ced at time point zero. In a companion study, 6 animals were treated with single intravenous bolus of 50 mg/kg MPL. Blood samples were taken serially from the cannula at 0.5, 1, 2, 4, 5, 6, 8, 10, 12, 24, 36, 48 and 72 h. Blood samples were stored at −80 °C and were used to determine plasma MPL concentrations. In our chronic study, 40 male ADX Wistar rats weighing 300-350 g were given 0.3 mg/kg/h infusions of MPL reconstituted in supplied diluent. The infusions were given using Alzet osmotic pumps (Model, 2001, fl ow rate 1 ul/h; Alza, Palo Alto, CA). The pump drug solutions were prepared for each rat based on its predose body weight. Pumps were subcutaneously implanted between the shoulder blades on the back. Four rats were sacrifi ced by aortic exsanguination at 6, 10, 13, 18, 24, 36, 48, 72, 96, and 168 h. Four control rats were implanted with a saline-fi lled pump and sacrifi ced at various times throughout the 7-day study period. Gastrocnemius muscles were rapidly excised, quickly frozen in liquid nitrogen, and stored at −80 °C.

Drug assay
MPL plasma concentrations were measured by a previously described normal-phase high performance liquid chromatography (HPLC) method (Haughey and Jusko, 1988). The quantitation limit of this method is 10 ng/ml.

RNA extraction
Frozen gastrocnemius muscles were ground into powder using a mortar and pestle chilled with liquid nitrogen. Total RNA extractions were carried out by a Trizol-chloroform based extraction method. About 100 mg muscle powder from each animal was added to prechilled Trizol Reagent (Invitrogen, Carlsbad, CA) at a ratio of tissue: Trizol of 1:10. Extractions were performed according to manufacturer protocols. Extracted RNA was further purifi ed by passage through RNUeasy mini-columns (Qiagen, Valencia, CA). Final extracted RNA samples were resuspended in nuclease-free water. Quantity of total RNA was determined by spectrophotometry and purity was assessed by agarose gel electrophoresis. Extracted total RNA preparations were stored at −80 °C.

Microarrays
Isolated muscle RNA from each individual animal was used to prepare biotinylated cRNA target according to manufacturer protocols. Then biotinylated cRNAs from the acute study were hybridized to 51 individual Affymetrix GeneChip ® Rat Genome U34A (Affymetrix, Santa Clara, CA), which contained 8799 probe sets. The target cRNAs from the chronic study were hybridized to 44 individual Affymetrix GeneChip ® Rat Genome 230A (Affymetrix, Santa Clara, CA), which contained 15967 probe sets. Oligonucleotide microarrays were utilized because of their high reproducibility between separate arrays. Both datasets have been submitted to the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus database (GSE490 and GSE5101) and are also available online at http:// pepr.cnmcresearch.org/ . Our approach to identifying probes of interest has been described in previously published articles on mining the datasets for muscle, liver and kidney from the acute MPL treated animals (Almon et al. 2004;Almon et al. 2005a;Almon et al. 2005c). Literature searches of the 653 probe sets identifi ed in the acute skeletal muscle dataset yielded 29 probe sets representing 22 genes whose expression is related to insulin resistance in skeletal muscles (Almon et al. 2005b). Suffi cient information from our own work and from the literature allowed us to develop four rational mechanistic models describing the regulation of the expression of six genes related to insulin resistance. Data from both acute and chronic studies were simultaneously modeled.

Kinetics
The kinetics of plasma MPL following intravenous drug injection and infusion was described by a two-compartment model with a zero-order input k 0 into the central compartment.
where (2) where A p and A t are the drug amounts in the plasma and tissue compartments and CL, V p , k 12 and k 21 are clearance, central compartment volume of distribution, and distribution rate constants. These parameters were fi xed according to previous literature values (Sun et al. 1999).

Receptor dynamics
Glucocorticoid receptor dynamics in skeletal muscle following single intravenous 50 mg/kg MPL treatment was previously described by a mechanistic receptor-gene mediated pharmacodynamic model ( Fig. 1) (Sun et al. 1999). The model was developed using the assumption that unbound MPL molecules diffuse through cell membranes and bind with cytosolic free receptors. Drug-receptor complex is translocated into the nucleus where it dimerizes and binds to glucocorticoid response elements (GREs) in the target DNA. The binding of drug-receptor complex to GREs enhances or inhibits the expression of target genes. Corticosteroids are known to inhibit the expression of their own receptors (Dong et al. 1988). After dissociation from DNA, GR receptors are recycled into cytosol, where part of receptors are degraded and the rest may be further activated by MPL (Sun et al. 1999). The equations describing the model are: where symbols represent cytosolic free glucocorticoid receptor density (GR), GR mRNA (GR mRNA ), cytosolic GR-drug complex (DR), GR-drug complex in nucleus (DR(N)), zero-order GR mRNA synthesis rate constant (k s_Rm ), and fi rst-order rate constants for GR mRNA degradation (k d_Rm ), GR synthesis (k s_R ) and GR degradation (k d_R ). Other rate constants include: second-order association rate constant of GR and drug (k on ), translocation rate constant of DR from cytosol to nucleus (k T ), and the overall turnover rate of DR(N) (k re ). The fraction of GR that is recycled from nucleus to cytosol (R f ) can be further activated by association of drug. Thus (1-R f ) is the fraction of GR that is degraded during one cycle of drug-receptor complex translocation. The IC 50_Rm is the concentration of DR(N) at which the synthesis rate of GR mRNA is reduced to 50% of its baseline level. Assuming the receptor system is at steadystate before administration of MPL and there is no endogenous glucocorticoid present in ADX rats, the following conditions are derived from Eq. (3) and (4).
where GR 0 mRNA and GR 0 are the baseline levels of GR mRNA and free cytosolic GR at time zero. In our previous study, these were measured in tissues from the same animals by quantitative Northern hybridization and [ 3 H]-dexamethasone ligand binding assays (Sun et al. 1999). In the present report, GR 0 mRNA and GR 0 were fi tted as well as other parameters in the receptor dynamics model. Values of DR 0 and DR(N) 0 were assumed to be zero in the absence of steroid. Parameters estimated from the modeling of acute data were used to simulate receptor mRNA and receptor density profi les for the chronic study.
Insulin receptor substrate-1 (IRS-1) IRS-1 belongs to the IRS family of adaptor proteins and is involved in insulin signaling. It has been reported that chronic CS treatment down-regulates the expression of IRS-1 in skeletal muscle (Giorgino et al. 1993). Additionally, interleukin 6 (IL-6) exerts long-term inhibitory effects on the transcription of IRS-1 while IL-6 receptor expression is enhanced by CS treatment (Rotter et al. 2003). A model is proposed according to these regulatory pathways (Fig. 2).
The equation describing the up-regulation of IL-6 receptor type 1 mRNA (IL6R1 m ) time profi le following MPL treatment is: where rate constants represent zero-order synthesis rate (k s_IL6R1m ) and fi rst-order degradation (k d_IL6R1m ) of IL-6 receptor 1 mRNA. The enhancement of gene transcription is dependent on DR(N) with a maximum stimulatory effect (S max_IL6R1m ) and SC 50_IL6R1m is the concentration of DR(N) that leads to 50% of maximum stimulation of mRNA transcription. The baselines of IL6R1m were fi xed to 1 in the modeling for both acute and chronic studies.
The IRS-1 mRNA (mRNA) was described as: where symbols represent number of transit compartments (N), the nth transit compartment (TC n ), the fi rst-order rate constant describing transit compartment turnover (k e ), and the zero-order synthesis rate (k s_m ) and fi rst-order degradation rate (k d_m ) of mRNA. The IC 50_DR(N) represents the concentration of DR(N) at which the production rate of mRNA is reduced to 50% of its baseline level under the inhibition of DR(N). The IC 50_IL6R1 refl ects the concentration of TC n at which the production rate of mRNA is one-half of its baseline level when it is only inhibited by IL-6. Assuming steady-state at time zero, Eq. 12 allows calculation of: where TC n 0 is the baseline value of TC n and is fi xed as 1.
Uncoupling protein 3 (UPC3) and pyruvate dehydrogenase kinase isoenzyme 4 (PDK4) Dexamethasone leads to increased expression of UCP3 in skeletal muscle in rats and mice (Gong et al. 1997;Sun et al. 2003). Similarly, PDK4 mRNA and protein are augmented following dexamethasone treatment and are associated with decreased glucose oxidation rates (Qi et al. 2004). It has been reported that the expression of UCP3 in muscle is enhanced by MyoD (Solanes et al. 2000;Solanes et al. 2003). Examination of the promoter regions of both UCP3 and PDK4 reveals several GREs and MyoD binding sites that are conserved in several species. Models were thus proposed based on assumptions that both UCP3 and PDK4 are regulated directly by CS and indirectly by MyoD (Fig. 3). The equation describing MyoD mRNA (MyoD m ) profi le is: where zero-order synthesis rate (k s_MyoDm ) and fi rstorder degradation rate constants (k d_MyoDm ) of MyoD mRNA occur. An intermediate biosignal (BS) is described by a production and degradation rate constant of k BS (Jin et al. 2003). The inhibition of gene transcription is dependent on DR(N) and BS with inhibitory coeffi cients (IC 50_MyoDm and IC 50_BS ) refl ecting corresponding concentrations leading to 50% mRNA transcription reduction. The baselines of MyoD m in acute (BmRNA a_MyoD ) and chronic (BmRNA c_MyoD ) studies were estimated in the modeling.
The UCP3 mRNA (mRNA) time profi les were described as: The PDK4 mRNA (mRNA) time profi les were described as: where γ MyoD is a power term representing the stimulatory effi ciency on UCP3 or PDK4 gene transcription by MyoD. The transcription of PDK4 is dependent on DR(N) with a maximum stimulatory effect (S max ) and a stimulatory constant (SC 50 ). For UCP3, assuming SC 50 is much higher than DR(N) concentrations, S max and SC 50 are combined as one parameter S DR(N) , which represents the linear coeffi cient describing the stimulatory effi ciency of DR(N) on UCP3 mRNA production rate. Assuming steady-state at time zero, Eq. 18 and 19 allow calculation of: where TC n 0 is the baseline value of TC n and is fi xed as 1.
Fatty acid translocase (FAT) and glycerol-3-phosphate acyltransferase (GPAT) These two genes are related to carbohydrate and lipid metabolism and are involved in insulininduced metabolic actions. They are all stimulated . Pharmacodynamic/pharmacogenomic model for CS effects on gene expression of MyoD and its regulated genes. Symbols and differential equations are defi ned in Eq. 14 to 18 and Table 3. by sterol regulatory element binding protein-1c (SREBP-1c) (Ericsson et al. 1997;Ntambi, 1999;Nguyen et al. 2000;Almon et al. 2005b). Models are proposed to describe their time profi les according to underlying mechanisms (Fig. 4). The equations describing SREBP-1c mRNA (SREBP m ) time profi le following MPL treatment are: where symbols represent zero-order synthesis rate (k s_SREBPm ) and first-order degradation rate (k d_SREBPm ) of SREBP-1c mRNA. The IC 50_SREBPm and IC 50_BS represent the concentrations of DR(N) and BS at which the corresponding production rates are reduced to 50% of baseline values. The baselines of SREBP m in acute (BmRNA a_SREBP ) and chronic (BmRNA c_SREBP ) studies were either fi xed to 1 or estimated in the modeling. The FAT and GPAT mRNA (mRNA) versus time profi les were described as: where IC 50_DR(N) represents the concentration of DR(N) at which the production rate of mRNA is reduced to 50% of its baseline level and γ SREBP is a power term representing the stimulatory efficiency on gene transcription by SREBP-1c. FAT and GPAT mRNA data from both i.v. bolus and infusion studies were simultaneously fitted. Assuming steady-state at time zero, Eq. 25 allows calculation of: where TC n 0 is the baseline value of TC n and is fi xed as 1. Endothelin-1 (ET-1) It has been found that ET-1 may enhance nitric oxide production via its receptor located in endothelial cells (Namiki et al. 1992). Conversely, nitric oxide can inhibit ET-1 synthesis in different cell types (Takada et al. 1996). A model is proposed assuming there is a feedback loop regulating the production of ET-1 (Fig. 5). The equations describing ET-1 mRNA (mRNA) time profi le following MPL treatment are: where gene transcription is dependent on DR(N) with a linear stimulatory coeffi cient (S DR(N) ). The IC 50_ET1 represents the concentration of TC n at which the production rate of mRNA is reduced to 50% of its baseline level under the inhibition of itself and γ is a power term describing the amplifi cation effect during the transit steps. Assuming steady-state at time zero, Eq. 27 allows calculation of: where TC n 0 is the baseline value of TC n and is fi xed as 1.

Data analysis
The gene array data were obtained from a 'giant rat' study design refl ecting naïve pooling of data from a group of animals (Sun et al. 1999;Ramakrishnan et al. 2002). Microarray data from both i.v. bolus and infusion studies were simultaneously fi tted. Since acute and infusion data were obtained from two different types of chips, a scaling factor (SF) was incorporated in the modeling to account for the different sensitivity of two chips to the same gene expression changes. The output functions of mRNA in acute (Y A (t)) and chronic (Y C (t)) studies were expressed as: where mRNA A (t) and mRNA C (t) represent mRNA expression at time t in acute and chronic studies.
Since relative sensitivity was considered in the present report, the scaling factor was only added to the chronic data.
Assuming that the errors from the observed and predicted responses are normally distributed, the ADAPT II program (D'Argenio and Schumitzky, 1997) with the maximum likelihood method was applied for all fi ttings. The following variance model was used: (33) Figure 5. Pharmacodynamic/pharmacogenomic model for CS effects on gene expression of ET-1. Symbols and differential equations are defi ned in Eq. 25 to 28 and Table 5.
where V i is the variance of the response at the ith time point, t i is the actual time at the ith time point, θ represents the systemic parameter vector for the respective pharmacodynamic (PD) model, σ is defi ned as the vector of variance parameters for the model, Y(θ, t i ) is the predicted response value at time t i from the model. Variance parameters σ 1 and σ 2 were estimated along with model parameters during fi ttings. Mean data from 3 or 4 animals at each time point were used in fi tting the models. The goodness-of-fi t criteria included visual inspection of the fi tted curves, estimator criterion value, sum of squared residuals, Akaike information criterion, Schwartz criterion, and coeffi cients of variation (CV) of the estimated parameters.

Drug kinetics
Plasma concentrations of MPL following i.v. bolus were described by a two-compartment model. The MPL PK profi le following drug infusion was generated using the reported parameters (Ramakrishnan et al. 2002) (Fig. 6). The estimated and fi xed parameters are summarized in Table 1.

Receptor dynamics
The integrated receptor-gene mediated model adequately captured both GR mRNA and cytosolic free receptor data following acute MPL treatment. The profi les of GR mRNA and cytosolic GR density versus time and the least-squares fi tted curves are shown in Figure 7. The GR mRNA concentration drops from the baseline level to the trough in 10 h and then slowly returns to the baseline. Free GR disappears quickly from cytosol and is not detectable between 0.5 to 2 h after dosing. The recovery of GR density shows an initial rapid phase between 2 to 20 h which is responsible for 40% of recovery and a later slower phase persisting until 72 h. Simulations of major components in the receptor dynamic model (data not shown) suggest that DR(N) increases after MPL treatment and reaches a peak at 0.1 h and declines thereafter. The DR(N) disappears 10 h after dosing. In the present report, the simulated profi le of DR(N) served as the driving force for stimulation or inhibition of downstream gene transcription. Since receptor mRNA and receptor density data were absent for the infusion study, the resultant parameters from the acute study were used to simulate receptor dynamics in muscle during the infusion and receptor dynamics were fi xed in the subsequent gene regulation modeling (Fig. 7). These simulations suggest that both receptor mRNA and free cytosolic receptor are maintained at very low levels during long-term MPL infusion, which is consistent with our measurements in the liver from these same animals (Ramakrishnan et al. 2002).
The fi xed and estimated parameters of the receptor dynamics model are listed in Table 1. This model is over-parameterized leading to the necessity of fi xing the k T value, which was derived from in vitro literature data. Translocation of GR from cytoplasm to nucleus is complete in 10 minutes when cells are incubated with 10 nM dexamethasone Figure 6. MPL pharmacokinetics following 50 mg/kg MPL i.v. bolus administration and 0.3 mg/kg/h 7-day infusion. Chronic PK profi le was simulated (Eq. 1 and 2) using pharmacokinetic parameters listed in Table 1. Original plasma MPL concentrations were obtained from Sun et al. (Sun et al. 1999). (Htun et al. 1996). Dexamethasone at 10 −6 M leads to complete translocation in 0.1 h (Hache et al. 1999). However, the translocation rate in those in vitro studies is determined by both receptor-steroid binding and fi rst-order cytosol-nucleus transloca-tion. According to literature data, the fi rst-order translocation rate k T has to be higher than 35 h −1 . Our k T was fi xed to 90 h −1 in the modeling of receptor dynamics. The GR mRNA 0 and GR 0 were estimated at 2.99 fmol/g tissue and 65.3 fmol/mg protein.  Table 1. Original GR mRNA and cytosolic GR density data were obtained from Sun et al. (Sun et al. 1999

IL6R and IRS-1
Simultaneous modeling using a basic indirect response model with direct stimulation by DR(N) well captured IL-6 receptor 1 mRNA following both i.v. bolus and infusion treatments. Figure 8 shows the normalized data from gene arrays and the fi tted curves. Following 50 mg/kg MPL treatment, expression of IL-6R1 mRNA exhibits a similar pattern as tyrosine aminotransferase (TAT) in the liver (Sun et al. 1998;Jin et al. 2003). It shows a rapid induction of expression which peaked at 5.5 h following drug treatment. The enhanced expression lasts for about 20 h and thereafter maintains baseline level. Compared to DR(N), the induction of IL-6 receptor mRNA reveals a 3 h time delay. Following the 0.3 mg/kg/h infusion, IL-6 expression is increased and maintained at a high level throughout the long-term treatment period. Table 2 lists the estimated PD parameters and their precision. The estimated k d_IL6R1m yields a degradation half-life of 2.26 h. The stimulatory constant SC 50_IL6R1m was estimated at 1.09 fmol/mg protein, suggesting high sensitivity to MPL action.
The transcription of IRS-1 was described by a more complex regulation with inhibitory effects derived directly by DR(N) and indirectly through IL-6 receptor. The proposed PK/PD model reasonably captured the down-regulation of IRS-1 message by MPL administration (Fig. 8). In the acute study, the expression profi le exhibits two inhibitory phases. The fi rst decline occurs immediately after drug treatment and reaches the nadir around 5.5 to 6 h. This phase was successfully described by direct inhibition by DR(N) under the assumption that the drug-receptor complex represses the transcription of IRS-1 mRNA. It returns to a value which is close to baseline by 12 h. The second decline phase starts at 18 h and exhibits a much slower decrease with a trough at 40 h. The second phase was described by augmented inhibition of IRS-1 mRNA transcription by the IL-6 receptor. Comparing the curves of IRS-1 and IL-6 receptor 1 mRNA reveals a delay time of 30 h between the second decline of IRS-1 and its driving force. This is consistent with a previous report that reduction of IRS-1 mRNA is seen 24 h following IL-6 treatment in adipose cells (Rotter et al. 2003). The long time delay necessitated the utilization of multiple transit compartments between the driving force and the regulated gene expression. Five transit compartments were adequate to capture the curve and give reasonable parameter estimates and CV values (Table 2). In the infusion study, IRS-1 mRNA declines following drug treatment and reaches a steady-state throughout the 7-day study.

MyoD, UCP3 and PDK4
MyoD mRNA expression profi les following i.v. bolus and infusion treatments were well captured by the dynamic model assuming a direct inhibition of transcription by DR(N) and another inhibitory effect from the biosignal (Fig. 9). In the acute study, MyoD mRNA exhibits a sharp decline immediately following MPL treatment. The trough occurs at 4 h, only 2 h after the peak of DR(N). Following infusion, MyoD transcription quickly declined and reached the fi rst trough at 6 h and was reduced to almost zero. Then it increased with a slow slope and regained 40% of its basal level at 36 h and then decreased thereafter. The second decline is described by the inhibitory effect from the biosignal. The estimated parameters are listed in Table 3. A high k d_MyoDm yields a short degradation half-life of 0.495 h −1 and is consistent with the rapid decline and short time delay. The parameters were estimated with reasonable precision with most CVs below 30% (Table 3). The proposed model was able to describe the time profi les of both PDK4 and two probe sets for UCP3. The original data and the fi tted curves are shown in Figure 9. In the acute study, both UCP3 and PDK4 mRNA reveal a sharp enhancement of production which is well captured by the direct stimulation of DR(N). The peaks occur at 5 h after drug administration. All three probe sets return to baseline by 15 h and thereafter continuously decline and reach steady-states after 18 h. The inhibitory phases were adequately captured by the reduced activation by MyoD. Following drug infusion, all three probes exhibit only enhanced transcription without the inhibitory phase. One transit compartment was adequate for modeling of both UCP3 and PDK4. The degradation half-life of UCP3 was estimated between 3.6 and 4.4 h while that of PDK4 was 2.1 h. Moderate stimulatory coefficients (S DR(N) ) around 0.15 (fmol/mg protein) −1 were obtained for UCP3. For PDK4, two parameters S max and SC 50 were required to capture both acute and chronic profi les simultaneously. The estimated low SC 50 suggests high efficiency of the stimulatory effect and also justifi es the requirement of using both parameters. Two probe sets for UCP3 exhibit similar expression versus time patterns and yield similar parameter estimates, suggesting good reproducibility of the oligonucleotide microarray assay. Most parameters were estimated with CV below 40% (Table 3).

SREBP-1c, FAT, and GPAT
The SREBP-1c mRNA time profi les were captured by the joint inhibitory effects on the transcription directly by DR(N) and indirectly by an unknown intermediate biosignal (Fig. 10). The expression of SREBP-1c is down-regulated after MPL treatment. Following its trough at 4 h in the acute study, it returned to baseline in a biphasic pattern with an initial fast return phase and a later shallow phase. A similar pattern was shown in the chronic study. The biphasic change suggests that two driving forces probably are involved. The mediator BS was incorporated in the model based on the assumption that the absolute change of a regulator level by DR(N) may elicit an inhibitory effect on the transcription rate of SREBP-1c mRNA (Jin et al. 2003). The estimated k BS of 0.0436 h −1 yields a relatively long mean transit time of 23 h, which is consistent with the prolonged repression of SREBP-1c until 72 h. The degradation half-life of SREBP-1c mRNA was fi xed to 1 to obviate over-parameterization. Both DR(N) and the biosignal display moderate inhibition on SREBP-1c transcription with an IC 50_SREBP of 66.7 fmol/mg protein and an IC 50_BS of 12.4 fmol/mg protein. The estimated parameters are listed in Table 4. The dynamics of SREBP-1c served as the driving force for FAT and GPAT and was fi xed in the following modeling.
Simultaneous modeling using the proposed model well captured the expression of FAT and GPAT (Fig. 10). They were down-regulated following 50 mg/kg MPL with the fi rst troughs between 5.5 to 8 h. The initial declines were well described by the direct inhibition of gene transcription by DR(N) for both genes. The suppressed activation of SREBP-1c on GPAT leads to a second decline while the same regulatory mechanism leads to a prolonged return phase for FAT. Following 7-day MPL infusion, the fi rst trough occurs around 20 h for both genes. Similarly to what was seen in the acute study, a second decline occurs only with GPAT but not with FAT. The degradation rate constants of the messages are around 0.16 h −1 . The DR(N) has a moderate to high inhibitory effect on the gene transcription with the IC 50_DR(N) between 2.27 to 9.97 fmol/mg protein. GPAT has a higher γ SREBP than that of FAT which may explain the stronger suppression of GPAT message under the control of SREBP-1c. Two to four transit compartments were required to fi t the data. The estimated parameters and their precisions are shown in Table 4.

ET-1
The time profi les of ET-1 mRNA following both acute and chronic treatments were captured by the proposed model (Fig. 11). Immediately following 50 mg/kg MPL dosing, it exhibits a sharp induction peak at 2-3 h and then decreases below baseline after 12 h. The abrupt stimulation and the short time shift between ET-1 mRNA and DR(N) can be explained by a fast mRNA degradation rate constant of 0.958 h −1 , yielding a half-life of 0.723 h. A short half-life might be necessary to obtain quickly induced or reduced effects on vasoconstriction. The power parameter γ was estimated at 8.23, representing a considerable magnifying effect during the signaling transduction steps. A model with 8 transit compartments was selected based on visual inspection of the fi tted curve and other model selection criteria. Following infusion, the change of ET-1 mRNA expression was fl at. This apparent lack of response of ET-1 to chronic infusion might be due to low sensitivity of the probe set on genechip 230A to the change of mRNA level, which can be seen from the estimated SF of 0.118. Table 5 summarizes the estimated parameters and their precision.

Discussion
Insulin resistance is a pathological state where peripheral tissues, particularly skeletal muscle, fail to respond to circulating insulin (Gual et al. 2005).
One of the earliest abnormalities observed in skeletal muscle is the reduced insulin-induced uptake of glucose. The 50 mg/kg MPL i.v. injection in ADX rats leads to a transient increase in the plasma glucose concentration that begins about 2 h and lasts for 9 h (Almon et al. 2005b). Plasma insulin concentrations are increased between 2 to 48 h after dosing, probably due to increased plasma glucose levels .
The development and progression of insulin resistance is not mediated by a single gene. Abnormalities in multiple molecules and tissues are involved in the pathogenesis of the disease. In such situations, research focused on single mediators may distort the true underlying molecular and cellular mechanisms and inadequately reveal overall regulatory factors and pathways. By analyzing joint temporal responses of many genes simultaneously, it should be possible to gain a better understanding of regulatory pathways and their contributions to this complex pathology. By using two different dosing regimens in the time series format, we not only obtained joint confi rmation but also augmented the elucidation of underlying mechanisms of regulation. In the present report we describe the development of four models that encompass six genes associated with insulin resistance. These models include both direct effects of CS on the genes as well as indirect effects due to the effects of CS on other transcription factors such as MyoD and SREBP-1c.
The fi rst model describes the regulation of IRS-1 which plays a central role in insulin signaling. IRS-1 is phosphorylated in response to insulin, insulin like growth factor-1, and cytokines and is preferentially involved in the metabolic actions of insulin (Gual et al. 2005). A decline of IRS-1 amount or function has been linked to decreased glucose uptake in insulin resistant animal models and type II diabetic patients (Sesti et al. 2001;Smith, 2002).
Strong correlations between local and circulating proinfl ammatory cytokines and insulin resistance have been reported (Senn et al. 2002). Among them, IL-6 has the strongest correlation with insulin  resistance and type II diabetes. In liver, IL-6 elevates hepatic glucose output and increases blood glucose by interfering with insulin-induced kinase cascades such as tyrosine phosphorylation of IRS-1 (Senn et al. 2002). Adipose tissue and muscle, two peripheral tissues that are responsible for glucose disposal, are important sites of IL-6 production (Senn et al. 2002). Thus it is possible that IL-6 produced in these tissues may act at the local site and interfere with insulin-induced signaling. It has been found that IL-6 exerts long-term inhibitory effects on the transcription of IRS-1 (Rotter et al. 2003). Although it is well accepted that CS also inhibit the expression of many cytokines including IL-6 in infl ammation, it is not known if CS exhibits similar effects in healthy animals. In our studies, the expression of IL-6 mRNA in muscle did not show signifi cant change. However, we observed a signifi cant up-regulation of IL-6 receptor type 1 transcription. Thus it is reasonable to assume that increased expression of IL-6 receptor potentially enhances muscle sensitivity to IL-6 and augments the inhibitory effect of IL-6 on IRS-1 transcription. The second model describes the regulation of two nuclear encoded mitochondrial genes that have been associated with insulin resistance, UCP3 and PDK4. Uncoupling proteins are a family of transporters that localize in the mitochondrial inner membrane and dissipate the transmembrane potential by transporting protons from the inter-membrane space back into the matrix (Garvey, 2003). UCP3 exhibits limited tissue-specifi c expression confi ned to skeletal muscle in humans. Accumulating evidence indicates that in skeletal muscle, UCP3 contributes to lipid uptake by mitochondria rather than uncoupling oxidative phosphorylation (Garvey, 2003). PDK4 is an enzyme that inactivates the mitochondrial pyruvate dehydrogenase complex. It thus reduces glucose utilization by skeletal muscle, preventing pyruvate (the product of glycolysis) from being used by the mitochondria (Sugden, 2003). The PDK4 gene is preferentially expressed in skeletal muscle. Both genes are associated with enhanced fatty acid oxidation, thus probably are involved in the preferential fuel utilization of lipids instead of glucose in insulin resistance states. Enhanced expressions of both UCP3 and PDK4 in skeletal muscle have been associated with obesity and type II diabetes in humans (Bao et al. 1998;Sugden and Holness, 2002).
The bi-directional profi les of UCP3 and PDK4 mRNA expression following MPL treatment suggest that more than one mediator is involved in the regulatory pathways. The down-regulation of MyoD by CS that we observed in our studies along with literature searches and examination of the promoter regions of the two genes indicate that  myogenic factor MyoD might be involved in a complex signaling network (Solanes et al. 2000;Solanes et al. 2003). MyoD belongs to the basic helix-loop-helix family of DNA-binding transcription factors. It is a master regulator of muscle lineage differentiation and is responsible for the preferential expression of skeletal muscle-specifi c genes (Solanes et al. 2000). At least two phylogenetically conserved MyoD responsive sites were found in the promoter regions of both genes. MyoD is required not only for UCP3 promoter activity but also for its sensitivity to activation by other ligands (Solanes et al. 2000;Solanes et al. 2003). Additionally, the muscle-preferential expression of UCP3 and PDK4 also supports an essential role of MyoD in transcriptional regulation. The third model describes the regulation of the expression of two genes involved in lipid metabolism in skeletal muscle. The fi rst is FAT, also called CD36, which is a membrane protein that facilitates fatty acid uptake and use by skeletal muscle (Heron-Milhavet et al. 2004). Muscle-specifi c over expression of FAT reverses insulin resistance and diabetes in animal models (Heron-Milhavet et al. 2004). The second is GPAT which catalyzes the initial and committed step in the biosynthesis of triglycerides and phospholipids (Ericsson et al. 1997). Expression of these genes indicates a shift of energy consumption in the skeletal muscle from glycolysis towards β-oxidation.
SREBP-1c is a transcription factor of the basic helix-loop-helix/leucine zipper family. It controls expression of genes which are related to adipogenesis and fatty acid metabolism as well as cholesterol metabolism (Nguyen et al. 2000). SREBP-1c is transcriptionally induced by insulin and its enhanced expression mediates the transcriptional effects of insulin in skeletal muscle (Guillet-Deniau et al. 2002). It has been observed that over-expression of SREBP-1c mimics the effects of insulin such as stimulation of glycolytic and lipogenic enzymes. It exerts a pivotal role in long-term muscle insulin sensitivity. The prolonged returning of FAT expression and the second trough of GPAT at 48 h raises the possibility that at least one intermediate biosignal is involved. Both genes are transcriptionally induced by SREBP-1c (Ericsson et al. 1997;Ntambi, 1999;Nguyen et al. 2000;Almon et al. 2005b). The addition of SREBP-1c as an intermediate driving force was able to capture the time profi les of FAT and GPAT following both acute and infusion dosing regimens.
ET-1 is a potent vasoconstrictor. It has been observed that many organs such as heart, lung and skeletal muscles synthesize ET-1, which may act locally to regulate blood fl ow (Sakurai et al. 1991;Guo et al. 1998). The ET-1 detected in this study might be muscle in origin, endothelial cell in origin, or both. Elevated muscle ET-1 may reduce blood fl ow to the musculature and thus reduce glucose disposal in skeletal muscle (Santure et al. 2002). Reciprocal regulation of ET-1 and nitric oxide, a vasodilator have been demonstrated and summarized by Rossi et al. (Rossi et al. 2001). The production of nitric oxide is activated by ET-1 while the elevated nitric oxide level is able to inhibit ET-1 synthesis (Namiki et al. 1992;Takada et al. 1996). This forms a negative feedback loop which may offer an oscillatory feature to ET-1 expression pattern once it is induced or inhibited by drug treatment. The fi tted curve in our acute study exhibited a fl uctuation around the baseline. Other studies demonstrate that ET-1 also exhibits similar time profi les following interleukin-1β and lipopolysaccharide treatments in human endothelial cells with an initial up-regulation and a later down-regulation at 24 h (Zhao et al. 2001;Zhao et al. 2003). Those results suggest that this oscillatory feature of ET-1 is not specifi c to CS treatment. ET-1 expression following drug infusion did not show signifi cant change, possibly due to reduced probe set sensitivity in the chip used for the chronic study.
In the present report, most modeled time profi les feature long time delays between the intermediate biosignal and the regulated changes. One or more transit compartments were incorporated in the models to account for the delay. Following drug treatment, the resultant changes in message have to translate into protein changes before downstream regulatory steps can continue. Other steps such as translocation in or out of the nucleus or induction of other mediators might also be involved.
Time series studies of gene expression following CS treatment shows that many genes such as UCP3, PDK4, and ET-1 exhibit both induction and suppression depending on time. Evaluation of drug effects at a single time point may not reveal the overall effects of the drug. In such situations, time series design and mathematical modeling provide a useful approach to explore the diverse effects at different times and the actual or potential underlying mechanisms.
It is interesting to see that genes with common regulators such as UCP3 and PDK4 display comparable expression versus time patterns. This may