Single cell suspensions of mouse splenocytes were prepared as follows. Spleens were harvested from 4-month-old female C57BL/6J mice (Strain 00064 from Jackson Laboratory, Bangor ME) and mechanically dissociated with a sterile syringe plunger. We performed a 1 min ACK lysis to remove any red blood cells. Splenocytes were then washed in PBS and exposed to a solution of 0.1M glycine in PBS (with adjusted pH at 4.0) for 1min on ice, followed by 2 washes with complete medium. Cells were then rested in complete medium for 2 hours at 37°C.

Cells were then exposed to a serial dilution of IFN-α2 or IFN-β (R&D systems, Minneapolis MN) in complete medium and incubated for varied amounts of times. At each time point, ice-cold 4% paraformaldehyde solution was added to the reaction culture (for a final concentration of 1.6% paraformaldehyde) and cells were left in fixative for 15min on ice. Cells were then spun, their supernatants were discarded, and cells were then permeabilized with a 15-min incubation in 90% methanol on ice. Cells were then washed twice with FACS buffer (4% FBS in PBS + 0.1% sodium azide) and stained with primary anti-phospho-STAT1 (clone 58D6 from Cell Signaling Technologies, Danveres MA) for 30min at room temperature. Cells were then washed once with FACS and stained with a cocktail of antibodies against surface markers and Fab anti-Rabbit secondary antibody (see Table S1 ) for 30min at room temperature. Cells were then washed once with FACS buffer and resuspended in FACS buffer with DAPI (1µMol) and acquired on a 20-channel Fortessa cytometer (BD Biosciences). Single cell gating for live CD19+ B cells and levels of phosphorylated STAT1 were exported as geometric mean for further analysis ( Figures S1, S2 ).

Figure 1A illustrates a streamlined model of IFN which includes receptor assembly, STAT phosphorylation, receptor internalization, and SOCS1 inhibition. The model aggregates various biochemical sub-steps, such as STAT association and dissociation with active receptor complexes, into single-step reactions. Since only pSTAT1 was measured experimentally, our model does not describe STAT1-STAT2 dimers explicitly but rather uses a single STAT variable which can be phosphorylated to an active pSTAT form. Since the JAK kinases are constitutively bound to their corresponding receptor subunits, we do not model them explicitly and instead consider the ternary complex to phosphorylate STAT to pSTAT (28). A more detailed model which explicitly models additional biochemical steps is presented in the SI and yields similar predictions for the pSTAT response ( Figure S3 ), indicating that the minimal model of Figure 1 encapsulates the salient features of the signaling network. We include basal and IFN-induced receptor internalization in our model as well as recycling of the receptor subunits back to the cell surface, assuming that the internalized ligand is degraded (36). The model includes negative feedback by SOCS1, which is produced downstream of pSTAT and binds IFNAR1 to block STAT phosphorylation. It is known that IFNAR1 and IFNAR2 surface expression levels vary between cells, and in practice we found it necessary to model this variation using a log-normal distribution of expression levels for IFNAR1 and IFNAR2 (62). All model predictions are made by taking 30 samples from these receptor distributions and averaging the predicted pSTAT responses. This is a simple way to model heterogeneity of receptor expression, a factor we explore more in depth below.

The complete mathematical description of our model using a system of ordinary differential equations (ODEs) is provided in Tables S2, S3 and is solved numerically using the PySB python package (63). We used ligand-receptor binding and unbinding rates measured for human IFNs and applied detailed balance to obtain in-membrane rates (see Table S2 ) (10, 38). A global scale factor of 1.5 was applied to convert the number of pSTAT molecules predicted by the model into units of MFI measured by experiment. Estimates of the values of model parameters which have not been directly measured were obtained using the systems biology package PyDREAM to perform Markov Chain Monte Carlo (MCMC) simulations ( Figure S4 ) (64). The mean absolute error (MAE) between the data and the average prediction of our best fit model is 40%. A bootstrap analysis of the fitting procedure using an 80-20 train-test split of our experimental data yielded a similar MAE across bootstrap batches, further validating our parameter fit and demonstrating the robustness of the fitting procedure.

