Quantitative modeling of STAT1 and STAT3 dynamics in IFNγ and IL-10 pathways uncovered robustness of anti-inflammatory STAT3 signaling

Abstract Information relay by signal transduction pathways occasionally involves sharing of functionally opposing transcription factors. For instance, IFNγ-STAT1 and IL-10-STAT3 elicit pro and anti-inflammatory cellular responses, respectively, but IFNγ mediated STAT3 and IL-10 mediated STAT1 activation is also observed. Here, through experiments at the cell population level, we studied the dynamics of STAT1 and STAT3(S/1/3) in responses to both IFNγ and IL-10 stimulation and subsequently trained a simplified mathematical model that quantitatively explained the S/1/3 signal-response at different doses of IFNγ and IL-10. Next, to understand the robustness of the canonical signaling axis in each pathway, we simulated a costimulation scenario (IL-10 and IFNγ applied simultaneously) which predicted STAT3 activation would remain IL-10 driven in presence of IFNγ; subsequent experiments validated the same. We next investigated how protein expression variability may plausibly influence the robustness of IL-10-STAT3 signaling at the level of individual cells. Simulating thousands of single cells and analyzing their responses to co-stimulation we could identify emergence of two new subpopulations; in one subpopulation co-stimulation dominantly activated STAT3 and suppressed STAT1 activation, and, vice versa in the other subpopulation. Analyzing the protein concentration from these reciprocal subpopulations we found the key proteins whose cell-specific expression could control S/1/3 responses in individual cells. Taken together, we present a quantitative model that captures the signaling dynamics of STAT1 and STAT3 in response to functionally opposing cues and through single cell simulations show how reciprocal responses could emerge at the level of individual cells.


Introduction
Information encoded in the dynamics of signalling pathways triggers a plethora of biological processes like cell growth, proliferation, apoptosis or developmental lineage scenario where the canonical IL-10 could inhibits the IFNγ mediated activation of STAT3. The single cell simulations suggest, in addition to the cells exhibiting responses similar to the population level experiments, two additional subpopulations with reciprocal STAT1 and STAT3 activation profiles could emerge as a function of cell-tocell variability. Finally, analyzing the distribution of protein concentrations from these reciprocal subpopulations we identified key proteins whose relative concentration would result in either a dominant STAT3 or STAT1 response in individual isogenic cells.
We further demonstrate that with targeted minimal perturbations S/1/3 responses can be tuned further, such that, occurrence of a desired subpopulation type can either be enriched or inhibited.

