A Model of Glucocorticoid Receptor Interaction With Coregulators Predicts Transcriptional Regulation of Target Genes

Regulatory factors that control gene transcription in multicellular organisms are assembled in multicomponent complexes by combinatorial interactions. In this context, nuclear receptors provide well-characterized and physiologically relevant systems to study ligand-induced transcription resulting from the integration of cellular and genomic information in a cell- and gene-specific manner. Here, we developed a mathematical model describing the interactions between the glucocorticoid receptor (GR) and other components of a multifactorial regulatory complex controlling the transcription of GR-target genes, such as coregulator peptides. We support the validity of the model in relation to gene-specific GR transactivation with gene transcription data from A549 cells and in vitro real time quantification of coregulator-GR interactions. The model accurately describes and helps to interpret ligand-specific and gene-specific transcriptional regulation by the GR. The comprehensive character of the model allows future insight into the function and relative contribution of the molecular species proposed in ligand- and gene-specific transcriptional regulation.


INTRODUCTION
Glucocorticoids contribute to the maintenance of homeostasis in almost all organs and tissues under basal and stress conditions in higher organisms. Many of these homeostatic functions are exerted directly by GC binding to the GR, which regulates the transcription of broad networks of target genes. The GR is a ligand-induced transcription factor that, upon GC binding, translocates to the nucleus and promotes the assembly of multiprotein regulatory complexes at genomic GREs Rogatsky et al., 2003;Luecke and Yamamoto, 2005).
A commonly accepted framework of nuclear receptor function is that ligand binding to the receptor induces the formation of a MTC that includes coregulator proteins and hormone response elements on the DNA. However, this framework falls short in explaining how nuclear receptor binding to similar hormone response elements is able to differentially regulate individual genes in a given cell environment (Bain et al., 2014). Currently accepted theoretical models of ligand-induced nuclear receptormediated gene expression predict that the amount of gene expression depends on the amount of ligand-activated receptor species present in the system (Kenakin, 2004;Chow et al., 2011). Therefore, the properties of the ligand dose-gene response curve provide a quantitative means to investigate gene expression. However, the cellular response to GR stimulation is not always a simple binary response to GC binding ) and numerous factors contribute to GC action at each step of the GR signaling cascade from ligand binding to its end-point result, the induction of gene transcription. Importantly, these factors appear to be tissue-, cell type-, and gene promoter-specific (Rogatsky et al., 2003;Bolton et al., 2007;Oakley and Cidlowski, 2011;Gertz et al., 2013).
The intermediate steps between ligand binding to the GR and the resulting gene promoter-specific induction include ligand binding and interaction with coregulator proteins that promote (coactivators) or inhibit (corepressors) gene expression and DNA binding (Meijer, 2002). Ligand binding per se is an important factor, as different ligands can induce specific GR interaction patterns with coregulators (Coghlan et al., 2003;Zalachoras et al., 2013). Further, specific coregulator binding to the activated GR plays a role in gene-specific induction (Lachize et al., 2009). Coregulators bind to the GR and modulate its transcriptional activity by modifying DNA structure. In particular, recruitment of transcriptional coactivators may destabilize chromatin by specific mechanisms including histone acetylation and contacts with the basal transcriptional machinery. In contrast, the recruitment of corepressors may stabilize chromatin by targeting histone deacetylases (Collingwood et al., 1999). Similarly, recent studies have demonstrated that coregulators can behave as allosteric modulators of the GR, affecting ligand interactions (Pfaff and Fletterick, 2010). These observations highlight the complexity of the molecular interactions taking place during the formation of a MTC.
Here, we aim to develop a model of GR interaction with coregulators and other components of a MTC that could describe and interpret differential ligand-specific transcriptional regulation of individual target genes. Our model is conceptually based on the cubic ternary complex receptor-occupancy model (Weiss et al., 1996), which has been extensively used before to describe the pharmacological behavior of multifactor complexes containing receptor proteins (Monczor et al., 2003;Fitzsimons et al., 2004;Tubio et al., 2010;Granja-Galeano et al., 2017). To experimentally test the model, we focused on gene transactivation, where the GR interacts with cofactors and DNA rather than acting by tethering mechanisms based on protein:protein interactions that influence the activity of other transcription factors without directly contacting the DNA (Newton and Holden, 2007). We used a set of previously identified GR-responsive genes in A549 human lung adenocarcinoma cells (Wang et al., 2004) and a MARCoNI. This array contained 54 coregulator-derived peptides representing nuclear receptor (NR)-boxes, or LXXLL motifs (Heery et al., 1997;Darimont et al., 1998), that interact with the activation function domain 2 (AF2) within the ligand-binding domain LBD of NRs. This assay can be used as a sensor for receptor conformation and activity status and allows the characterization of nuclear receptor binding to coregulators (Hur et al., 2004;Moore et al., 2004;Wu et al., 2004;Estébanez-Perpiñá et al., 2005;Teichert et al., 2009).
We found significant GR ligand-dependent differences in the relative efficacy and potency of induction of three GR-responsive genes in A549 cells and in the binding of coregulators to GR. Based on the behavior of GR ligands in A549 cells and in vitro coregulator recruitment, we developed a model of transcriptional regulation by the GR including ligand binding and interaction with coregulators that interprets gene-induction potency observed in living cells. The model includes a unique parameter δ describing the allosteric interaction between the components of the multifactorial complex. Using the nonindependent action indicated by the model, we were able to interpret gene-specific transcriptional inhibition by partial agonists in the presence of a full agonist of the GR in A549 cells. This model is supported by previous observations on GR-mediated gene-specific transcriptional activation and could be used to understand and interpret the pharmacological action of ligands that selectively modulate GR-dependent transcriptional activity.