To fit the phenomenological models of anti-viral and anti-proliferative activity to the data in Section Proximal Receptor Signaling Maps to Biological Activity, we used standard nonlinear function fitting procedures provided by the scipy Python package which are based on the Trust Region Reflective algorithm. The fitting procedure for Eq. 5-6 yielded a root mean square error of 11%.

Many aspects of IFN signaling are captured by the relatively simple equilibrium description of receptor complex formation obtained from the steady state of our ODE model, providing qualitative insight into how various components of this system affect signaling dynamics. The equilibrium solution for the number of active ternary complexes formed in response to a ligand can be found by solving the chemical kinetic equations along with the detailed balance condition K ₁ K ₃ = K ₂ K ₄ (see below), where K ₁ and K ₂ are the equilibrium dissociation constants for IFN binding to each of IFNAR1 and IFNAR2, and K ₄ and K ₃ are the in-membrane dissociation constants for recruitment of IFNAR1 and IFNAR2 respectively to the ternary complex (see Figure 1A ) (53). The result is an expression for the equilibrium number of ternary complexes (53, 65):

where Ξ ( I ) = 1 + K 4 R T ( I + K 2 ) ( I + K 1 ) I   K 1 , I is the extracellular IFN concentration, Δ = R 1 T − R 2 T is the difference between the (two-dimensional) cell surface densities of R ₁ (IFNAR1) and R ₂ (IFNAR2), R T = R 1 T + R 2 T is the total surface density of receptors of both types. The detailed balance condition arises from the fundamental physics relating the dissociation constants to the binding energies of the ligand with the receptor subunits, E ₁ and E ₂. In the simplest case, K 1 = 1 v e − E 1 k T and K 2 = 1 v e − E 2 k T , where v is a small volume roughly proportional to the average volume of the reacting molecules; for three-dimensional binding affinities, it is customary to set v = 1M ^–1 ≈ 1.75nm ³ (1, 66). The in-membrane affinities K ₃ and K ₄ obey the same physics; in the ideal case where the ligand binding to one of the subunits does not affect its binding to the other one, K 4 = 1 a e − E 1 k T and K 3 = 1 a e − E 2 k T where a is the average area of the reacting molecules (67). The detailed balance condition also holds for more complicated dependencies of the dissociation constants on the molecular interaction energies as long as the binding-unbinding reactions involved in the receptor dimerization do not require direct input of metabolic energy (68).

The ternary complex dose-response curve is a non-monotonic function of the ligand concentration. The maximal response T_max is attained at I m = K 1 K 2 . Crucially, weak affinity ligands saturate to a lower maximum receptor activity than high affinity ligands. At concentrations below the maximal response, the dose-response curve can be approximated with a Michaelis-Menten type function:

where K_eff is the effective equilibrium dissociation constant of ternary complex formation, with K_eff given by:

where D = R T K 4 + 1 + 10 X 1 2 + X + 3 R T K 4 B 2 − 1 + ( Δ / R T ) 2 , B = 1 + ( K 4 / R T ) ( 1 + X ) 2 and X = K ₂/K ₁. It is important to note that Equations (1-3) predict contours of constant K_eff and T_max as a function of receptor binding strengths, as shown in Figure 2A .

The product of the dissociation constants K ₁ × K ₂ has been empirically observed to correlate well with K_eff for a wide range of IFNs (39). This product is therefore expected to be approximately proportional to related biological IC ₅₀ values and other quantities related to K_eff , although there is a systematic bias in using K ₁ × K ₂ in place of K_eff (see also Figure S5 ).

Ternary complex formation leads to the production of phosphorylated STAT1-STAT2 dimers at rate k p + , which are constitutively de-phosphorylated at rate k p − . In the limit that intracellular STAT equilibrates with the active receptor complexes faster than STAT is phosphorylated, the equilibrium copy number of pSTAT1-pSTAT2 dimers is well approximated by a Michaelis-Menten equation:

where S_T is the total number of STAT molecules in the cell, K p = k p − / k p + , and A is the cell surface area. Thus, the maximal magnitude of the pSTAT response is p S T A T max = S T T max T max + K p / A and the half-maximal response of pSTAT is pSTAT EC ₅₀ = K_eff /(1 + T_maxA/K_p ). Unlike the effective affinity of the receptor binding K_eff , the EC ₅₀ of the pSTAT response explicitly depends on the cell area and not only on the surface densities of the receptor sub-units.

To validate our model at the level of receptor assembly (Equation 1), we compared the predicted number of ternary complexes for a variety of IFNα2 mutants to previously reported measurements that used single-molecule tracking of receptor subunits, as reported in (10). IFNAR1 affinity has been shown to be the main differentiating factor between these IFNs (39, 43, 69). We demonstrate in Figure 2B that differences in the level of maximum ternary complex formation between IFNs are well explained in our model of receptor assembly using only the differences in IFNAR1 binding strength between these IFNs [i.e., holding IFNAR2 binding strength constant at the measured IFNα2 wild type level (38)]. This agreement between the model and data for a wide variety of affinity-altering mutations of the IFN molecule demonstrates that our model captures the essential features of receptor assembly kinetics in Type I IFN signaling.

Further model validation is provided by looking at the effect of two mutations to IFNα2 on the pSTAT dose-response curve. The mutants IFNα2-R120A and IFNα2-M148A have reduced affinities for IFNAR1 and IFNAR2 respectively and have been used to explore the dependence of the pSTAT response on IFN-receptor binding strength (10). Mutations of the IFN ligand alter the binding energy and therefore the dissociation rates of all reactions involving IFN-receptor interaction due to the detailed balance condition described above. In particular, reduction in IFNAR1 affinity alters both K ₁ and K ₄, and a reduction in IFNAR2 affinity alters K ₂ and K ₃ (see Mathematical Methods). The IFNα2-R120A ligand, with a 60-fold reduced affinity for IFNAR1, was shown to alter both pSTAT_max and the pSTAT EC ₅₀ while the IFNα2-M148A ligand, with a 50-fold reduced affinity for IFNAR2, only affected the pSTAT EC ₅₀. In Figure 2C we show that our equilibrium model recapitulates the effect of these changes in binding energy on the pSTAT dose-response curves. The difference in the effects of these two mutations is not surprising considering we already showed in Figure 2A that, in the range of K ₁ and K ₂ affinities observed for Type I IFNs, T_max – and therefore pSTAT_max – is largely insensitive to changes in IFNAR2 affinity while K_eff – and therefore pSTAT EC ₅₀ – is sensitive to both IFNAR1 and IFNAR2 affinity. The difference in pSTAT_max between IFNα2 and IFNα2-R120A, shown in Figure 2C , demonstrates that ligand affinity is sufficient to enable amplitude-based absolute discrimination. This effect has also been observed in other signaling pathways (55, 70, 71). Additionally, the results in Figures 2A–C suggest that pSTAT_max and the pSTAT EC ₅₀ can be tuned independently by ligand binding strengths. We will return to this point in the Discussion.

We next focused on the pSTAT response to IFNα2 and IFNβ in particular, explicitly including the effects of short-term and long-term negative feedbacks in order to capture changes in the response over time. To this end, we fit the free parameters of our ODE model (the rates of STAT phosphorylation/dephosphorylation, of SOCS1 production, of receptor internalization, and receptor subunit expression levels, summarized in Table S3 ) to our measurements of primary mouse B cell response to IFN stimulation (see Methods for details). Figure 2D shows that our model can be fit closely to our experimental data. The MCMC fitting did not tightly constrain the free model parameters, which is to be expected due to “sloppiness” observed in most systems biology models as a result of multiple compensatory feedbacks in these systems (72, 73).

Both the EC ₅₀ and pSTAT_max of the pSTAT dose-response curve depend on the stimulation time but converge at long times to steady state values for each IFN, and the values predicted by the ODE model are in agreement with the experimentally determined ones ( Figure 2E ). Interestingly, we do not observe a significant difference in the pSTAT_max between IFNα2 and IFNβ. This is somewhat surprising because there is a noticeable difference in the expected T_max between these IFNs (as shown in Figures 2A, B ), and furthermore our equilibrium model naïvely predicts an approximately 15% lower pSTAT_max according to Eq. 4, for reasonable choices of parameters (R_T = 2000 subunits, A = 450 μm ², affinities and K_p from Table S3 ). The fold-change in pSTAT_max between IFNs can be estimated from Equation 4 as

