Supplementary Information Supplementary Note 1 Jak-stat Model Supplementary Note 1.1 Steady-state Approximation for Jak-stat Kinetics

Our mechanistic understanding of JAK-STAT inhibition is motivated by the biochemistry of small molecule chemical inhibitors. In the following sections we present the three mechanisms of inhibition that we tested; noncompetitive (Supplementary Figure 2), uncompetitive (Supplementary Figure 4), and competitive inhibition (Supplementary Figure 6) [1, 2]. In the main text, and the proceeding models of pSTAT5 inhibition, we invoke a steady state approximation for our quantitative analysis. Based on experimental considerations delineated here, we chose the 15 minute time point, Supplementary Figure 1, as maximally suited for this approximation for two reasons. First, there is a minimal difference in the average abundance of pSTAT5 between the 15 minute and 30 minute time points for drug doses spanning the concentrations presented in the main text. Secondly, a 15-min timespan is sufficiently short an amount of time that complex gene regulation should not be taken under consideration. There exist additional responses to JAK inhibition (e.g. block of cell proliferation or induction of apoptosis) that may alter the signaling patterns of the cells, but these occur on much longer timescales (e.g. multiple hours) and can safely be excluded from our analysis of signaling responses.


Supplementary Note 1.1 Steady-state approximation for JAK-STAT kinetics
In the main text, and the proceeding models of pSTAT5 inhibition, we invoke a steady state approximation for our quantitative analysis. Based on experimental considerations delineated here, we chose the 15 minute time point, Supplementary Figure 1, as maximally suited for this approximation for two reasons. First, there is a minimal difference in the average abundance of pSTAT5 between the 15 minute and 30 minute time points for drug doses spanning the concentrations presented in the main text. Secondly, a 15-min timespan is sufficiently short an amount of time that complex gene regulation should not be taken under consideration. There exist additional responses to JAK inhibition (e.g. block of cell proliferation or induction of apoptosis) that may alter the signaling patterns of the cells, but these occur on much longer timescales (e.g. multiple hours) and can safely be excluded from our analysis of signaling responses.

Supplementary Note 1.2 Noncompetitive inhibition
Noncompetitive inhibitors bind to their target enzyme irrespective of the presence of the enzymes substrate.
Here, we consider the enzyme substrate to be the measurable abundance of the STAT5 protein. Although, while fluctuations on the ATP levels could potentially affect the enzymatic activities of the kinases under consideraion, we expect ATP levels to be maintained at homeostasis very tightly, such that we could safely assume, without lack of generality, that ATP levels were fixed in our biological system. To develop a quantitative description of noncompetitive inhibition in living cells we start by considering Where we identify the equilibrium reaction parameters as k + for the association of JAK and STAT5, k − represents the dissociation of the JAK-STAT5 complex, and k Ij ± represents the inhibitor forward and reverse rate constants. The enzymatic reactions of the phosphorylation and dephosphorylation are parameterized by rate constants Γ and λ, respectively. Conservation of mass dictates that, In addition we simplify the dynamic equations by using a quasi-steady state approximation, which assumes that the binding / debinding kinetics of the equilibrium reactions are much faster than the catalytic reactions. This separation of time-scales between the equilibrium and the catalytic reactions simplifies the dynamics by setting all the time derivatives of the equilibrium reactions to zero. Using the quasi-steady state approximation in conjunction with our approximated constraints we obtain, Here Incorporating this approximation and noting that Γ sets the time-scale, we change to dimensionless units, dτ = Γdt and Λ = λ Γ , so that In our experimental measurements, we measure pSTAT5 abundance after the system has reached steady state (Supplementary Figure 1). In addition, the measurements are proportional to the abundance of pSTAT5, and not precise concentrations as assumed in our model development. We incorporate these experimental facts by setting the time derivative to zero, introducing a parameter φ to appropriately scale our measurements, and rearranging Eq. 5 we obtain, Here y = [pSTAT5]/k m and represents our pSTAT5 measurement, φx = [STAT5 total ]/k m is the adjusted total STAT5 measurement, and α = [JAK total ] Λkm representing the our pSTAT5 amplitude measurements. To test the mechanism of inhibition we fit our data to the closed form solution Eq. 7 is used to fit our experimental data of the response of pSTAT5 for all levels of total STAT5 and inhibitor dosages. Supplementary Table 1 lists the parameters found to minimize the sum of squared residuals (Supplementary Figure 3). We then used our model and the inferred parameters to derive new coordinates that transform our data to a straight line with a y-intercept of zero and slope of unity, We see excellent agreement between our transformed data and the predicted values (Main text Fig. 2D).  The uncompetitive inhibitor is a small molecule that binds exclusively to its target enzyme bound to its cognate substrate -in this case the JAK-STAT5 chemical complex. Indeed, this mechanism of inhibition manifests in the equilibrium equations of the drug-enzyme interaction and the conservation equation for [JAK total ]. Apart from the inhibitor's dynamic equations our analysis remains the same as the noncompeti-tive scenario, The following mass conservation constraint follows as: Then applying our quasi-steady state approximation as before we obtain an expression for the dynamics of pSTAT5, Where . Solving for the dimensionless variables describe in Equation 6, we obtain This representation provides a transparent understanding of the difference between the uncompetitive and uncompetitive inhibitors. Unlike the noncompetitive inhibitor which reduces the response amplitude exclusively, the uncompetitive inhibitor effectively reduces both the response amplitude and the half effective inhibition concentration. We conclude by presenting the closed form solution used to fit the data (Supplementary Figure 5A) Competitive inhibitors compete with the substrate for binding to the target enzyme. Apart from the difference in the dynamic equations describing the unique mechanism of action, the proceeding analysis utilizes the same approximations as the previous sections. Incorporating the inhibitor mechanism of action into the dynamics 6 As before we approximate the mass conservation equations of constraint as (15b) Applying our quasi-steady state approximation with . As before we proceed by solving for the unitless variables describe in Eq. 6 and arrive at From this equation we can identify the unique functional difference between the noncompetitive and uncompetitive inhibitors with the competitive inhibitor. Unlike the former two equations (Eqs 6,12) representing the JAK inhibitor action, the competitive inhibitor increases the half effective concentration of STAT5.
The closed form solution used to fit our data for JAK inhibitor model selection is Supplementary Figure 7A,B, shows the large systematic error when this model is applied to our experimental data.