Cell Culture and Treatments
A549 human lung carcinoma cells were obtained from the American Type Culture Collection (ATCC: Manassas, VA, United States) and cultured in complete medium (DMEM with 14.5 g/L glucose, supplemented with 10% fetal bovine serum and 1% penicillin and streptomycin, all from Invitrogen). Cells were cultured at 37 • C in a humidified 5% CO 2 atmosphere. For cell passaging or plating, cells were first washed out with 1X phosphate buffered saline (Invitrogen) and then trypsinized using 1X trypsin-EDTA (10X 0.5% Trypsin, Invitrogen). For all ligand treatments, A549 cells (200,000 cells/well) were seeded into six-well plates (Corning International, NY, United States) and cultured for 24 h in complete medium and maintained in steroid-free medium prepared as complete medium but using charcoal-stripped serum (Adams et al., 2003) for 24 h before ligand or vehicle (ethanol) treatments. DEX, RU486, and CYP (Sigma) stock solutions (1 mM) were dissolved in 100% ethanol and kept at −20 • C and further diluted in steroid-free medium before usage.

RNA Isolation and cDNA Synthesis
Total cellular RNA was extracted using the TRIzol R reagent (Invitrogen) following the supplier's manual (Invitrogen). Total RNA was dissolved in RNAase free water, denatured for 5 min at 65 • C and RNA concentration was quantified by spectrophotometric OD260 measurement using the Bioanalyzer (Agilent Technologies, Palo Alto, CA, United States). RNA samples were stored at −80 • C until further use. 1 µg of total RNA was used for cDNA synthesis. In order to remove genomic DNA carry-over, RNA samples were treated with 1.5 u of DNAase I (Invitrogen) for 15 min at 25 • C. DNAase I treated samples were then incubated at 65 • C for 10 min following addition of 25 nM of EDTA (Invitrogen). Finally, they were reverse transcribed using the iSCRIPT TM cDNA Synthesis Kit according to the manufacturer's instructions (Bio-Rad). From each DNAase I treated RNA sample, a non-reverse transcribed (-RT) sample was similarly generated (reverse transcriptase was replaced with water). cDNA as well as -RT samples were kept at −20 • C.

Quantitative Polymerase Chain Reaction (qPCR)
Forward and reverse primer pairs against reference (b-actin) and GR were generated using the primer3 Input on line software 1 and designs were based on publicly available human mRNA sequences. Primers were designed to have approximately 50% G/C content and to generate 150-250 bp amplicons. Primer pair specificity against target sequence was checked in the NCBI GenBank database using BLAST 2 . The sequences of the primers used to detect glucocorticoid-induced leucine zipper (GILZ), Solute Carrier 19A (SLC19A), and thrombomodulin (THBD) were provided by Dr. J. C. Wang and have been used before to detect gene expression in A549 cells (Wang et al., 2004). The sequences of the primers used to detect secretory leukocyte protease inhibitor (SLPI) and nuclear receptor subfamily 0 group B member 1 (NR0B1) were as follows: SLPI forward 5 -TCAAATGCCTGGATCCTGTTGA-3 ; SLPI reverse 5 -GCATCAAACATTGGCCATAAGTC-3 ; NR0B1 forward 5 -TGCTCTTTAACCCGGACGTG-3 ; NR0B1 reverse 5 -GCGTCATCCTGGTGTGTTCA-3 . In all cases, primers were supplied by Isogen Life Sciences (Netherlands) and dissolved in water according to the supplier's instructions and kept at −20 • C until use. qPCR monitoring and analysis was performed using the LightCycler R Carousel-Based Detection System 2.0 (Roche). PCR reactions were performed in a total volume of 10 µl containing 2 µl of LightCycler R FastStart DNA Master PLUS SYBR Green I master mix (Roche), 2 µl undiluted cDNA and 1 µl of each forward (5 pmol/µl) and reverse primer (5 pmol/µl). Every PCR reaction mix was filled in the LightCycler glass capillaries which were subsequently closed and centrifuged using the LC Carousel Centrifuge 2.0 (Roche). Cycling conditions were a single preincubation step at 95 • C for 10 min followed by 45 cycles of 10 s at 95 • C, 10 s at 60 • C and 10 s at 72 • C. To verify that the primer pairs used yielded single PCR products, a dissociation protocol was added after thermocycling, determining dissociation of the PCR products from 65 to 95 • C for 15 s. Finally, a cooling step was set for 20 s at 40 • C. 1 http://primer3.ut.ee/ 2 http://www.ncbi.nlm.nih.gov/BLAST To estimate the efficiency of the amplification reaction, serial half logarithm unit dilutions of cDNA from the A549 cells were used and standard curves were generated. The linear slope of the standard curve for each primer pair was estimated using GraphPad Prism 4 software and the efficiency was calculated based on the following equation (1).
Additionally, the -RT samples and a water-template were included in the analysis to confirm the absence of any residual DNA or contamination. All cDNA samples were analyzed in triplicates. Finally, the following equation (2)