For large K_p /A, the difference in response is expected to be the same as the difference in T_max . In contrast, for small K_p /A no difference in pSTAT_max is expected despite a difference in T_max . As we will show below, additional factors besides the cell surface area and phosphorylation constant (A and K_p ) also affect the difference in pSTAT_max . We note here that the lack of a region of absolute discrimination between IFNα2-induced and IFNβ-induced pSTAT response in our experimental results in primary mouse B cells is consistent with some of the measurements performed in a variety of human cell types exposed to IFNα2 and IFNα2-YNS, an IFNβ mimic (39). On the other hand, there are also reports of differences in maximal pSTAT response between IFN subtypes in other settings (35, 39, 42). This is reasonable since different cell types have different values for K_p , A, and other feedback rates which may tune the difference in pSTAT_max . Thus, in the sections below we explore additional factors besides ligand affinity which can regulate the pSTAT response and allow for functional plasticity.

Receptor density on the cell surface has been suggested as an important factor affecting the specificity of response (35, 74, 75). As shown in Figure 3A , the theoretically predicted separation between pSTAT_max for ligands of different affinities (i.e., the region of absolute discrimination) remains nearly constant except at very low receptor expression levels – because the dominant effect of a change in the number of the surface receptors is the overall increase in the number of ternary complexes. For this reason, total receptor expression is not expected to affect amplitude-based absolute discrimination.

By contrast, as shown in Figure 3B , the theoretically predicted pSTAT EC ₅₀ decreases monotonically with increasing total number of receptor subunits and can shift by several orders of magnitude over a biologically realistic range of receptor expression levels [e.g. 10² versus 10³ copies, (18)]. This suggests that modulation of total IFNAR1 and IFNAR2 expression levels could make a cell preclude weaker binding IFNs from activating downstream responses that require high intracellular pSTAT levels. This is conceptually similar to absolute discrimination but without complete decoupling of ligand concentration from the pSTAT response. For physiological numbers of receptor subunits (e.g., ~10³ total receptor subunits), preventing a high pSTAT-dependent response to IFNα2 would require IFN concentrations to remain well below 1 nM at all times. Note, however, that this does not provide discrimination between two distinct ligands requiring different cellular responses because the weaker binding ligand is effectively never detected. Rather, this effect could provide a way of selectively responding only to high affinity ligands as a function of receptor expression levels. This could be used to selectively target specific cell types, a point we return to in the Discussion.

Figure 3C shows that asymmetric expression of IFNAR1 and IFNAR2 is predicted to be an important factor controlling the difference in response for different ligands. The greatest degree of separation in pSTAT_max is achieved for equal levels of IFNAR1 and IFNAR2, and this separation decreases monotonically and symmetrically for a growing imbalance in expression levels. The symmetry in the effect of expression imbalance may seem counter-intuitive given that IFNAR1 affinity is the major discriminant between IFN subtypes. However, Eq. 1 shows that the response only depends on IFNAR1 abundance relative to IFNAR2 abundance (i.e., through the terms R_T and Δ), so the separation in pSTAT_max between different IFNs remains nearly constant unless an imbalance between IFNAR1 and IFNAR2 is introduced.

Although we could not directly probe the effect of tuning receptor subunit expression on the pSTAT response in primary mouse B cells, this effect has been measured in WISH cells by siRNA knockdown experiments (45). In Figures 3D, E we compare the measured and simulated pSTAT1 response after siRNA transfection and subsequent treatment with either IFNβ analog IFNα2-YNS or IFNα2. Both model and data exhibit a reduction in pSTAT level that varies approximately linearly with the reduction in IFNAR1 levels for both IFNs. In comparison, the reduction in pSTAT is lower at equivalent levels of IFNAR2 knockdown in both model and experiment. The good agreement between our model and data, especially considering the difference in cell type between the siRNA experiments and those used to parameterize our model, demonstrates that our computational model captures essential factors regulating IFN signaling.