IFNγ
We took the peritoneal macrophages obtained from BALBc mice and stimulated them with increasing doses (0. 5  . STAT3 phosphorylation upon IFNγ stimulation exhibits a distinct profile: at M dose STAT3 is rapidly phosphorylated to a high amplitude and gets rapidly inhibited as well, but at L and H doses STAT3 phosphorylation is relatively negligible (Figure1C, 2 nd panel). Such bell-shaped dynamics of signalling intermediates is also observed elsewhere [30]. In response to IL-10, the STAT1 phosphorylation peaks at M dose and a slightly inhibited response can be observed at the H dose ( Figure 1D, 1  induced/activated as a function of applied signal [11,14,31,32]. SOCS1 being a commonly induced negative regulator in IFNγ signalling [19, 20, 22 - of SOCS3 induction was also captured which is a target gene downstream to both the stimulation types ( Fig 2E, 2 nd panel). SOCS3 in principle can also inhibit IFNγ signalling, but the relative inhibitory strength of SOCS3 is observed to be negligible compared to SOCS1 [33], hence w.r.t IFNγ stimulation we studied SOCS3 only as a target gene.

Quantitative modelling captures dynamics of STAT3 and STAT1 activation at different doses of IL-10 and IFNγ ligand
We used the observed dynamics of S/1/3 to calibrate a mathematical model comprising both IFNγ or IL-10 pathways.  Figure 2A schematically shows the minimal models of IFNγ pathway. The model has a simplified step of receptor activation where details of the interaction between IFNγ activation process [23]. As the ligand binds to the receptor the active receptor complex (IFN-LR) phosphorylates the transcription factors STAT1 and STAT3. Both the STATs undergo dimerization and forms transcriptionally active complexes [21] which in turn induces target genes such as SOCS1 and SOCS3. The time delay in transcriptional induction is captured with Hill functions [14]. Dephosphorylation of S/1/3 are assumed to be carried out by phosphatases such as SHP2 [20,36,37]  IL-10 pathway model Figure 2C shows the schematics of IL-10 receptor-mediated phosphorylation of STATs and the transcriptional induction of the SOCS. Similar to the IFNγ pathway, the explicit steps of IL-10 receptor1(IL-10R1) and receptor 2(IL-10R2) binding to JAK1, Tyk2 kinases leading to the formation of active signaling complex [38] is simplified to one step activation process [25]. STAT1/3 phosphorylation steps are explicitly modeled and dephosphorylation of STAT1/3 is assumed to be carried out by a phosphatase [36,37].
Negative regulation of IL-10 receptor through degradation is considered and we named the inhibitor as IL-10Ri; such negative regulators are observed to act by sequestering the IL-10 receptor 1(IL-10R1) that subsequently leads to ubiquitination and degradation [22]. Figure 2D shows the STAT1/3 dynamics in the model upon IL-10 stimulation with L, M, and H doses. As both the pathways are built as one model with common signaling intermediates, we calibrated both IL-10 and IFNγ models simultaneously to their respective datasets. Hence, during the model calibration, the common signaling intermediates and biochemical parameters of both the pathways were constrained to have a common value that fits the S/1/3 activation in response to both stimuli.

Mechanism controlling dose dependent kinetics of the STATs
The observed bell-shaped dose response of STAT3 in both IFNγ and IL-10 stimulation, as well as the graduating pattern in activation of STAT1 for both M and H dose of IFNγ . CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/425868 doi: bioRxiv preprint first posted online Sep. 24, 2018; stimulation, is captured simultoneously by our simplified model. To capture the signalresponse of all the training datasets simultoneously the model required stoichiometric inhibition of IFN-LR by SOCS1 and IL-10-LR by IL-10Ri, respectively [20][21][22][23][24][25]. SOCS1 is a commonly reported feedback regulator in the IFNγ pathway but our experiments show negligible SOCS1 induction upon stimulation (Fiure1E, top panel, dotted line) and the basally present SOCS1 acts as the stoichiometric inhibitor of IFNγ signaling. Negative regulation of signal at the receptor level in both the pathways and activation of a S/1/3 phosphatase in the IFNγ pathway [20-25, 39, 40], are suggested by the model as the key regulatory processes controlling the observed S/1/3 responses.

Prediction and validation of a co-stimulation scenario: STAT3 dynamics remain robustly IL-10 driven in presence of IFNγ stimuli
As SOCS1 is strongly induced upon IL-10 but not upon IFNγ stimulation, in a scenario when both these functionally opposing pathways are simultaneously activated(costimulation) an IL-10 driven negative regulation of the IFNγ signaling can be expected.
To quantitatively understand the consequence of IL-10 induced SOCS1 on the IFNγ pathway we next simulated both the pathways simultaneously. We chose the M doses of both the ligand types as both the canonical and cross-activated STATs are strongly activated in M dose( Figure 2B and 2D). Figure 3A depicts the co-stimulation scenario and the predictions for dynamics of STAT3 (Fig 3B), STAT1( Fig 3C) and SOCS1( Fig 3C) are shown(representative immunoblots in figure S2B and S2C). To achieve robust predictions we used 40 independently fitted models with comparable goodness of fit(see in methods for details). The models predict: IL-10 induced SOCS1 would impart additional inhibition of IFNγ signaling as the excess SOCS1 would inhibit the active IFN-LR, blocking its access to S/1/3. The simulation suggest, cells subjected to both AIF (IL-10) and PIF(IFNγ) signal would robustly exhibit IL-10 driven signal response overriding the IFNγ specific activation of both the STATs. Subsequent experiments quantitatively show that STAT3 dynamics indeed remain strongly IL-10 driven as predicted by the model ( Figure 3B, black filled circles). Notably, the experimentally observed STAT1 amplitude is higher than the predicted amplitude ( Figure S2A) which is quantitively not in line with the initial model prediction. However, we found, despite the differences in amplitude STAT1 dynamics remain strongly comparable between the prediction and validation datasets as a quantiative match betweeen model and data was obtained by multiplying the model trajectory with a scaling factor( Figure 3C). Figure 3D shows the SOCS1 expression dynamics in co-stimulation which is closely comparable to IL-10 only stimulation scenario (compare figure 3D with figure 2D, 3 rd . CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/425868 doi: bioRxiv preprint first posted online Sep. 24, 2018; row, 1 st column). The SOCS1 expression dynamics is also predicted using the 40 independent models with similar goodness of fit obtained using local multistate optimization [41].
Next, we asked how cell-to-cell variability in protein expression [38,39] might influence the robustness of IL-10-STAT3 signal response. It is usually observed that isogenic cells have heterogeneity in their protein expression that subsequently results in heterogeneity in signal-response [42][43][44][45]. Such heterogeneity in signaling dynamics has phenotypic consequences; for instance, individual cells with a certain signalling state can become cancerous and in rare cases can even become drug resistant [46,47]. Hence signal processing at the level of individual cells remain as an important topic of investigation both computationally [48] as well as experimentally [14,44].  Figure 4 shows the normalized STAT3 dynamics from C0, Rc-t1 and Rc-t2 subpopulation (for 1000 cells in each subpopulation) with their respective median and standard deviations.
To mechanistically understand how the reciprocal subpopulation emerges we analysed the protein concentration distribution in the three subpopulations to uncover the sensitivity of single-cell variables (protein concentrations) specific to a given subpopulation [14]. Figure 4B compares the distribution of the protein concentration for C0, Rc-t1 and Rc-t2 cell types. The analysis show Rc-t1 subpopulation has a high level of IL-10R and low level of IL-10Ri compared to their respective values in . CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. Subpopulation responses can be altered by tuning only one/two critical single-cell variable Experimental studies show protein concentration mixing [43] due to processes like cell division and the inherent noise associated to expression/degradation of proteins in individual cells [49][50][51] leads to stochastic changes in the cellular states of individual cells. To understand the significance of single-cell variables with distinct distribution in the three subpopulations, we next performed a set of perturbation studies where we perturbed the signalling states of cells in one subpopulation by systematically replacing the single cell variables from a subpopulation with different signalling state.
As the cells in different subpopulations emerge as a function of cell-to-cell variability, this analysis is designed to systematically understand the effect of expression noise on the most sensitive proteins of the pathway. For instance, we took one C0 cell and replaced a single cell variable (like STAT1) with its counterpart from an Rc-t1 cell( C0 -> Rc-t1). To obtain good statistics a variable in one cell of C0 subpopulation is replaced with its counterpart from 1000 other cells from Rc-t1 and this was repeated for each cell of C0 subpopulation; we show 100 representative cells in the heatmaps ( Figure 5). We studied the effect of such perturbations in all pairwise combinations of the subpopulations (C0 -> Rc-t1 , Rc-t2 ->C0 , C0 -> Rc-t2 , Rc-t2 -> Rc-t1 , Rc-t1 -> Rc-t2). In Figure 5A each row shows the median fraction of Rc-t1 cells transformed to a C0 cell types when a given variable in the Rc-t1 cell is replaced by its counterpart from 1000 C0 cells. The analysis show, emergence of the two reciprocal subpopulations Rc-. CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/425868 doi: bioRxiv preprint first posted online Sep. 24, 2018; t1 an Rc-t2 can be attributed to the differences in the values of a small number of single-cell variables. The protein concentration distributions already show a low basal SOCS1 concentration (SOCS1 B ) in the Rc-t2 subpopulation ( Figure 4A, 2 nd row, 1 st column), hence, in Rc-t2 ->C0 relatively more frequent occurrences of C0 cells are observed when the SOCS1 concentration in Rc-t2 is replaced by its counterpart from C0 cells( Figure 4C, 2 nd row). The analysis suggests that SOCS1 B is a key determinant for the emergence of Rc-t2 subpopulation. The perturbation analysis also show significance of IFNR for the Rc-t2 -> Rc-t1 transformation ( Figure 5A, 2 nd row, 2 nd column), which is not evident only from the protein concentration distributions (Figure4B). Next, as the ratio ILR0R/IL-10Ri is distinct in Rc-t1 we studied the effect of replacing the ratio from the other two subpopulations; figure 5B(1 st row) shows the dramatic increase in the median frequency of Rc-t1 -> C0 transition when both ILR0R and IL-10Ri were simultoneously replaced from C0 to Rc-t1. Similarly, when both IFNR and its inhibitor SOCS1 B is replaced between the pair of reciprocal subpopulation, the frequency of Rc-t2 -> Rc-t1 or Rc-t1 -> Rc-t2 transitions dramatically increased (compare figure 5B, 2 nd row and figure 5A, 2 nd row 2 nd column; figure 5B, 3 rd row and figure 5A, 3 rd row 2 nd column).
Hence our single cell analysis shows, a combinatorial change of key single-cell variables driven by the cell-to-cell variability can result in cell specific high STAT1(low STAT3) or a high STAT3(low STAT1) signalling and the robustness of IL-10-STAT3 signaling axis as observed in the cell population level would be lost, indicating a potential to generate opposing cellular responses within a population of isogenic cells.

Discussion
IL-10-STAT3 and IFNγ-STAT1 are observed to be functionally opposing a spectrum of biological processes; activation of macrophages is enhanced by STAT1 and inhibited by STAT3, cell proliferation is inhibited by STAT1 and promoted by STAT3, and in Th differentiation, STAT1 promotes Th1 responses and STAT3 inhibits Th17 response [52][53][54][55][56][57][58]. Transcriptionally active STAT1 induces death receptor expression to promote apoptosis and it negatively regulates the expression of several oncogenes [59,60]; in contrast, constitutive activity of STAT3 is essential for the survival of many primary tumor cells [59]. Transcriptional targets of STAT3 are also many anti-apoptotic genes that promote tumor cell proliferation [59,61]. Further, dynamic response of the STATs are shown to be critical in determining the cellular responses: IFNγ-STAT1 signaling or IL6-STAT3 signaling is transient and pro-inflammatory but anti-inflammatory IL-10-. CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/425868 doi: bioRxiv preprint first posted online Sep. 24, 2018; STAT3 signaling is sustained [16,23] in nature. However, despite the ability to activate functionally opposing cellular fates, experimental investigations unravel that STAT1 is also activated during IL-10 signaling and STAT3 is also activated during IFNγ signaling [ 18,19,[62][63][64]. In this study, we firstly asked (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/425868 doi: bioRxiv preprint first posted online Sep. 24, 2018; population level led us to investigate if similar robustness in STAT3 responses could also also be observed at the level of the single cells. Recent studies show, the signalresponse at the level of single cells can determine an individual cell's fate in response to the external signal which often results in cell specific gene expression [50,65] or phenotypic outcomes [14].
To computationally investigate the consequences of protein expression heterogeneity in S/1/3 signaling in co-stimulation we next introduced protein expression heterogenity in the population-level model [14,28] (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.  Lysates were centrifuged (10,500 rpm, 10mins) and supernatants were collected.
Protein was quantified by using the Bradford reagent (Pierce) and an equal amount of protein was run on SDS-PAGE. Resolved proteins were blotted to PVDF (Millipore) and then blocked with 5% nonfat dried milk in TBST [25 mM Tris (pH 7.6), 137 mM NaCl, and 0.2% Tween 20]. Membranes were incubated with primary antibody at 4˚C overnight, washed with TBST, and incubated with HRP-conjugated secondary Ab.
. CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

Mathematical model building:
The pathway schemes in Figure 2A and 2C and 3A were be converted to a set of ordinary differential equations (Supporting Material) which captures the dynamics of signalling in both FNγ and IL-10 pathways. We have both the pathways built within one model where we preferentially switch on either the IFNγ or IL-10 pathways for individual pathway stimulation scenarios, or activate both the pathways simultoneously in a costimulation scenario. Below we explain the model rections specific to each of pathway as well as reactions common to both the pathways.

A. IFNγ pathway
The ligand (IFNγ) bindss to the receptor (IFNR) to form an active signaling complex IFN-LR. We lumped several steps of receptor activation [21,22] into a one step receptor activation process [23] assuming reversible kinetics .
. CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

B. IL10 pathway
Similar to the IFNγ pathway the IL10 receptor also binds to ligand and become functionally active. In our simplified model of the pathway the explicit receptor activation deactivation steps are simplified to a one step activation and deactivation process .
The acive signaling complex IL10-LR phosphorylates and activates STAT3 and STAT1(explained in section 'STAT1 and STAT3 signaling', below). A receptor level inhibitor that is observed to act by targeting the IL10R1 [25] is considered in our model as the negative regualtor of IL10-LR. . CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/425868 doi: bioRxiv preprint first posted online Sep. 24, 2018; As observed experimetnally [25], IL10Ri production and degradation is implememted in our model as stimulation independent proesses; reaction (12) and (13) depicts the basal production and degradation of IL10Ri. The differential equations below captures the dynamics of IL10 receptor activation and its inhibition by IL10Ri.

II. STAT1 and STAT3 signaling
In addition to equations [2], [3], [7] and [8] representing the activation of IFNγ and IL-10 specific activation STAT1 and STAT3, in the third condition, co-stimulation, STAT1 and STAT3 are activated by both the stimulus types simultoneously. is implimented in our model. In the same lines, we implimented competition of STAT1 and STAT3 for access to IL10 (see differential equations for x8 and x10). STAT1_act and STAT3_act are the activated/phosphorylate froms of STATs which subsequently undergo dimeriazation and in turn become transcriptionally active [21] STAT1_act + STAT1_act -> STAT1_act_Dm STAT3_act + STAT3_act -> STAT3_act_Dm

III. Transcriptional induction of SOCS1 and SOCS3
Transcriptional induction of SOCS1 and SOCS3 were experimetnally tested upon both IFNγ and IL10 signaling. Both the SOCS are induced relatively strongly upon IL10 signaling compared to IFNγ signaling ( Figure 2B and 2D), especially SOCS1 induced upon 1L10 signaling is ~3 fold higher compared to IFNγ signaling. . CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

Model equations
The differential equations below captures the information propagation in both IFNγ and IL10 pathways. The xs' are model species and the ps are model paramters, names of the species and paramter as well as their bestfict values are detailed in supplementary table TS1. = x5 * p2 -x1 * x2 * p1 * p44 + p46 * p47 * p44.
. CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
. CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The model comprising both the pathways is calibrated to the experimental data ( Figure 2), and further, the calibrated model is used for making predictions that we validated experimentally( Figure 3). Details of model calibration and validation steps can be found in supplementary file S1.

Single cell simulations
To convert the calibrated population average model to single cell level we adopted an ensemble modelling approach in which a population of single cells are generated by sampling the protein concentration assuming protein expression noise at the level of single cells follow a log-normal distribution [14,43]. Several studies show the distribution of protein concentration in individual cells were drawn from lognormal distributions [28,29,43,46], where, for a given protein its population average value from the best-fit model is used as the median of the generated single cell population.  represents an inhibitor of IL-10 signaling which acts by sequestrating and degrading . CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.    (A) Protein concentrations subjected to cell-to-cell variability were systematically replaced between the 6 pairs of subpopulations. In each heat map the x-axis direction shows the 100 representative cells from a subpopulation subjected to such parameter replacement. The notation Rc-t1 → CO, for instance, would mean single-cell variables in Rc-t1 cells are replaced with their counterparts from C0 cells, following which, we calculate the fraction of CO cell emerging from the Rc-t1 subpopulation. Each column in a row, for example in the 1 st row 1 st column of the heat map titled "Rc-t1 → CO", shows the number of C0 cells emerging per Rc-t1 cell when IL-10Ri is replaced from . CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/425868 doi: bioRxiv preprint first posted online Sep. 24, 2018; 500 other C0 cells; the result is normalized to show the fractions where fraction = 1 means replacement of IL-10Ri from 500 C0 cells to Rc-t1 would result in 500 C0 cells.
For all the heat maps the same procedure is followed, the title of the heat map shows the pair of subpopulation studied and the y-direction shows the protein concentrations that were replaced between the pairs. For each protein concentration, 100    Figure 3C. (B) Representative immunoblots for co-stimulation is shown. The experimental procedure is same as described in figure S1A or S1B. (C) Kinetics of SOCS1 and SOCS3 expression on stimulation with M dose of IL-10 or IFNγ.
. CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.