Cell Transfection With siRNA
Transfections with siRNAs against the GR were performed as described before (Fitzsimons et al., 2013). Briefly, a total of 100 pmol of a GR targeting (Hs_NR3C1_6_HP validated siRNA, Qiagen and siGENOME NR0B1 siRNA, Dharmacon) or a nontargeting control siRNA (AllStrars Neg. siRNA AF 546, Qiagen or siGENOME Non-Targeting Control siRNAs #1, Dharmacon) were transfected into A549 cells using the Nucleofector I (Lonza) and the Cell Line Nucleofector R Kit T for the A549 cell line, according to the supplier's instructions (Lonza) using the U-29 Nucleofector program. The control siRNA was tagged with a red fluorophore, to monitor transfection efficacy, which was always higher than 90%. The cell medium was refreshed 1 day after transfection and all ligand treatments were performed 3 days after transfection.

Western Blotting
Total cellular proteins were extracted from A549 cells with fresh RIPA lysis buffer (200 µl/well) on ice, transferred into 1.5 ml eppendorfs, mixed and kept at −20 • C until use. Total protein concentration was measured using the BCA TM Protein Assay Kit according to the supplier's guidelines (Pierce). 3X of sample buffer was added to 10 µg of total protein, denatured for 5 min at 95 • C, spinned shortly and then loaded onto 10% SDSpolyacrylamide gels. Samples were run through stacking gel at 100 V for 10 min and separating gel at 200 V for approximately 1 h. Protein transfer onto methanol-activated Immobilon TM -P SQ Transfer membranes (Millipore) was performed overnight at 4 • C at 125 mA.
For the detection of the GR protein levels, the blots were incubated in blocking buffer, consisted of 5% low fat milk powder in TBST solution for 1 h at RT (10 ml/membrane) and subsequently with a primary antibody against the GR (H-300 rabbit polyclonal IgG, Santa Cruz), or against NR0B1 [Anti-NR0B1/Dax1 antibody (EP13786) -N-terminal (ab196649), abcam], or against alpha-tubulin antibody (clone DM1A, Sigma), or against GAPDH antibody (Santa Cruz). Primary antibodies were added in blocking buffer (1:2000 and 1:1000 dilution, respectively) for 1 h at RT (5 ml/membrane). Following 3X washing with TBST, the blots were probed with species-specific horseradish peroxidase-conjugated secondary antibodies (Santa Cruz) in blocking buffer (1:5000 dilution) for another 1 h at RT (10 ml/membrane) and finally washed 5X with TBST. For all incubations and washes rolling shakers were used. For luminescent signal detection, membranes were incubated with 10 ml of luminol solution, supplemented with 100 µl of enhancer solution and 3.1 µl of 30% H 2 O 2 for approximately 1 min at RT in the dark. Following film exposure, development and fixation, GR protein levels among samples were quantified relatively to a-tubulin signal using the ImageJ software 3 (Rasband, W.S., ImageJ, United States National Institutes of Health, Bethesda, MD, United States, 1997-2012).

Peptide Interaction Profiling
Interactions between the GR-LBD and coregulator NR-box peptides were determined using a MARCoNI assay (PamChip no. 88011; PamGene International) as described before (Koppen et al., 2009;Zalachoras et al., 2013). Each array was incubated with a reaction mixture of 1 nM Purified Glucocorticoid Receptor Recombinant Human Protein, Ligand Binding Domain, (Thermo Fisher Scientific, cat # A15668), ALEXA488-conjugated anti-GST antibody and buffer F (PV4689, A-11131, and PV4547; Invitrogen). For ligand induced peptide interaction profiling experiments 1 µM DEX, RU486, CYP or solvent (2% DMSO in water) were added. Incubation was performed at 20 • C in a PamStation96 (PamGene International). GR binding to each peptide on the array, reflected by fluorescent signal, was quantified by image analysis using BioNavigator software (PamGene International).

Nuclear Translocation Assay
GR translocation to the nucleus was studied using a YFP-GR construct kindly supplied by Dr. Cidlowski (National Institute of Environmental Health Sciences, National Institutes of Health) as previously described, with some modifications (Fitzsimons et al., 2008). Briefly, the previously described protocol was scaled down to a 96-wells plate format and semi-automated. 6000 A549 cells/well were plated 24 h prior to transfection. The cells were transfected with YFP-GR plasmid using a Nuclefector I (Lonza) as described before. Complete medium was refreshed 24 h after transfection and 48-h after transfection cells were incubated for 6 h in steroid free medium. After this procedure, the GR localized to the cytosol in all cells as described before (Fitzsimons et al., 2008). The test compounds or vehicle (ethanol) were manually dispensed into the corresponding wells and cells were incubated with the compounds for 30 min. Subsequently, cells were fixed with 80% acetone in water and stained with Hoechst 3342 (1:10.000) for nuclear staining. All experiments were performed in triplicates. Three non-overlapping images were taken from each well by using a Zeiss Axiovert 200/200M inverted microscope and 10× magnification. DAPI and FITC filters with excitation wavelength 409 and 487 nm were used to excite Hoechst (blue emission) and YFP (green emission), respectively. All images were collected with the same settings in Microsoft Window's BMP format. For post-acquisition image analysis Images were opened in ImageJ software and image backgrounds were subtracted by using the built in subtract background command with Rolling ball radius parameter set at 50. Subsequently, the images were opened with the CellProfiler software 4 to automatically identify the nuclear and cytoplasm compartments of cells and the green fluorescence intensity in each compartment (Carpenter et al., 2006).

Statistical Analysis
For comparisons between groups a two-tailed Student's t-test was applied using GraphPad's Prism 5 Software. For multigroup comparisons, a one-way ANOVA test with a Tukey's post test was performed using the same software package. Dose response curves were fit to a sigmoidal (four-parameter logistic) curve using GraphPad's Prism 5.

Dexamethasone (DEX) and RU486 Induce the Expression of GR Responsive Genes With Different Pharmacological Parameters
To model ligand-specific effects on GR-mediated gene transcription, we used three GR ligands with different pharmacological characteristics: DEX is a well-characterized GR agonist, while RU486 is usually described as a partial agonist (often used as antagonist), and CYP as a passive antagonist (Rousseau and Baxter, 1979;Honer et al., 2003;Meijer et al., 2005;Matthews et al., 2009). We studied the effects of increasing concentrations of DEX on the expression of three GR responsive genes in A549 cells. Glucocorticoid Induced Leucine Zipper (GILZ) is an important mediator of the anti-inflammatory effects of GCs (Ronchetti et al., 2015), THBD is an endothelial cell surface glycoprotein that controls thrombosis by downgrading thrombin-mediated fibrin generation and promoting protein C activation (Loghmani and Conway, 2018) and SLC19A2 is a thiamine transporter associated with the thiamine-responsive megaloblastic anemia syndrome (TRMA) (Aoyagi and Archer, 2011). In this cellular system, DEX dose-dependently induced the expression of the three genes tested, albeit with different pharmacological parameters. We found significant differences in calculated effective concentration 50 (EC 50 ) and maximal response (R max ) values. Specifically, DEX was significantly more potent in inducing GILZ than THBD or SLC19A2, and significantly less efficacious in inducing SLC19A2 than THBD or GILZ ( Figure 1A and Table 1).
These differences are difficult to explain from the standpoint of theoretical models of ligand-induced gene expression because the ligand, the amount of receptors and all other components of the cellular environment were the same in all experiments. One possible explanation could be in the structure of the FIGURE 1 | DEX and RU486, but not CYP, induce a dose-dependent and gene-specific response in A549 cells. (A) Gene expression response of three responsive genes, GILZ, SLC19A2, and THBD, to increasing concentrations of DEX measured by qRTPCR. Results are mean ± SEM of three independent experiments performed in triplicates. Fitted parameters are detailed in Table 1. (B) Binding affinity predictions corresponding to GREs present in proximal 5 UTR regions to transcription start sites (TSS) corresponding to GILZ, SLC19A2, and THBD. (Continued)

FIGURE 1 | Continued
Results are expressed as GR binding scores (GBS), as previously described (Datson et al., 2011). (C) Gene expression response of three responsive genes, GILZ, SLC19A2, and THBD, to increasing concentrations of RU486 measured by qRTPCR. Results are mean ± SEM of three independent experiments performed in triplicates. Fitted parameters are detailed in Table 2. (D) Gene expression response of three responsive genes, GILZ, SLC19A2, and THBD, to increasing concentrations of CYP measured by qRTPCR. Results are mean ± SEM of three independent experiments performed in triplicates.  gene promoter and its GRE composition (Wang et al., 2004;So et al., 2008), suggesting that response magnitude depends on the strength of GR binding to GREs. However, numerous previous studies argue to the contrary (So et al., 2008;Pfaff and Fletterick, 2010;Dougherty et al., 2012;Chen et al., 2013), and previously published chromosome immunoprecipitation (ChIP) data analyzing these promoters do not predict all the differences in sensitivity to the GR agonist DEX (Wang et al., 2004). The promoter sequence data predicted GILZ to be strongly bound by the occupied GR with two high GR binding score (GBS) sites in its promoter, followed in order of predicted responsiveness by SLC19A, with one high GBS sites, and THBD with one low GBS site (Figure 1B), while DEX dose response curves showed a different order in both response efficacy (THBD∼GILZ>>SLC; Figure 1A) and potency (GILZ>>THBD∼SLC; Figure 1A). Similarly, RU486 induced significant changes in GILZ, SLC19A2 and THBD, expression in A549 cells, albeit only at high doses ( Figure 1C and Table 2). Importantly, the observed order of gene-transcription potency (SLC19A∼THBD>GILZ; Table 2) and efficacy (GILZ>SLC19A∼THBD; Table 2) induced by RU486 were different than those induced by DEX (Table 1). Noteworthy, although CYP alone did not have any detectable effect on the expression of the three genes analyzed (Figure 1D), DEX, RU486 and CYP induced GR translocation from the cytosol to the nucleus, indicating that the three ligands promote active changes on GR behavior and therefore cannot be ascribed as passive (or inactive) antagonists (Table 3). To understand the differential, ligand-and gene-specific, biological behavior observed in A549 cells in more detail, we aimed to develop an MTC model that could integrate ligandreceptor-DNA-coregulator interactions.

Development of a MTC Model for Ligand-Receptor-DNA-Coregulator Interactions
An alternative hypothesis to explain the gene-specific response to DEX observed in A549 cells involves the differential recruitment of specific coregulators to the MTC active at each gene promoter, which could induce gene-specific expression. To model this situation, we postulated a theoretical equilibrium model of receptor action that explicitly includes the receptor (R), the ligand (L), the coregulator(s) (C), and the DNA (D) and four parameters that govern receptor species equilibria: α represents the effect of ligand binding on the binding of coregulator, β describes the effect of coregulator binding on the binding of receptor to DNA and γ the effect of ligand binding on the binding of receptor to DNA. In turn, δ represents the extent to which the joint effect of any two of ligand binding, coregulator binding or DNA binding varies conditional on the level of the third. A detailed description of all these parameters is presented as Supplementary Material and schematically in Figure 2. Our model assumes that the GR can spontaneously couple to the coregulators or the DNA even in the absence of the ligand (Power et al., 1991;Mani et al., 1994). Thus, the model is fully described by three basic equilibrium constants that account for ligand binding, coregulator coupling and DNA binding and four parameters that illustrate the interaction effect between them (Figure 2A and Supplementary Information S1).
Assuming that the MTC formed by the receptor, the ligand, the coregulator(s) and the DNA is responsible for the final induction of gene transcription, dose-response curves can be simulated with equations describing how the relative composition of the different components of the MTC affects ligand-dependent gene induction (Figures 2B-E). The model predicted a series of system characteristics that could be experimentally validated. According to the model, the concentration-response curve to a ligand can be modified adjusting the values assumed for each component of the MTC (Supplementary Information S1). Indeed, a first validation of the model using siRNAs to reduce GR expression in A459 cells resulted, as indicated by our model, in significant changes in R max but not EC 50 (Figure 3). Transfection of A549 cells with increasing siRNA concentrations resulted in concomitant decreases in GR expression (Figures 3A,B) and in GILZ's R max to DEX, without affecting its EC 50 ( Figure 3C and Table 4).
Importantly, when changes in the α, β, or γ parameters are simulated, the results predict changes in either R max or EC 50 , but never both simultaneously (Figures 2B-D and Supplementary  Information S1). This prediction leaves out the possibility that differences in coregulator recruitment or in DNA binding alone could explain the simultaneous change in EC 50 and R max observed in DEX concentration-response curves in A549 cells (Figure 1). However, when the δ parameter is taken into consideration, simultaneous variation in both EC 50 and R max can be explained theoretically using our model ( Figure 2E). δ describes how two binding events impact on a third within the MTC proposed by the model, reflecting interaction between all the components of the complex. Therefore, our model suggests that non-independent binding events between the different components of the multifactor complex could be responsible for the simultaneous differences in EC 50 and R max to DEX observed in A549 cells ( Figure 1A and Table 1).
The sequence of binding events leading to transcriptional activation by the GR has been intensively studied but not fully established yet (Rogatsky et al., 2003;John et al., 2009;Dougherty et al., 2012;Chen et al., 2013). It is considered to be a dynamic process and the most parsimonious hypothesis is that transcriptional activation by the GR, and NRs in general, involves multiple factors that act in both a sequential and combinatorial manner to reorganize chromatin templates (Pollard and Peterson, 1998;Glass and Rosenfeld, 2000;Voss et al., 2006;Stavreva et al., 2012). Moreover, the temporal order of the events leading to the formation and composition of the MTC that leads to transcription activation can take place in a gene-and cell-specific manner (Gronemeyer et al., 2004). Potential dissimilarities in the structure of promoters for GILZ, SLC19A, and THBD genes do not explain the differential expression patterns obtained with DEX and RU486 ( Figure 1A vs. Figure 1C). This divergence can be explained by our MTC model, considering that each ligand induces a specific pattern of coregulator binding to GR forming the ternary complex LRC, which in turn may display differential affinity for DNA. Therefore, using a DNA-free system represented by the MARCoNI peptide array, we focused on understanding how the interactions between ligand and GR could cooperate to modulate coregulator binding.

Ligand-Independent and Ligand-Specific Receptor-Coregulator Binding Events
Our model includes the existence of ligand-specific parameters governing multifactor complex formation Atucha et al., 2015) that can be tested experimentally in DNA-free conditions (Supplementary Information S1). Ligand affinity constant Ka (specific for each ligand), coregulator affinity constant Kc (specific for each individual coregulator NR-box), and the parameter α (characteristic of each ligand/coregulator pair) indicate how ligand and coregulator affects each other's binding to the receptor and are key components of the model when only receptor, ligand, and Ligand-independent and -dependent coregulator binding profiles were measured using a MARCoNI peptide array. This array contained 53 coregulator-derived peptides representing a wide range of coregulator NR-boxes known to interact with FIGURE 3 | GR knockdown affects GILZ response to DEX. (A) Western blot depicting the effect of varying concentrations of a previously described specific siRNA targeting the GR (Fitzsimons et al., 2008) on GR protein levels in A549 cell lysates. The image shown is representative of five independent blots. (B) Quantification of the effect of varying concentrations of the specific siRNA used in (A) on GR expression. Results are expressed mean ± SEM of five independent blots. Statistically significant changes were identified using Student's t-test. * p < 0.05; * * p < 0.01. (C) DEX dose-dependent effect on GILZ expression at decreasing GR expression levels induced by siRNA-induced GR knockdown. The calculated EC 50 values were not affected by GR knockdown, while GILZ maximal response to DEX was significantly attenuated by increasing GR knockdown ( Table 4). DEX-induced GILZ expression was detectable even at maximal GR knockdown (100 pmol siRNA), indicating GILZ induction is robust even at low levels of GR expression.
the GR and other nuclear receptors, which sequence details, Gene Name and UniProt Knowledge Base accession numbers are shown in Supplementary Table S2. First, NR box peptides were titrated against the GR LBD in the absence of ligand using a customized MARCoNI array (Figure 4). A concentration series of each NR box peptide was immobilized on the array and incubated with 3 nM of the GR LBD and binding isotherms were calculated from these data. We observed a significant basal  binding of the GR LBD to several NR-box peptides in the absence of any ligand ( Figure 4A). This basal coregulator recruitment in the absence of ligand is predicted by the model (RC species in Figure 2) and is in agreement with similar observations done using fluorescence polarization assays (Pfaff and Fletterick, 2010).
Twelve NR box peptides showed positive ligand-independent binding to the GR LBD, although with variable binding profiles (Supplementary Table S3). Ten of these NR box peptides showed lower affinity binding profiles, while two NR box peptides showed higher affinity binding. These higher affinity binding peptides corresponded to the coactivator PRGC1 and the corepressor NRIP1 ( Figure 4B and Supplementary Table S3). Secondly, we used a MARCoNI array in which 1 mM of each NR box peptide was immobilized and incubated with 3 nM of the apo GR LBD in the presence of a receptor-saturating concentration of three selected GR ligands or a vehicle control, to study the effect of GR ligands on the basal binding of GR LBD to NR box peptides. We analyzed the coregulator binding profile induced by DEX, RU486, and CYP. These three ligands induced characteristic coregulator binding profiles, in some cases favoring and, in some others, disfavoring the basal binding of the apo GR LBD to NR box peptides observed in the absence of ligand ( Figure 5). Notably, the NR box peptide binding profile induced by RU486 resembled for some peptides the effect induced by DEX, although RU486's effects were substantially weaker in many cases, reflecting its partial agonist activity (Figures 5A,B).
Interestingly, GR LBD binding to NR box peptides from the coactivators NCOA2, NCOA3 and NR0B1 was strongly favored by DEX, while it was disfavored by RU486 (Figures 5A,B), suggesting that this differential effect on coregulator binding is involved in the differential pharmacological effects of the two ligands. Supporting this preliminary conclusion, we did not find any example of the opposite pattern in our dataset, this is, binding to NR box peptides that were disfavored by DEX but promoted by RU486. In contrast, the NR box peptide binding profile induced by CYP was more divergent from that induced by DEX (Figures 5B,C). We did not find any example in our dataset of NR box peptides for which binding to GR was favored by CYP but disfavored by DEX ( Figure 5B). Interestingly, from the 10 NR box peptides whose binding was disfavored by CYP, 7 were favored by DEX (Figure 5B), suggesting that these differences in coregulator binding may be crucial to understand the pharmacological differences between DEX and CYP.
In summary, the pharmacological behaviors observed for the three GR ligands used in experiments in A549 cells (Figure 1) can be interpreted using our MTC model considering the interactions between receptor-ligand-coregulator-DNA binding events. Our model suggests that there exist an allosteric phenomenon involving the joint effect of the three GR partners within the MTC (L; C and D), reflected by the δ factor, i.e., a specific ligand induces the recruitment of a specific set of coregulators, that differentially affects the expression of a particular gene. In fact, our experimental observations in vitro with the MARCoNI, support the induction of ligand-specific binding profiles between the GR LBD and NR box peptides Atucha et al., 2015) that can explain the differential gene-specific behaviors of GC ligands. The model indicates that the distinctive ligand-specific binding of GR to coregulators would differentially impact on the expression of specific genes when cells are coincubated with DEX and RU486 or CYP. Indeed, Figure 6 shows that RU486 can block DEX-induced GILZ expression without affecting SLC and THBD expression significantly, while CYP specifically blocked DEX-induced SLC and THBD expression without affecting GILZ maximal expression levels (Figure 6).
According to the MTC, each ligand may induce a specific response by recruiting a defined array of cofactors. Among the coregulators that were differentially bound by GR ligands, we chose NR0B1 as a proof-of-concept experimental validation. NR0B1 displayed a maximum difference between DEX-favored and RU486-disfavored binding to the GR LBD in vitro (Figure 5). SLPI, a protease inhibitor expressed by cells at mucosal surfaces, is stimulated by IL-1β. IL-1 β-induced SLPI expression is increased by DEX and inhibited by RU486 in A549 cells (Ito et al., 2001). We found that SLPI expression is strongly induced by DEX but not by RU486 in non-transfected A549 cells (Figure 7A). While NR0B1 downregulation using specific siRNAs (Figures 7B-D) significantly diminished the maximal DEX-induced expression of SLPI ( Figure 7E) and simultaneously increased DEX's EC 50 for SLPI induction, it did not affect RU486induced SLPI expression (Figure 7F), demonstrating its specific role in DEX-induced SLPI expression in A549 cells.