To further investigate the dependence of pSTAT_max and pSTAT EC ₅₀ on the receptor subunit numbers, we used the cell size as a proxy for the receptor number. We grouped B cells from our flow cytometry experiments by their forward scatter size as a proxy for cell size. Cells in the top 20^th percentile in forward scatter (i.e., large cells) were experimentally observed to have a greater pSTAT MFI in response to both IFNs compared with those in the bottom 80^th percentile ( Figures 3F, G ). To compare experimental data with our model quantitatively, we varied the cell size in our model while keeping the surface concentration of receptor subunits and the cytoplasmic concentration of STAT molecules constant (76). The experimentally observed separation in the pSTAT_max in small and large cells was explained in our model by assuming fixed effective cell radii of 6.5 μm and 8 μm (corresponding to roughly two-fold difference in the cell volume) for the small and large cells respectively. These cell groups may correspond to cells in different phases of the cell cycle or other factors which affect cell size (77). Importantly, none of the rate constants for the biochemical reactions in our model needed to be adjusted from our original fitting procedure (in Figure 2D ).

Both receptor internalization and feedback by SOCS1 have been discussed in the literature as factors relevant for functional plasticity (36, 38, 39, 46–49). Our computational model allows us to compare the role of these two short-term feedback mechanisms in Type I IFN signaling. Since both receptor internalization and SOCS1 upregulation happen rapidly in response to IFN stimulation, there may be partial redundancy in their roles as negative feedback mechanisms (36, 78).

It has been shown that the rate of receptor internalization may be significantly increased in response to IFN stimulation and that this effect is much stronger in response to IFNβ than to IFNα2 (36, 39). These works additionally showed that the rates at which IFNAR1 and IFNAR2 are recycled back to the cell surface may also differ. As a result, each IFN induces a ligand-specific change in the levels of IFNAR1 and IFNAR2 present on the cell surface post-stimulation. This suggests that receptor internalization could be a mechanism for regulating functional plasticity, since receptor expression asymmetry is expected to affect signal separation by amplitude-based absolute discrimination (Section Receptor Expression Levels Can Regulate Signaling Specificity and Figure 3C ). In contrast to receptor internalization, SOCS1 inhibits Tyk2 phosphorylation of STAT1, reducing the number of ternary complexes which are actively phosphorylating STAT. We expect this to have a similar effect to reducing total receptor expression, which was predicted not to regulate signal specificity (Section Receptor Expression Levels Can Regulate Signaling Specificity, Figures 3A, B ).

To test if receptor internalization and SOCS1 play non-redundant roles in regulating response specificity, we used our computational model to query how these feedbacks each regulate the course of the pSTAT response. In Figure 3H we plot the predicted dose response curve for IFN signaling without SOCS1 or receptor internalization. By then looking at the isolated effect of SOCS1, we find that both the IFNα2 and IFNβ dose-response curves are reduced by similar amounts by the action of SOCS1. The difference in pSTAT_max between IFNs remains proportionally the same because SOCS1 inhibits STAT phosphorylation rather than ternary complex formation, so that feedback on the system is not dependent on ligand identity. This suggests that SOCS1 does not play a role in the emergence of functional plasticity via amplitude-based signal discrimination.