Supplementary Note 1.5 Model selection
To find the most likely mechanism of action of AZD1480 in living cells we measured the ability of each model to explain our data (D = {([I jak ], x = STAT5, y = pSTAT5)}). First, we found the optimal parameter set (θ i ) for each model, f i where i ∈ {Noncompetitive, Uncompetitive, Competitive }, by minimizing the total residuals ( ),θ where In this equation N is the number of inhibitor doses and M is the number of CCVA determined STAT5 data points. Once θ i is found for all i models, we selected the optimal model by, In Fig. 2B of the main text we present the total residuals for each model used in Eq. 21. Other metrics (e.g. AIC, BIC) could be used for selecting the optimum model. These techniques are essential for comparing models with different number of parameters. However, no such technique is required for our analysis, because each model has the same number of parameters.

Supplementary Note 2 Coarse grained model of TCR signaling
We developed a coarse-grained model founded upon previous models of T cell receptor (TCR) activation and phosphorylation of ERK [4,5,6]. In our model we simplify the network to five components: 1) peptide bound MHC and TCR complex which we designate "L − R", 2) the SRC family kinase Lck, which is bound the the intracellular domain of CD8, designated "SRC", 3) the active receptor complex as "SRC * ", 4) phosphorylated MEK designated "pMEK", and 5) phosphorylated ERK designated "ppERK". The components interact according the following mass action differential equations: In this model all equilibrium rate parameters are represented using k θ +,− notation, where (+, −) represent the forward and reverse rates, respectively, and θ provides the corresponding reaction identity. Specifically, no superscript represents the L − R and SRC interaction, 'Is' superscript is the SRC and I SRC interaction, and lastly 'Im' superscript is the pMEK and I MEK interaction. The feedbacks are parameterized by the half effective concentration constants for the forward reaction k m+ and the reverse reaction k m− , the hill coefficient of the positive feedback n, and the positive (Γ + ) and negative (Γ − ) feedback strengths. Lastly, all dephosphorylation rates, λ, are identical.
Analysis of the equations are broken into two portions -first the dynamics of [SRC * ], followed by those of [pMEK] and [ppERK]. To analyze the dynamics of [SRC * ] we first define the constraints originating from mass conservation as In our experiments we administer saturating doses of antigen bound to peptide MHC on antigen presenting cells, a reasonable approximation. Furthermore, we apply a quasi-steady state approximation in which all equilibrium reactions are much faster than enzymatic reactions. Together we see that with k src I = k Is − /k Is + . We continue by normalizing our data by k m+ allowing us to rewrite the dynamics of [SRC * ] by the unitless variable x. Specifically we take x = [SRC * ]/k m+ , α = k m− /k m+ , s t = [SRC total ]/k m+ , Γ − /k m+ = Γ − , and I s = I/k src I to obtain Throughout our experiments we do not change the quantity of antigen or antigen presenting cells allowing us to group Γ+[L−R] Γ−k d into a single dimensionless parameter A 1 . Furthermore, we take the feedback amplitudes to be much larger than 1, resulting in 1 Γ+ , 1 Γ − → 0 and a simpler expression for x, We continue by rewriting the expression with respect to dimensionless variables and obtain Where where γ e = Γekm+ k−Γ− and y = [ppERK]/k m+ . By solving the differential equation for m we can simplify our system of three equations to two. Two simple dynamic equations is ideal, because the system properties can be easily interpreted in a single twodimensional phase diagram (Figures 4 and 5 in the main text). The m equation can be described as an integral equation The integral equation representation provides transparency to the systems behavior at the initial conditions, τ = τ o , designated A, and the resulting dynamics, τ > τ o , designated B. This construction is useful to understanding the stream plots in Figure 4 of the main text. In this figure we propose a hypothetical experiment. In this experiment we imagine that it is possible to perfectly manipulate the system, which results in our ability to place it at any coordinate pair (x, m, y). Our question then is, what will be the initial velocity of our system be if we allow it to evolve in accordance to its nature. Mathematically, this hypothetical experiment results in the exclusive evaluation of τ = τ o , which results in our setting B = 0. Substituting for m and γ m γ e /Λ 2 = A 0 , the initial fluxes are described as The Equations 31 represent the instantaneous flux and are used for the stream plots of the [ppERK] versus [SRC * ] (y and x) phase planes in the main text. We proceed by finding the properties of our equations that are experimentally measured -namely the locations and the respective stabilities of the stationary values, 0 = A 1 x n x n + 1 We evaluated the fixed points of x and their corresponding stability numerically in accordance to Eq 32a. The stability of each fixed point was assessed by numerical evaluating the direction of the flux, positive or negative, for small displacements about the respective fixed point. Given our chosen parameter set (Supplementary Table 2), we found our system to elicit bistable behavior [7].      Trametinib_Tube_002.fcs 57.4