DISCUSSION
In this work, we have developed a mathematical model describing the interactions between the GR and other components of an MTC that control the transcription of GR-target genes. This model is based on basic equilibrium equations governed by three constants and four parameters accounting for ligand, coregulator and DNA binding events and accommodates plausible nonindependent effects among them. Here, we show proof-of concept for the model using A459 cells, a well characterized system to study GR-mediated activation of responsive genes , and three extensively characterized GR ligands; DEX, RU486, and CYP. Future research should be aimed to extensively validate the model using multiple cell lines, ligands, coregulators, and responsive genes.
Factors that regulate gene transcription are assembled in multicomponent complexes by combinatorial interactions (Britten and Davidson, 1969;Gierer, 1974;Yamamoto et al., 1998). In this context, the GR provides a well-characterized and physiologically relevant study system in which the effects of ligand dose and chemistry, treatment duration and kinetics, interaction with coregulatory factors and DNA binding site sequence and structure have all been intensively studied Collingwood et al., 1999;Rogatsky et al., 2003;John et al., 2009;Reddy et al., 2009;Stavreva et al., 2009). Interestingly, many of these factors induce allosteric changes on the GR (Meijsing et al., 2009;Pfaff and Fletterick, 2010;Wang et al., 2012;Watson et al., 2013), suggesting that allosteric conformational changes induced by components of the MTC are crucial to understanding GR-mediated gene transcription. In fact, a model based on such a "conformational ensemble" has been proposed to explain ligand-induced switch between ER-α-mediated genomic and non-genomic effects (Norman et al., 2004).
One of the significant challenges in modeling GR-mediated gene transcription is to understand how a multiple step FIGURE 5 | DEX, RU486 and CYP induce specific GR LBD-to-NR box peptide binding profiles. (A) Characteristic binding profiles induced by 1 × 10 −7 DEX (red), RU486 (blue) or CYP (green) obtained using the quantitative in vitro assay, MARCoNI. Modulation Index (MI) > 0 suggests ligand-favored binding, while MI < 0 suggests ligand-disfavored binding of a peptide compared to DMSO. (B) Heatmap depiction of details of ligand-induced binding of coregulator peptides using MARCoNI. (C) Venn diagrams showing the number of peptides whose binding was favored (left), unfavored (center) or unchanged (right) by GR ligands. In all cases, statistically significant changes relative to DMSO were identified by Student's t-test. * p < 0.05, * * p < 0.01 or * * * p < 0.001. reaction sequence can result in dose response curves that follow a Michaelis-Menten function as we show here in A549 cells. Explaining this experimental observation requires careful consideration of parameters such as EC 50 and R max . Consequently, a theoretical framework able to predict changes in shape, position and R max of ligand dosedependent response has been developed (Ong et al., 2010). This model suggests that full dose-response curves need to be considered to achieve complete mechanistic insight into GR-mediated gene transcription. The model we describe here is in agreement with this concept and is conceptually based on the Cubic Ternary Complex Receptor-Occupancy Model (CTC) originally proposed by Weiss et al. (1996). The CTC has been successfully used to describe interactions between receptors and their interacting partners, and to predict doseresponse curves for many ligand/receptor pairs (Weiss et al., 1996). Indeed, due to its theoretical completeness, the application of the CTC framework has allowed us to experimentally characterize low-abundance receptor states predicted by the model in other systems (Monczor et al., 2003;Fitzsimons et al., 2004;Tubio et al., 2010). Following the law of mass action, our model describes the interactions between a ligand, a receptor and accessory molecules, including coregulators and DNA sequences resulting in eight receptor species that coexist at equilibrium; the free receptor (R) and receptor species bound to the ligand (LR), to the accessory coregulators or DNA (RC and RD, respectively), and bound to any two or three factors (LRC, LRD, RCD, LRCD) (Figure 2 and Supplementary Information S1). There are several mechanisms by which GR modulates gene expression. We decided to experimentally validate our model measuring GR mediated gene transactivation because the model's assumptions reflect better transactivation mechanisms, since GC-induced gene transactivation is mediated by the direct interaction of GR with DNA and cofactors, while gene repression is largely based on by protein:protein interactions (Ratman et al., 2013). FIGURE 7 | NR0B1 knockdown affects SLPI response to DEX but not to RU486. (A) DEX but not RU486 induces SLPI expression in A549 cells, as expected, coincubation with RU486 blocked DEX-mediated effects, * * p < 0.01. (B) Effect of cell transfection with siRNA targeting NR0B1 on mRNA and (C,D) protein levels in A549 cells measured by qRTPCR and Western-blot, respectively. The expression of the glyceraldehyde 3-phosphate dehydrogenase (GADPH) was used as internal control and normalization in (C,D), * * * p < 0.001. (E) Transcriptional response of SLPI to increasing concentrations of DEX measured by qRTPCR in cell transfected with specific siRNAs against NR0B1 or a scramble siRNA control. The fitted parameters are detailed in Table 5. (F) Transcriptional response of SLPI to RU486 measured by qRTPCR. Results are expressed as mean ± SEM of three independent experiments performed in triplicates.
The experimental validation of predictions of this model in A549 cells, a well-characterized system to study GRdependent gene expression , and in vitro using a system of real time quantification of GR-coregulator interactions, suggests that the model can help to interpret specific ligand dose-dependent effects on gene expression, and could be used to explain the pharmacological characteristics (e.g., the simultaneous variation in relative potency and efficacy) of ligand-specific and gene-specific transcriptional regulation by the GR that are difficult to explain using alternative models based only on ligand potency. One of the distinctive features of our model is that receptors may spontaneously bind any partner without binding the ligand. This feature is particularly interesting for the exploration of pharmacological properties of GR ligands that seem to behave as neutral antagonists in vivo and as partial agonists in vitro such as RU486 and the non-steroidal CP-472555 (Weigel and Zhang, 1998). We verified the gene-specific behavior indicated by our model for RU486 and DEX in A549 cells. Interestingly, SLC19A2 and THBD showed 50% gene expression response (EC 50 ) values for DEX in the low nM range (Table 1), in agreement with previous reports (Reddy et al., 2009). However, GILZ robustly responded to significantly lower concentrations of DEX (Table 1). At least one other gene, PER1, with similar high sensitivity to DEX has been shown to have similar characteristics and its increased responsiveness could not be explained by GRE composition or Pol II occupancy in its promoters (Reddy et al., 2009), suggesting that differential binding to coregulators maybe involved. The order of gene induction potency and efficacy observed for the partial agonist RU486 was different from the one observed for the full agonist DEX, a crucial observation that cannot be explained using conventional pharmacological models based only on ligand potency. We have previously explored the possible physiological relevance of downregulating the GR using RNA interference in the mouse brain, showing that this experimental approach is technically feasible in vivo (Fitzsimons et al., 2013). The results presented here using siRNAs to downregulate GR expression in A549 cells indicate that GR downregulation may result in a decrease of the maximal response to GR without affecting its EC 50 , which we exemplify using DEX and one GR responsive gene, GILZ. If these results could be generalized they may be relevant to further understand the physiopathological relevance of the expression of dominantnegative GR isoforms, such GR-beta (Yudt et al., 2003;Lewis-Tuffin and Cidlowski, 2006), which has been associated with decreased GC responsiveness and susceptibility to develop autoimmune diseases (Tait et al., 2008). Using the in vitro system provided by the MARCoNI chip, we observe active binding of the GR-LBD to both coactivators and corepressors in the absence of ligand. This observation is compatible with similar ones made before using the GR-LBD (Pfaff and Fletterick, 2010). Possibly, full-length GR species that spontaneously bind full-length coregulators may exist at very low levels, perhaps undetectable or destabilized in most cell-based assays and thus, their existence in vivo will require further experimental validation. Furthermore, the GR amino terminal end, not included in the recombinant peptide used in our MARCoNI assay, also contains binding sites for several coregulators (Godowski et al., 1987;Eickelberg et al., 1999;Siriani et al., 2003). Importantly, the notion that the GR may exist in native conformation ensembles capable to assuming active and inactive behaviors is compatible with the conformation ensemble model proposed to explain cellular behaviors mediated by the ER-alpha (Norman et al., 2004). Several studies suggest that the GR is refractory to ligandindependent activation (reviewed in Weigel and Zhang, 1998). Others have suggested that the GR can be activated in the absence of hormone (Cenni and Picard, 1999;Eickelberg et al., 1999). Mutagenesis studies demonstrated that single amino acid mutations and phosphorylation events can render the GR constitutively active in the absence of hormone (Godowski et al., 1987), reinforcing the idea that GR conformational states allowing ligand-independent activity exist, albeit in low abundance with respect to the inactive forms. Supporting this hypothesis, some studies have demonstrated the ability of the GR to regulate gene expression through non-hormonebinding forms of the receptor in overexpression systems (Siriani et al., 2003). Therefore, the CTC framework provides a theoretical environment to understand and interpret GR behaviors that may be mediated by low abundance receptor states (Tubio et al., 2010).
Here, we introduce the application of the most parsimonious version of the CTC to understand GR-mediated cellular behaviors. The model (Figure 2) includes one GR, one ligand, one coregulator and the DNA. This is a representation including all the thermodynamic equilibria between the four partners involved in the MTC and results in a convenient cubic depiction of the eight receptor species discussed above. This representation cannot accommodate simultaneous interactions with more than one coregulator, which could be expected given the wellcharacterized complexity of the MTC (McKenna and O'Malley, 2002). However, the basic cubic structure of the CTC can be easily extended to include more than one accessory species, such as multiple coregulators. This is done by simply joining two cubes, differing only in their accessory species, by their accessory species-free face resulting in a model represented by a double cube. If more than one accessory species is to be added, the model cannot be visualized in three dimensions but can still be modeled and analyzed mathematically (Weiss et al., 1996). Similarly, the GR is the product of a single gene from which multiple transcriptional and translational isoforms are generated through alternative splicing and alternative translation initiation. Multiple studies have shown before that the specificity of GR signaling may arise, at least in part, from this molecular diversity, because some of these isoforms have differential ligand affinities and efficacy to induce the expression of target genes (Cain and Cidlowski, 2015). This complexity seems difficult to grasp, however, in terms of our model each GR isoform could be modeled using a different CTC with specific affinity constants as described in Figure 2. Although admittedly practically cumbersome, this procedure does not impose theoretical limitations for our model.
With respect to the validity of the model predictions presented in the manuscript, our in vitro experiments were done with the LBD of the human GR-alpha, which is the most abundant splice variant of the GR in most tissues (Pujols et al., 2002). For the cellular studies that we performed in A549 cells, the full-length GR mRNA expressed there may present N-terminal translation variants (Lu and Cidlowski, 2005), and these go undefined in the vast majority of studies. However, in our case this diversity in N-terminal translation variants may be irrelevant, as the coregulator interactions that we model concern the LBD (AF-2). Therefore, by including all possible thermodynamic equilibria between the species involved, our model considers the nonindependent effect of multiple factor binding events on GR ligand dose-response curves.
An important prediction of our model is that a basal GRcoregulator binding profile is differentially modulated in a ligandspecific manner. All the receptor species present in a ligand-free environment, i.e., R, RC, RCD and RD are induced to form LR, LRC, LRCD and LRD in the presence of a ligand and this induction is governed by the ligand-specific constant and parameters Ka and α, γ and δ (Figure 2). Our observations support the view that the changes induced by different ligands on a basal profile of GR binding to coregulators could play a significant role in cell type and tissue-specific ligand actions, such as those observed before with RU486 and other selective SGRMs Atucha et al., 2015). We performed a proof-of-concept validation by downregulation NR0B1 in A549 cells (Figure 7). Downregulation of NR0B1 levels resulted in a significant decrease in DEX-induced SLPI expression, without affecting SLPI's lack of response to RU486. These data indicate that gene-specific response to GR ligands can be, at least partially, regulated by coregulator abundance.
Our model predicts a number of complex interactions between different molecular species engaged in the MTC, however, we have not been able to validate all of them experimentally. Importantly, all the interactions in the MTC may be changed upon posttranslational modifications of the GR or the coregulator proteins, that may affect the affinity of the protein:protein interactions or the localization of the proteins involved. The interactions between DNA and other components of the MTC are particularly challenging because their full understanding would require characterization in living cells, where these interactions are most relevant. DNA has been proposed to act as an allosteric ligand of the GR in cell-based assays (Pfaff and Fletterick, 2010), we face technical limitations of the MARCoNI array to experimentally measure interactions with DNA. In view of these limitations, in the present work we focused on interactions between the GR and coregulators. There are other factors, such as variations in DNA binding motives, number of GREs, the role of chromatin structure and epigenetic factors, which may influence DNA binding of the GR or other components of the MTC. In our model, these are variations on the theme of specific sequence and structure of individual GREs -the parameters in the cubic model will depend on specific sequence-be it genomically or epigenetically determined. Thus, these variations are automatically incorporated -per GREdependent process, in the model. Our experimental validation incorporates cell line-specific effects of the GR assuming that they will, at least in part, result from differences in coregulator expression. In this respect, we used A549 cells in our primary validation of the GR-responsive genes used in this work because these were previously extensively characterized in this cell line, including their GRE composition and sequence, thus leaving out factors that would only introduce uncertainty in our validation steps such as variations in DNA binding motives, number of GREs (Wang et al., 2004) and epigenetic factors (Reul et al., 2009), mentioned before. In conclusion, our theoretical model based on the mathematical framework of the CTC is able to accurately interpret a variety of GR behaviors in different experimental setups, ranging from interaction with coregulator binding in vitro to differential effects on gene expression in A549 cells and may be used to characterize and interpret ligand-and gene-specific effects on transcriptional activity. Our observations may explain previous reports of RU486 and CYP having (partial) agonistic activities and justifies their classification as SGRMs, in the sense that they interact differentially with subsets of the GR functions induced by the full agonist (Gronemeyer et al., 2004), probably in a gene-specific manner as indicated by our model and the coregulator binding profiles they induce in vitro.

DATA AVAILABILITY
All datasets generated for this study are included in the manuscript and/or the Supplementary Files.

AUTHOR CONTRIBUTIONS
CF and FM developed the model, designed the study, and wrote the manuscript. AC, CZ, and RH performed the experiments. OM supervised the experiments and corrected the manuscript. All authors analyzed and interpreted the data and revised the manuscript critically for important intellectual content and approved the final version.