In contrast to the effect of SOCS1, the effect of receptor internalization on the difference in pSTAT_max between IFNs could be significant for amplitude-based signal discrimination. This can be seen in Figure 3H by comparing the response predicted by our computational model without any feedback to that with internalization. Receptor internalization rates can play a role in the emergence of functional plasticity because they are specific to the type of bound IFN. It has been reported that receptor internalization is greater for IFNβ than IFNα2, and the parameter fitting procedure for our ODE model also yielded greater receptor downregulation in response to IFNβ than IFNα2 (36). Since in the absence of feedback the pSTAT_max of the dose-response curve at 60 minutes is predicted by our ODE model to be greater for IFNβ than IFNα2 (see Figure 3H ), and since receptor internalization decreases each pSTAT dose-response curve, greater receptor internalization following IFNβ stimulation effectively acts to reduce the difference in pSTAT_max between these IFNs. In particular, faster IFNAR2 recycling in response to IFNα2 relative to IFNβ [as has been reported in (36)] diminishes the difference in pSTAT_max which might occur due to differences in receptor affinity. This means that receptor internalization may decrease the specificity of the pSTAT response in Type I IFN signaling. This effect is not an inherent property of receptor internalization, but rather is the result of the particular internalization rates observed in the Type I IFN system. These results also demonstrate that SOCS1 and receptor internalization have non-redundant roles as feedback mechanisms in Type I IFN signaling.

Although we have not observed evidence of amplitude-based absolute discrimination in the pSTAT dose-response curves, absolute discrimination based on the time course of Type I IFN signaling must also be considered. The most common type of time course signal discrimination is based on distinguishing transient from sustained signaling dynamics (57–59), although more complicated inference is possible (79, 80). We investigate time course discrimination as a possible mechanism for distinguishing IFNα2 from IFNβ by comparing early and late time levels of pSTAT. A sustained response exhibits similar pSTAT levels at both early and late times, while a transient response exhibits a different (usually diminished) pSTAT level at late times relative to early times. In Figure 4A we show that the character (transient or sustained) of the pSTAT response in our experiments appears to depend more on the concentration of IFN rather than the IFN identity, with a sustained response for low concentrations and a transient response for high concentrations. However, in either concentration regime, both IFNs seem to follow the same pattern, which suggests that the ligand discrimination does not occur via time course discrimination.

To investigate further, in Figure 4B we re-plot our pSTAT measurements in mouse B cells with the response at early times on the horizontal axis and the response at late times on the vertical axis (at different IFN concentrations). A sustained response would produce points on the diagonal (i.e., the line pSTAT_early = pSTAT_late ) in such a plot, while a transient response would produce points away from the diagonal. Ligand discrimination would present itself in such a plot as a separability of points based on ligand identity. However, as Figure 4B shows, the data points corresponding to two different IFNs are not distinguishable (within the error bars) except perhaps for a narrow window of intermediate concentrations. Thus, we conclude that time course discrimination is unlikely to be important for response specificity in Type I IFN signaling although further studies might be required.

Our results thus far suggest that functional plasticity in Type I IFN signaling does not arise from absolute discrimination of the initial pSTAT response within the first 60 minutes. However, differential refractoriness in the pSTAT response mediated by the protein USP18 offers another potential mechanism to achieve absolute discrimination (52, 81). It has been established that the upregulation of USP18 in response to Type I IFN stimulation over the course of 24-48 hours leads to somewhere between a 15- and 60-fold increase in the in-membrane dissociation rate of IFNAR1 from the ternary complex, with the magnitude of this effect varying between cell types (10, 81). Inhibition by USP18 leads to significantly reduced response to further IFN stimulation, which is relevant for clinical uses of Type I IFN (10, 42, 81). Importantly, IFNα2 typically exhibits much greater refractoriness than IFNβ at equivalent doses, although the degree of deactivation can be partially compensated by increasing the stimulation dose (42, 46). It is therefore possible that differential signaling mediated by USP18 over long timescales could lead to ligand specific differences in the pSTAT_max , enabling absolute discrimination between the two ligands.

Since the action of USP18 occurs on the intracellular domain of the receptor, its effects are represented in the model via a change in K ₄, the in-membrane dissociation constant of IFNAR1 binding to the binary complex IFN : IFNAR2 (see Mathematical Methods). Figure 5A shows that the experimentally observed decrease in the numbers of ternary complexes in cells primed with IFN compared with unprimed cells can be accounted for by a 15-fold increase in K ₄ (10). Crucially, the USP18-induced change in the in-membrane dissociation constant K ₄ not only shifts the dose response curve towards higher concentrations but also decreases the plateau saturation level pSTAT_max , raising the possibility of amplitude-based absolute discrimination. The model also reproduces the effect of USP18 refractoriness on the mutant IFNα2-M148A, as shown in Figure 5B . The agreement between the model and the data validates the simple approach of modeling the presence of USP18 indirectly as a shift in K ₄, at least for high expression levels of USP18.