Supplementary Note 3.3 pMEK (S221) response to MEKi
To ensure the model is correct, and that MAPK is truly a separate and unidirectional sub-network, we measured the response of phosphorylated MEK (pS221; monoclonal antibody (166F8) from Cell Signaling Technologies) to MEK inhibition by PD325901 (Supplementary Figure 12). We quantified the data in accordance with the analysis of ppERK in the main text. Specifically we applied a two component Gaussian mixture model to quantify the means of the pMEK high and pMEK low components and the fraction of cells occupying the pMEK high state. In Supplementary Figures 15 and 16, we show that pMEK abundance and the fraction of cells occupying the high, and consequently the low, states are invariant to MEK inhibition by PD325901.

Supplementary Note 4 Proliferation
We extend our signaling findings to the proliferative response of cells to MEK and SRC inhibition by analyzing the dilution of amine reactive dyes (CFSE or CTV), cell size, and the abundance of CD8 per cell. Analogous to the binary activation of ppERK signaling, we classified individual cells as either inactive or active by using a binary gate that separated small and CD8 low expressing cells from large and CD8 high expressing cells, respectively. To correlate a reduction in the average abundance of ppERK to a functional decrease in proliferation, we computed a normalized average number of divisions among active cells: every time a cell divides from generation n to n + 1, it contributes an additional sister to the count of n + 1, and this doubling in cell number must be taken into account when computing the average number of divisions a cell has undergone. Hence, we introduced two observables, the fraction of activated cells and the average number of divisions Divisions .
We computed the normalized average number of divisions as: with N total = nmax n=0 N n /2 n . In this expression the 2 n accounts for the contribution of sister cells for the n th division. We then computed the fraction of active cells and noticed that the effectiveness of our binary gate (CD8 low vs CD8 high ) did not report accurately the dynamics of cell proliferation and drug inhibition. Specifically some of the cells gated as inactive (i.e. CD8 low cells are capable of diluting the CTV dye and as a result should have been counted as active (cells cannot divide unless they have been activated). To accurately count these dividing cells as active, we devised a strategy for the robust estimation we derived a correction for the gating strategy. Ultimately, we reported Divisions active and the Robust Fraction Activated in Figure 6C of the main text: these observables best encapsulated the dynamics of lymphocyte response in our system.