Interestingly, comparison of Figure 2C with Figure 5B shows that the effect of USP18 (which only alters K ₄) on the pSTAT dose-response curve is recapitulated by the IFNα2-R120A mutant that has lower binding strength to IFNAR1 (which alters both K ₁ and K ₄), as was first described in (10). The similarity of these two perturbations on the dose-response curve indicates that ternary complex formation is highly sensitive to the rate of in-membrane receptor subunit association.

We further examined how USP18 may affect the ligand discrimination specifically between IFNα2 and IFNβ, as shown in Figure 5C , using our computational model. It demonstrates that the action of USP18 is predicted to reduce the amplitude of the IFNα2 dose-response curve much more than the corresponding amplitude of the IFNβ dose-response, in agreement with previously reported experimental results (42) and our own experimental analysis (not shown). The reason for this differential refractoriness is that the high affinity of IFNβ for IFNAR1 still allows formation of ternary complexes even in the presence of USP18, which is not the case for IFNα2 due to its lower affinity to IFNAR1. Most importantly, this differential effect of USP18 increases the separation in pSTAT_max values between the two interferons and leads to the emergence of a region of absolute discrimination, which is expected to be apparent at long times (10-48 hours) when sufficient USP18 has accumulated in the cell. This fact will be important for the potential explanation of the differences in the anti-viral and anti-proliferative activities in Section Proximal Receptor Signaling Maps to Biological Activity.

Interestingly, the time scale for USP18 upregulation is similar to that for the emergence of anti-proliferative activity (~24 hours) (11, 48, 82). By contrast, anti-viral genes are typically upregulated very quickly (less than 8 hours) in response to IFN stimulation, regardless of IFN subtype (11, 83). Motivated by this observation, we now discuss how differential signaling at the level of pSTAT may be translated into differences in physiological activity, focusing on the hypothesis that the differential activity may be linked to the regulation of differential pSTAT response by USP18 mentioned in the last section and Figure 5 .

A multitude of signaling reactions and gene expression downstream of receptor signaling are involved in anti-viral and anti-proliferative responses. Not all of these reactions are known, making detailed modeling of these processes unfeasible. Rather, in this paper we use phenomenological mappings that translate the pSTAT response into anti-viral and anti-proliferative activity, as measured by common experimental functional assays, to gain insights into the emergence of functional plasticity based on the specificity of the receptor and STAT phosphorylation levels (38–40). Despite substantial variability in the results of anti-viral and anti-proliferative assays between different studies, our conclusions remain qualitatively robust across cell types and experimental variation.

Our aim is to capture the following notable features distinguishing anti-viral and anti-proliferative activity. First, the anti-viral IC ₅₀ (the IFN concentration required for half-maximal anti-viral activity) is substantially lower than the pSTAT EC ₅₀ [ Figure 6A (left)]; this is not the case for anti-proliferative activity [ Figure 6A (right)] (38, 39). Second, complete inhibition of viral replication is achievable by all IFN subtypes that have been investigated while the maximum inhibition of cell proliferation differs between IFN subtypes by as much as ~20% (38–40).

The observed discrepancy between the IC ₅₀ of the anti-viral response and the pSTAT EC ₅₀ indicates that the components of the JAK-STAT pathway downstream of pSTAT that mediate anti-viral gene transcription require very little pSTAT to become activated. We found that despite its unusual properties, the dependence of the anti-viral response on the pSTAT levels can be reproduced by a simple Michaelis-Menten type dependence:

where pSTAT(IFN) is the intracellular concentration of pSTAT as a function of the extracellular IFN concentration, and K_M is the effective Michaelis-Menten constant for the intra-cellular processes leading to anti-viral activity. As shown in Figure 6A , anti-viral activity described by Eq. 5 saturates at much lower IFN concentrations than the corresponding pSTAT dose-response curve (c.f. Figure 2D ) for sufficiently small K_M (i.e., K_M ~pM); we used K_M = 0.8pM in Figure 6A .

The anti-viral IC ₅₀ can be found by substituting Eq. 4 for the pSTAT response into Eq. 5 and solving for the IFN concentration which gives half-maximal activity. In the limit of K_M /[S_T ] ≪1, the IC ₅₀ is proportional to (pSTAT EC ₅₀) × K_M /[S_T ], recalling that [S_T ] is the concentration of intracellular STAT (our model used [S_T ] = 0.7 nM). It follows from this proportionality that a small K_M guarantees that the anti-viral IC ₅₀ is lower than the pSTAT EC ₅₀. This also explains why the anti-viral IC ₅₀ in this regime is essentially the same for all IFNs because differences in pSTAT EC ₅₀ between IFNs are also shrunk by the factor K_M /[S_T ].

We now turn to anti-proliferative activity. A phenomenological mapping for anti-proliferative activity must preserve differences in maximal anti-proliferative response between IFN subtypes, as apparent in the high concentration regime of Figure 6A (right), as well as the difference in the anti-proliferative IC ₅₀ between IFNβ and the other IFN subtypes [also Figure 6A (right)]. These differences between IFNs are only observed in anti-proliferative assays, not anti-viral assays. Since the anti-proliferative activity emerges on longer time scales we include the effect of USP18 on the pSTAT response, denoted as pSTAT_primed (IFN) (29, 31, 52). Unlike the anti-viral response that could be modeled by a single Michaelis-Menten type function, the anti-proliferative response required a more complex function of the following form to explain the data in Figure 6B :

where γ ₁ and γ ₂ are phenomenological Hill coefficients and K_{M 1}, K_{M 2} are the effective anti-viral EC ₅₀’s for each contributing term. Fitting Eq. 6 to the data, using our model of pSTAT induction presented in Section Validation of the Model and the Effects of Ligand Affinity on Response Specificity and measured IFN binding affinities (39), results in the best fit parameter values K_{M 1} = 72 pM, K_{M 2} = 158 pM, γ ₁ = 0.9, γ ₂ = 3.7. For sufficiently large values of K_{M 1} and K_{M 2} – on the order of 100 pM – differential refractoriness expressed in the different saturation levels between the strong binding IFNα2-YNS and the weaker binding variants results in the differential saturation of anti-proliferative activity at high concentrations in Figure 6A (right). Thus, absolute discrimination on the pSTAT level shown in Figure 5C translates into functional plasticity on the level of anti-viral and anti-proliferative responses. This result is robust to variation in the parameters of Eq. 5-6, not simply the result of a choosing specific parameter values.

To demonstrate the robustness of our phenomenological mappings, in Figure 6B we compare the predicted IC₅₀ values of both the anti-viral and anti-proliferative responses to previously reported measurements for several IFNα2 mutants spanning four orders of magnitude in K_eff (39). We find that our phenomenological models semi-quantitatively recapitulate both anti-viral and anti-proliferative IC₅₀’s for most IFNs, although agreement was somewhat worse for the quadruple mutant IFNα2-YNS-L153A as well as for the weakest binding IFNs IFNα2-L30A and IFNα2-R149A. The product K ₁ × K ₂ using affinities measured for each IFN is provided in Figure 6B for context, and provides a good estimate of each IC ₅₀ other than the high affinity mutant IFNα2-YNS (39). This quantity has been used as an approximation for K_eff in previous work, and we show in Figure S5 that this approximation is reasonable for a realistic range of K ₁ and K ₂. The advantage of our phenomenological mappings is that they can predict the biological activity at any concentration, and the maximal biological activities, for each IFN subtype; the quantity K ₁ × K ₂ alone cannot provide such estimates.

The order of magnitude difference between K_Mi for anti-proliferative activity (~100 pM as obtained from the best fit values for Eq. 6) and K_M for anti-viral activity (less than 1 pM), as well as the Hill coefficients greater than 1 in Eq. 6, suggest that different molecular processes are responsible for the activation of these two types of cellular response. We return to this point in the Discussion.