In Silico Modeling of Itk Activation Kinetics in Thymocytes Suggests Competing Positive and Negative IP4 Mediated Feedbacks Increase Robustness

The inositol-phosphate messenger inositol(1,3,4,5)tetrakisphosphate (IP4) is essential for thymocyte positive selection by regulating plasma-membrane association of the protein tyrosine kinase Itk downstream of the T cell receptor (TCR). IP4 can act as a soluble analog of the phosphoinositide 3-kinase (PI3K) membrane lipid product phosphatidylinositol(3,4,5)trisphosphate (PIP3). PIP3 recruits signaling proteins such as Itk to cellular membranes by binding to PH and other domains. In thymocytes, low-dose IP4 binding to the Itk PH domain surprisingly promoted and high-dose IP4 inhibited PIP3 binding of Itk PH domains. However, the mechanisms that underlie the regulation of membrane recruitment of Itk by IP4 and PIP3 remain unclear. The distinct Itk PH domain ability to oligomerize is consistent with a cooperative-allosteric mode of IP4 action. However, other possibilities cannot be ruled out due to difficulties in quantitatively measuring the interactions between Itk, IP4 and PIP3, and in generating non-oligomerizing Itk PH domain mutants. This has hindered a full mechanistic understanding of how IP4 controls Itk function. By combining experimentally measured kinetics of PLCγ1 phosphorylation by Itk with in silico modeling of multiple Itk signaling circuits and a maximum entropy (MaxEnt) based computational approach, we show that those in silico models which are most robust against variations of protein and lipid expression levels and kinetic rates at the single cell level share a cooperative-allosteric mode of Itk regulation by IP4 involving oligomeric Itk PH domains at the plasma membrane. This identifies MaxEnt as an excellent tool for quantifying robustness for complex TCR signaling circuits and provides testable predictions to further elucidate a controversial mechanism of PIP3 signaling.

T cells are key mediators of adaptive immune responses. Through a plasma-membrane anchored TCR, they recognize pathogen-derived peptides bound to Major Histocompatibility Complex proteins (pMHC) on the surface of antigen-presenting cells. TCR engagement triggers activation, proliferation and effector functions in peripheral T cells that then kill pathogeninfected cells and control immune responses. During T cell development in the thymus, somatic mutation of the antigenbinding TCR a/b subunit genes creates a thymocyte repertoire with random TCR specificities. However, many of these TCRs are non-functional or interact with the body's self-antigens with high affinity, causing autoimmune disorders if the respective T cells were allowed to mature. To prevent this, thymic selection processes eliminate thymocytes carrying TCRs that fail to interact with, or interact with too strong affinity with self-peptide-MHC (pMHC) complexes. The latter process is known as negative selection, a key mechanism of central tolerance. Only those thymocytes whose TCR generates mild signals are positively selected to mature into T cells, which then populate peripheral organs. Balanced positive and negative selections are critical for generating a diverse but self-tolerant T cell repertoire [8][9][10]. Recent experiments provided a more complex picture of thymic selection, where certain high affinity peptides can 'agonist select' distinct regulatory T cell types [11,12].
We and others identified ItpkB as essential for thymocyte positive selection [2,6,7]. ItpkB 2/2 DP thymocytes show intact proximal TCR signaling but defective IP 4 production, Itk PIP 3binding, signalosome recruitment and activation with ensuing reduced PLCc1 activation, DAG production, and, Ras/Erk activation [2]. The ability of soluble IP 4 to bind to the Itk PH domain and in low mM doses promote PIP 3 binding, and the ability of the Itk PH domain to oligomerize suggested that IP 4 might promote Itk recruitment to membrane-PIP 3 through a cooperative-allosteric mechanism. In this model, IP 4 -binding to one PH domain in an oligomer allosterically increases the ligand affinities of the other PH domains in the same oligomer [2]. IP 4 promoted Itk activation appears to be required for sufficient Itk activation to ensure positive selection, because an exogenous DAG-analog restored positive selection of ItpkB 2/2 thymocytes [2]. However, high-dose IP 4 inhibited Itk PH domain binding to PIP 3 in vitro [2]. Whether it does so in vivo is unknown [14]. In neutrophils, NK cells and myeloid progenitors, IP 4 competitively limits Akt PH domain binding to membrane PIP 3 [18][19][20]. Which PH domains are positively versus negatively controlled by IP 4 , and what determines whether IP 4 promotes or inhibits PH domain binding to PIP 3 or leaves it unaffected are important open questions [14,21]. In particular, the Itk PH domain might be bimodally regulated by IP 4 . However, the detailed molecular interactions between Itk, PIP 3 and IP 4 in vivo are not well characterized. This leaves room for multiple alternate hypotheses/mechanisms. For example, one could also propose that the binding affinity of PIP 3 and IP 4 for Itk changes from a low to a fixed high value above a threshold IP 4 concentration. Such a mechanism implies that the interaction of Itk with IP 4 and PIP 3 after the threshold IP 4 concentration is reached does not involve a positive feedback. The situation is further confounded by elusive results from experiments probing Itk oligomerization [2,[22][23][24][25][26][27][28].
The current lack of a mechanistic understanding of how IP 4 controls Itk PIP 3 -interactions and whether Itk PH domain  oligomerization is physiologically relevant arises from difficulties  in quantitatively measuring the interactions between Itk, IP 4 and  PIP 3 , and in generating soluble Itk PH domain preparations for  biophysical studies and non-oligomerizing Itk PH domain mutants for genetic analyses. Additional limitations arise from difficulties in measuring membrane recruitment of Itk in cell population based assays. It is also difficult to measure PIP 3 bound Itk or phosphorylation of PLCc1, a substrate of PIP 3 bound Itk, in large numbers of individual cells using flow cytometry techniques due to limited antibody quality. In vitro and cell-based studies based on ectopic Itk expression suggest the existence of several different monomeric and oligomeric Itk species, including head-to-head and head-to-tail dimers [2,[22][23][24][25][26][27][28]. Andreotti and colleagues [22] showed that Itk molecules can self associate via their SH2-SH3 domains into auto-inhibitory oligomers. This is hindered by SLP-76 interactions with the Itk SH2-SH3 domains. It was suggested that Itk molecules might exist as auto-inhibited multimers in the cytosol, but after plasma membrane recruitment, Itk monomers might mediate downstream activation [22,26]. Other experiments [27,28] employing fluorescence complementation showed that formation of Itk head-to-head and head-to-tail dimers requires the PH domain and may primarily occur at the plasma membrane, although low-abundance cytoplasmic dimers have not been excluded. Here, monomeric Itk was proposed to be primarily cytoplasmic and autoinhibited [27]. At least head-to-head dimerization is unaffected by mutations in the other (SH2/SH3) domains [28]. We found that the Itk PH domain can oligomerize with other Itk PH domains or full length Itk [2]. Thus, the PH domain is well suited to contribute to at least certain modes of Itk oligomerization, some of which could have positive or a combination of positive and negative functions regulated by IP 4 / PIP 3 . This could account for the limited activity-enhancing effect of disrupting SH3/SH2-domain mediated Itk dimerization [26].
Altogether, whether Itk PH domain dimerization has a physiological function, whether it promotes or inhibits Itk activation, whether IP 4 controls Itk function through positive or negative feedback, or both, and whether IP 4 has additional unrelated functions in thymocytes, are all contentious questions in the field. Resolving them is very important, because PI3K is a paramount regulator of signaling from many receptors in most cells. PIP 3 hyperactivity is a major contributor to immune, metabolic and other diseases including cancers [29,30]. IP 3 3kinases are broadly expressed and IP 4 has been found in many cell types. Thus, IP 4 regulation of PIP 3 function could be broadly important and elucidating the precise molecular mechanisms through which IP 4 controls PIP 3 signaling improves our understanding of a very fundamental and important signaling pathway with great therapeutic relevance [14].
To further explore how the presence or absence of Itk PH domain oligomerization, of positive or negative IP 4 feedback or both, or of specific molecular modes of association of Itk, PIP 3 and IP 4 impact TCR signaling, we constructed seven different molecular models (Table 1 and Figure S1B). We used a Maximum Entropy (MaxEnt) [31][32][33] based approach to quantify the robustness of each model against variations in rate constants and protein expression levels at the single cell level. Each model was constrained to reproduce the Itk activation kinetics of an entire cell population measured in biochemical experiments. We found that those models involving dimeric Itk molecules with IP 4 mediated competing positive and negative feedbacks are most robust. As in many other cell signaling systems [34], the actual signaling kinetics in thymocytes are likely to be robust against such variations, while retaining their sensitivity to small variations in antigen affinity or dose. On this basis, our simulations best support biphasic Itk regulation by IP 4 in thymocytes. Future testing of this exciting hypothesis will require the so far unsuccessful generation of nonoligomerizing Itk PH domain mutants and their expression in Itk 2/2 mice, along with currently impossible single-cell measurements of IP 4 levels in large cell populations.

Results
Multiple Molecular Models can be Constructed to Probe Itk, IP 4 , and PIP 3 Interactions in silico We constructed seven different molecular models (Table 1, Figure S1B) based on available details about interactions between Itk, PIP 3 and IP 4 from the biochemical studies described above. Including Itk kinase domain activation by Lck only caused qualitative changes in the relative robustness of the models (Fig.  S17, Tables S9-S15). Therefore, for simplicity, we considered models that do not contain Itk activation by Lck explicitly. We also did not consider Itk autophosphorylation explicitly in the models as it does not affect Itk catalytic activity. In addition, the role of Itk autophosphorylation in PLCc1 activation remains unclear [22]. Since we aimed to elucidate general characteristics of the kinetics of PIP 3 binding to Itk, we used a simplified modeling scheme ( Fig. 1) and did not consider the detailed molecular composition of the TCR and the LAT associated signalsome. The models also do not investigate different mechanisms for formation of Itk oligomers. Rather, they probe the functional consequences of having Itk PH domain dimers versus monomers and how these can affect interactions between Itk, PIP 3 and IP 4 in the presence or absence of IP 4 mediated positive feedback. The kinetics of PIP 3 production due to signal-dependent recruitment of PI3K are not considered explicitly as PIP 3 is produced at a much faster time scale (in seconds, [35] [36] [37]) than the time scales of PLCc1 activation (up to 60 min, Figure S18). The concentrations of LAT bound Itk and of PIP 3 were considered approximate markers for the strength of the stimulation generated by antigen-TCR interactions. Therefore, we considered fixed initial concentrations of Itk and PIP 3 in the models. We approximated the production of IP 4 from PIP 2 by a single one-step reaction to simplify the models further.
The models can be broadly classified into two types: (i) Models M1-M4 and M7 containing IP 4 mediated positive feedbacks. (ii) Models M5 and M6 lacking IP 4 positive feedback. In each type, we further considered models that contained Itk dimers (models M1-M3, M5, M7), or monomers (models M4, M6). In models M1-M3, each of the two PH domains in the Itk dimer can independently bind to either IP 4 or PIP 3 with a weak affinity when the other PH domain is unoccupied. However, once a PH domain is bound to an IP 4 molecule, it allosterically increases the affinity of the other PH domain for PIP 3  In Silico Modeling of Itk Activation Kinetics PLOS ONE | www.plosone.org IP 4 concentration is generated. We assumed that the small threshold level of IP 4 is generated at a time scale much smaller than the timescale (min) of robust Itk activation and did not consider the kinetics generating the threshold level of IP 4 explicitly in M5 and M6. The models are summarized in Table 1, Figure  S1B, and Tables S1-S8.

The Shape of Transient Itk Activation Kinetics Depends on Specific Molecular Wirings and Feedbacks in the Different Models
We studied the kinetics of Itk binding to PIP 3 using deterministic mass-action kinetic rate equations described by ordinary differential equations (ODE) for all the models, ignoring stochastic fluctuations in the copy numbers of signaling proteins occurring due to the intrinsic random nature of biochemical reactions [38]. Including such fluctuations did not change the kinetics qualitatively ( Figures S2-S3). In all seven models, the kinetics of PIP 3 bound Itk showed a transient behavior ( Fig. 2A); PIP 3 bound Itk started with a low concentration, reached a peak value at an intermediate time, and then fell back to a small concentration at later times. We found that initially few Itk molecules were bound to PIP 3 . With increasing time, more Itk molecules became associated with PIP 3 molecules due to the binding reactions between Itk and PIP 3 . This produced the rise in the Itk-PIP 3 concentration. However, as the concentrations of PIP 3 bound Itk molecules increased, they also induced increased production of IP 4 molecules. IP 4 competed with PIP 3 for binding to the Itk PH domain, and when the number of IP 4 molecules exceeded that of PIP 3 molecules, most of the Itk molecules were sequestered to the cytosol by forming stable complexes only with IP 4 . This reduced the rate of PIP 3 association of Itk and eventually resulted in the decrease of the PIP 3 bound Itk molecules. IP 4 outnumbered PIP 3 at later times because the number of PIP 2 molecules, the source of IP 3 and IP 4 in a cell, is considered not limiting in contrast to PIP 3 [37,39]. We emphasize that the results of our models do not depend on the cytosolic nature of Itk-IP 4 complexes, but on the model assumption that Itk (or Itk oligomers) Figure 1. Relevant basic interactions between Itk, PIP 3 and IP 4 . Following TCR-pMHC binding, Itk molecules are bound by the LAT signalosome via SLP-76 (not shown). Itk molecules (monomers or dimers, blue diamonds), bind the membrane lipid PIP 3 with low affinity through their PH domains. PIP 3 bound Itk phosphorylates and thereby activates LAT-bound PLCc1. Activated PLCc1 then hydrolyzes the membrane lipid PIP 2 into the soluble second messenger IP 3 , a key mediator of Ca 2+ mobilization. IP 3 3-kinase B (ItpkB) converts IP 3 into IP 4 (red filled circle). For our in silico models, we simplified this series of reactions, encircled by the orange oval, into a single second order reaction where PIP 3 bound Itk converts PIP 2 into IP 4 . In models M1-M4 and M7, IP 4 modifies the Itk PH domain (denoted as Itk C , purple diamonds) to promote PIP 3 and IP 4 binding to the Itk PH domain. At the onset of the signaling, when the concentration of IP 4 is smaller than that of PIP 3 , IP 4 helps Itk C to bind to PIP 3 (left lower panel). However, as the concentration of IP 4 is increased at later times, IP 4 outcompetes PIP 3 for binding to Itk C and sequesters Itk C to the cytosol (right lower panel). In models M5/M6, IP 4 and PIP 3 do not augment each other's binding to Itk. However, IP 4 still outcompetes PIP 3 for Itk PH domain binding when the number of IP 4 molecules becomes much larger than that of PIP 3 molecules at later times. doi:10.1371/journal.pone.0073937.g001 In Silico Modeling of Itk Activation Kinetics PLOS ONE | www.plosone.org bound to IP 4 at every PH domain does not induce any PLCc1 activation.
We characterized the 'shape' of the temporal profile of PIP 3 bound Itk in terms of (i) the largest concentration of PIP 3 bound Itk in the entire temporal profile (peak amplitude value A); (ii) the time taken for PIP 3 bound Itk to reach the peak value (peak time t p ); and (iii) the time interval during which the PIP 3 bound Itk concentration is greater than or equal to half of the peak value (peak duration t w , Fig. 2B). A dimensionless variable quantifying the asymmetry in the shape of the kinetics, denoted as the asymmetry ratio R = t w /t p (Fig. 2B), turned out to be a useful indicator for differentiating temporal profiles of concentrations of PIP 3 bound Itk in simulations and experiments. R also quantifies if the time scale for the decay of the concentration of PIP 3 bound Itk after the peak value is reached is larger than or comparable to t p (the timescale for producing the peak value A). E.g.,when R > 1, it implies that the t p is comparable to decay time. R . 1 indicates a more persistent signal with long decay times. Differences (transient vs. persistent) in the shapes of kinetic profiles of signals downstream of Itk activation have been observed to influence thymocyte decision outcomes [2,40]. Therefore, R, which characterizes the persistent or transient nature of Itk activation, also contains details directly relevant for thymic selection outcomes. We found that the shape of the transient kinetics of PIP 3 binding to Itk varied substantially depending on the feedbacks and the molecular wiring of the networks. Since the reaction rates used in the models are difficult to measure in vivo for thymocytes, we estimated the rates based on interaction strengths measured in vitro between PH domains and inositol phosphates in other cells, and from temporal profiles of PLCc1 activation measured in experiments with T cells (Tables S1-S7). Previous work demonstrated the essential role of phosphorylated PLCc1 and its kinetics in regulating thymocyte positive, negative and agonist selection [41,42]. Phospho-PLCc1 is also known to mirror other indicators of T cell activation such as TCRf-or The shapes of the temporal profiles can be characterized by the parameters peak time (t p ), peak width (t w ), and peak value or amplitude (A). The dimensionless asymmetry ratio R = t w /t p quantifies how symmetric the shape of the time profile is. A larger R value indicates larger asymmetry. (C) Variations in R in models M1-M7 for different initial concentrations of Itk and PIP 3 . Color scales for R values are shown on the right of each panel. doi:10.1371/journal.pone.0073937.g002 LAT-phosphorylation, or Erk-activation [2,40]. Therefore phospho-PLCc1 is a relevant marker for functional T cell responses.
We studied variations in the kinetics of PIP 3 bound Itk for different initial concentrations of Itk and PIP 3 . This probed how different ligand doses or affinities affected the PIP 3 binding of Itk. We found that the peak concentration of PIP 3 bound Itk increased in a graded manner with increasing initial Itk and PIP 3 concentrations in all models ( Figure S4). However, the peak time t p ( Figure S5), and the asymmetry ratio R (Fig. 2C), were affected differently in different models. Among the feedback models, M1-M3 and M7 containing Itk dimers generated smaller values (varied between 2 to 6) of R compared to monomeric model M4 which produced a much larger range of R (,20 -120) (Fig. 2C). The models lacking positive feedbacks (M5 and M6) generated large values (,100-700) of R compared to feedback models with Itk dimers (Fig. 2C). In the feedback models, the initial low affinity binding-unbinding interactions between Itk and PIP 3 /IP 4 are converted into high affinity interactions due to the positive feedback. Therefore, a large part of t p is spent in building up the positive feedback interactions controlled primarily by the weak affinity binding-unbinding rates (or K D ). The small values of R in models M1-M3 and M7 occurred because stronger negative feedbacks resulted in much smaller timescales for substantially reducing the concentration of PIP 3 bound Itk after it reached its peak value compared to the other models. In the models lacking positive feedback (M5, M6), concentrations of PIP 3 bound Itk decreased at a much slower rate than the peak time leading to large values of R. In the monomeric model, the relatively weaker strength of positive and negative feedbacks resulted in larger decay time scales for the PIP 3 bound Itk, producing large values of R. These results are analyzed in detail in the web supplement and Figures S6-S11. We will show below how the ability of feedback models with Itk dimers to produce R values within a small range leads to higher robustness of these models against parameter variations at the single cell level.

Models Containing Dimeric Itk and IP 4 Mediated Dueling Positive and Negative Feedbacks are the Most Robust Models
Quantification of robustness in in silico models. The reaction rates describing non-covalent primary and secondary interactions between Itk, PIP 3 and IP 4 can depend on specific properties of the local cellular environment, such as local membrane curvature [43], molecular crowding [44,45], and the presence of different lipid molecules in the proximity [46]. Since these factors can vary from cell to cell, the reaction rates can vary at the single cell level. In addition, protein expression levels can vary between cells. Such variations are also known as extrinsic noise fluctuations [47,48]. The IP 4 production rate depends on the concentrations of ItpkB, Calmodulin (CaM), and released calcium [3]. Hence, the IP 4 production rates in our models which approximate all such dependencies with a one-step reaction will vary between individual cells as well. The above variations are capable of producing differences in the shapes of temporal profiles of activation of signaling proteins in individual cells [49]. In the coarse-grained or approximate models we have constructed, many molecular details have been approximated. For example, multiple phosphorylation sites or SH2/SH3 binding sites of Itk, LAT, SLP-76 and their regulation via TCR induced signaling are not considered explicitly [3,22,50,51]. These detailed molecular signaling events can depend on the concentrations of proteins, enzymes, and lipids, and can thus be regulated differently in different cells due to extrinsic noise fluctuations. Consequently, the rates in our in silico models that effectively describe those detailed signaling events can vary from cell to cell. Consistent with this view, our simulations with the ODE models showed that the shape of the kinetics of PIP 3 bound Itk, characterized by, A, t p , and R, changed significantly as the rate constants and initial concentrations in a model were varied (Figures S12-S14). Thus, activation kinetics of a marker molecule (e.g. PLCc1) measured in experiments (e.g., immunoblots) assaying a large cell population represent averages over a range of temporal profiles with different shapes occurring at the single cell level.
We found that for some ranges of the reaction rates, multiple different in silico models can produce the same values of A, t p , and R ( Figures S12-S14). This implied that more than one in silico model could reproduce the mean temporal profile measured in cell population assays. However, it is possible that each model could show a different degree of robustness to variations in reaction rates and initial concentrations at the single cell level. Robustness of time dependent responses in a cell population against variations at the single cell levels has been observed in several systems, e.g., oscillations in adenosine 39,59 cyclic monophosphate (cAMP) concentrations in a population of Dyctostelium [52,53], or damped oscillations of protein 53 (p53) in a population of human breast cancer epithelial cells [54]. Robustness of cellular functions against variations in external conditions and cell-to-cell variability has been proposed as a required design principle for a wide range of biochemical networks [55][56][57][58]. We therefore decided to ask: Which model(s) can accommodate the largest variation in reaction rates and initial concentrations, while reproducing the mean temporal profile of PIP 3 bound Itk measured as generation of phosphorylated PLCc1 in cell population experiments? We postulate that the answer to this question will point us to the molecular circuitry most likely to be the relevant model, in the sense that it robustly produces a specific temporal response at the cell population level despite variations in the kinetics in individual cells.
To identify the most robust model(s), we quantified robustness using a method based on the principle of Maximum Entropy (MaxEnt) [31][32][33]. MaxEnt provides a mechanism for estimating the probability distribution of the rate constants and initial Itk and PIP 3 concentrations under constraints derived from experimental data ( Fig. 3A-B, 4A-B). Here, we used the experimentally obtained values t p expt and R expt as the constraints. It is difficult to directly relate the amplitude (in units of number of molecules in the simulation box) in the in silico models to experiments, where amplitudes are calculated from the fold change of the immunoblot intensities upon stimulation. Therefore, the experimental values of A can be related to the number of activated molecules, at best, through a proportionality constant dependent on specific protocols used in an assay. Because of these issues we chose a value of A expt , representing A in experiments, where every in silico model produced amplitudes at A expt for a set of parameters within the range of variations considered here. We then varied A expt to investigate the change in robustness of the models and address the arbitrariness in the choice of A expt . We constructed a relative entropy measure (Kullback-Leibler distance, D KL , calculated on the log 10 scale) [59] that measures the deviation of the constrained MaxEnt distribution from the unconstrained MaxEnt distribution, in which all values of the rate constants and initial concentrations are equally likely (uniform distribution). Thus D KL is being used as a measure of how ''close'' each model can get to one which is completely indifferent to the values of the rate constants, given the experimental constraints. We then compared D KL across our models in order to find the most robust model compatible with experimental results. Note that the minimum value of D KL is 0, with smaller values indicating greater robustness. We have also analyzed D KL for different models when t p expt and R expt were constrained but the amplitude A expt was not constrained ( Figure  S23). The results are qualitatively similar to that of the case when t p expt , R expt , and A expt were constrained. This indicates that the robustness of the temporal shape of Itk membrane recruitment kinetics rather than the amplitude contributes substantially toward the increased robustness of the feedback models with Itk dimers.
Experimental analysis of Itk activation kinetics in mouse thymocytes. To determine which Itk activation profile predicted by models M1-M7 produces maximum robustness while reproducing experimental data, we analyzed Itk activation kinetics in mouse CD4 + CD8 + double-positive (DP) thymocytes, the developmental stage where positive selection occurs [9,60]. To generate a homogeneous cell population in which every cell expresses the same TCR and in which the TCR has not been stimulated by endogenous ligands prior to in vitro stimulation, we used OT1 TCR-transgenic, RAG2 2/2 , MHCI(b2m) 2/2 mice. Their DP cells express exclusively the transgenic OT1 TCR, which recognizes the ovalbumin-derived peptide ligand OVA and recently identified endogenous peptide ligands presented by MHCI molecules [40,61]. In MHCI 2/2 mice, no endogenous ligands are presented to OT1 TCR-transgenic T cells and their development is blocked at the DP stage due to impaired positive selection. In vitro, OT1 TCR transgenic DP cells can be stimulated with MHCI tetramers loaded with OVA peptide [40,61]. Due to its high affinity for the OT1 TCR, OVA stimulation generates strong TCR signals and induces DP cell deletion. A number of OVA-derived altered peptide ligands (APL) have been generated which carry single or multiple amino acid substitutions compared to OVA. In the peptide series OVA.Q4R7.Q4H7.G4, such substitutions progressively reduce OT1 TCR affinity and signaling capacity [40]. Consequently, OVA and Q4R7 cause OT1 TCRtransgenic DP cell negative selection, whereas Q4H7 and G4 trigger positive selection.
We used MHCI tetramers presenting either one of these four peptides to stimulate RAG2 2/2 MHC 2/2 OT1 TCR-transgenic DP cells for various time points. We analyzed PLCc1 phosphorylation at Y 783 , normalized to total PLCc1 protein levels, as a measure for Itk activation [2] (Fig. 3A). Stimulation by all peptides induced fast PLCc1 phosphorylation already at 1 min which peaked at 2 min and then decreased over the next 60 min to very low levels which, however, were still above background levels in unstimulated cells. The decrease was fastest between 2 and 5 min and then progressively slowed down. As expected, overall levels of PLCc1 phosphorylation progressively decreased with decreasing peptide affinity/signaling capacity in the order OVA.Q4R7. Q4H7.G4. An asymmetric peak shape with an extended right flank was preserved across all signal intensities. We calculated the peak durations (t w ), peak times (t p ) and asymmetry ratios R = t w /t p in Table 2 for stimulation with OVA, Q4R7, Q4H7 and G4, respectively. Consistent with preserved peak asymmetry, all ratios R were .1.
Comparison between experiments and conclusions. The phospho-PLCc1 levels (representing active Itk) for different affinity peptides peaked at t p = 2 mins with R values from 1.9-4.3 (Table 2). Therefore, we fixed t p expt = t p = 2 mins (the bar indicates average over the cell population) for quantifying robustness in the in silico modeling. The low, medium and large initial Itk and PIP 3 concentrations represent stimulation by weak (G4), moderate (Q4R7, Q4H7) and high affinity (OVA) ligands, respectively. Analyzing D KL (Fig. 3C) showed that for large initial PIP 3 and Itk concentrations (representing OVA stimulation) the feedback models incorporating Itk PH domain dimers (M1-M3, M7) were substantially more robust (Smaller D KL values) than the models lacking feedbacks (M5, M6) for small values of R (,3). Monomeric feedback model M4 produced large D KL values (1.5-5). M5, M6 and M4 produced much larger ranges of R ( Figure  S14) as the parameters were varied compared to the feedback models with Itk dimers where the values of R were clustered around R expt ,2. This behavior contributed substantially to the increased robustness of the feedback models with dimers as these models could accommodate for larger ranges of parameter variations while being able to maintain the constraint imposed by R expt . The relative robustness of the feedback versus feedbackfree models showed similar qualitative trends for the other ligands, Q4R7, Q4H7, and G4 ( Figure S15-S16). This suggests that the models containing feedbacks and Itk dimers are substantially more robust than models with Itk monomers or lacking feedbacks.
Evaluation of robustness in polyclonal thymocytes. The molecular wiring of Itk, PIP 3 and IP 4 interactions is unlikely to depend on the clonal nature of the T cells. Thus, the feedback models with Itk dimers should also be more robust than the other models when used to describe the kinetics of PLCc1 activation in polyclonal DP thymocytes expressing many different TCRs with different ligand specificities, stimulated by antibodies against the common TCR subunit CD3 alone or with co-ligation of the common coreceptor CD4. Stimulation of non TCR-transgenic MHC 2/2 DP cells with 1 mg/ml or 5 mg/ml of aCD3 or combined aCD3/aCD4 antibodies produced different R expt and t p expt values than the OT1 system above ( Fig. 4A-4B, Figure S18, Table. S16). Calculation of the robustness constrained by R expt , t p expt and A showed that feedback models M1, M2, M3 and M7 are again substantially more robust than the other models ( Fig. 4C-4F, Figure S19, S20). Large variations of R in M4, M5 and M6 as parameters were varied again made these models substantially less robust than the feedback models with Itk dimers.

Discussion
Here, we used in silico simulations combined with a novel Maximum Entropy (MaxEnt) based method and cell population averaged measurements of PLCc1 activation kinetics to distinguish between multiple models constructed to elucidate different mechanisms of Itk activation in TCR signaling. Our analysis quantified the robustness of seven different models employing monomeric or dimeric Itk PH domains with or without positive and negative IP 4 feedback against variations of parameters (rates and concentrations) at the single cell level. MaxEnt has been widely used in diverse disciplines ranging from physics [62] via information theory [63] to biology [64][65][66][67] to estimate probability distributions of variables subject to constraints imposed by experimental data [33,65]. However, to our knowledge these methods have not been used for evaluating the robustness of dynamic models in cell signaling or gene regulatory systems. Using thymocyte positive selection as a physiologically important model process, our results show the usefulness of MaxEnt methods for such studies. We are currently working on extending the methods to include additional information from experiments (such as variances), and also evaluating their performance in comparison with closely related approaches such as Bayesian analysis [68].
Our simulations predict that the models containing IP 4 feedbacks and Itk dimers are most robust. This is consistent with our previously proposed model of cooperative-allosteric regulation of Itk-PIP 3 interactions via IP 4 -binding to oligomeric Itk PH domains [2]. Thymocyte selection critically depends on TCR induced signals. Small differences in antigen peptide concentration or affinities for the same TCR can produce opposite (negative vs. positive) selection outcomes [40]. Thus, we consider it plausible that for a fixed antigen dose and affinity (or average initial concentrations of Itk and PIP 3 in our models), TCR signaling in thymocytes should be robust against cell-to-cell variations of protein/lipid concentrations, rate constants and local environment. But TCR signaling should retain sensitivity to small variations in antigen affinity or dose. A direct experimental validation of this assumption will require to test the probability distributions of t p , R, and A in cell populations where PLCc1 activation kinetics are measured in individual cells. However, we were unable to perform such single cell comparisons due to the insensitivity of FACS-based PLCc1 signaling assays. This indicates the importance of studying the effects of network architecture, rate constants, protein and lipid concentrations on system robustness in DP thymocyte selection in detail in the future. Thymocytes are an excellent in vivo model to probe the exquisite dependency of cell fate decisions on the affinity of TCR ligands with important physiological and pathological implications. This provides a valuable addition to the experimental and theoretical investigations of robustness in synthetic systems or transformed tissue culture cells in vitro.
On the basis of robustness, our simulations support bimodal positive and negative Itk regulation by IP 4 in thymocytes. They make a supportive argument that Itk PH domain oligomerization and IP 4 feedback are physiologically important, consistent with the severely defective TCR signaling, IP 4 production, Itk/PLCc1 activation, positive selection and resulting immunodeficiency in ItpkB 2/2 mice, the ability of IP 4 to bimodally control Itk PH domain binding to PIP 3 in vitro, and the reported Itk PH domain oligomerization [2,6,7,28]. They do, however, not exclude the possibility that IP 4 also has additional, unknown functions in DP cells [14].
Testing this exciting hypothesis will require currently impossible single-cell measurements of IP 4 levels in large cell populations. Moreover, the physiological roles and modes of Itk oligomerization, the specific PH domain contributions to Itk oligomerization, whether Itk oligomerization occurs in the cytoplasm or at the plasma membrane or both, whether it exclusively inhibits or can also promote Itk activation, and whether IP 4 promotes or inhibits Itk PH domain binding to PIP 3 or does both depending on its local concentration are all matters of active debate [2,[22][23][24][25][26][27][28]. Their conclusive elucidation requires quantitative biophysical studies of full length Itk with or without mutational perturbation of individual and combined interactions among the different Itk domains implicated in its monomeric and oligomeric selfassociation, and the reconstitution of Itk 2/2 mice with these mutants at endogenous expression levels. Unfortunately, difficulties to produce sufficient quantities of soluble full-length Itk or Itk PH domain protein, and a tendency of Itk and its PH domain to aggregate in vitro have precluded more quantitative analyses of Itk PH domain oligomerization and IP 4 /PIP 3 interactions, as well as the generation of non-oligomerizing Itk PH domain mutants. Despite progress regarding SH2/SH3/proline-rich domain interactions [22][23][24][25][26] and some evidence for PH domain involvement [2,27,28], formation of several different homotypic Itk dimers with differing subcellular localization and functions further complicates such analyses and their interpretation. Our in silico results suggest that by enabling competing positive and negative IP 4 induced feedback, Itk PH domain oligomerization could render Itk signaling in DP thymocytes much more robust to parameter fluctuation between individual cells than could be achieved without Itk dimers, or without IP 4 feedback. Models M1-M3 and M7 involving Itk dimers and IP 4 feedbacks showed substantially larger robustness than models lacking feedbacks (M5-M6) or containing only monomeric Itk (M4). M1-M3 and M7 can describe the experimentally observed PLCc1 kinetics with similar robustness. They differ only at the level of secondary Itk/ IP 4 /PIP 3 interactions. Similar robustness and the inherent variability of experimental data preclude the identification of one of these dimeric Itk feedback models as the only one operative in vivo thus far.

Signaling Kinetics in the in silico Models
We constructed ODE based models. The ODEs described kinetics of concentrations of proteins and lipids in two well-mixed compartments representing plasma membrane and cytosol ( Figure  S1A). The biochemical signaling reactions for each model are shown in Tables S1-S7. The details regarding the construction of the ODEs and the parameters are given in the web supplement and Figure S1. We use the rule based modeling software package BioNetGen [69] to generate time courses for the species kinetics for the signaling networks described by models M1-M7. This program produces a set of ODEs corresponding to the mass-action kinetics describing biochemical reactions in the networks and solves them numerically using the CVODE solver [70]. The ODEs for each model are listed in the supplementary material.

Quantification of Robustness Based on the Maximum Entropy Principle
When a variable x can assume multiple values and is distributed according to a probability distribution p(x), then the uncertainty associated with the distribution can be quantified by the entropy (S) defined as, S is non-negative and is maximized when x is distributed according to a uniform distribution (i.e., x can take any value within a range with equal probability). Suppose p(x) is unknown, but we do know the average value of a variable, f, that is a function of x, i.e., f = f(x). We can then maximize S under the constraint The constrained MaxEnt distribution is given by p(x) / exp(2lx), where the constant l, also knovn as the Lagrange multiplier, is determined by solving Eq. (2) for l when the above MaxEnt distribution for p(x) is used in Eq. (2). The method can be Table 2. Values of peak time, peak width, and asymmetry ratio R calculated from the PLCc1 activation kinetics in Fig. 3 for different ligands.  easily generalized to accommodate multiple variables and constraints. We used the constraints imposed by t p expt , R expt , and A expt , or, t p expt and R expt that are measured over a cell population. Therefore, the MaxEnt distribution of the parameters in our calculation is given by, p({k i }) / exp(2l 1 t p ({k i }) -l 2 R({k i }) 2 l 3 A({k i })), where l 1 , l 2 and l 3 denote the Lagrange's multipliers, and {k i } denote the values of rate constants and initial concentrations in individual cells. The Lagrange multipliers can be calculated from the constraint equations, The MaxEnt distribution thus describes how t p , R, and A, in individual cells are distributed over a cell population. The distribution also produces an estimation of the probability distributions for the rate constants and initial concentrations that regulate t p , R, and A, through the functions t p ({k i }), R({k i }), and, A({k i }), respectively. The specific relationship between the parameters, {k i }, and the observables (t p , R, and A) is dependent on the molecular details of the models, M1-M7. In all the models prior to the MaxEnt calculation, the rate constants were chosen from a uniform distribution with lower and upper bounds equal to 1/10 and 10 times, respectively, of the base values shown in Tables S1-S7. Similarly, the initial concentrations of proteins (e.g., Itk) and lipids (such as PIP 3 ) were varied within a 35% [71] range from uniform distributions centered at the base values shown in Table S8. The joint uniform distribution in the parameters is given by q({k i }). We then used these MaxEnt distributions to quantify relative robustness of the models by calculating the Kullback-Leibler distance [59] D KL~X That is, for each model, we first find the probability distribution for the rate constants and initial concentrations that maximizes the entropy (robustness) for that model under the experimental constraints, giving the model a kind of ''maximum benefit of the doubt.'' We then compare the resulting MaxEnt models with one another to evaluate their relative robustness to variation in the rate constants, in order to select the model(s) most likely to correctly represent the actual kinetics. When p({k i }) is equal to q({k i }), D KL assumes the minimum value 0; as the distribution p({k i }) starts deviating from the uniform distribution, say by becoming sharply peaked around a particular value, D KL increases. Thus maximizing the entropy S, is equivalent to minimizing D KL in Eq. (4). The calculations of D KL were done at a specific antigen dose which fixed the average values of initial concentrations of Itk and PIP 3 . Therefore, the robustness calculations did not exclude the sensitivity of PLCc1 activation to changes in PIP 3 concentrations resulting from antigen dose variations. We calculated p({k i }) by by minimizing the D KL subject to the constraints imposed by Eq. (3). We used D KL to rank order the models for a particular measured value of t p expt , R expt , and, A avg . All the calculations were carried out using MATLAB. Additional details can be found in the supplementary material ( Figures S12-S15). Note that D KL is unaffected by inclusion of additional parameters that do not influence the experimentally measured variables ( Figure S21, Table S17). Thus having extra variables in a model does not in and of itself affect the relative robustness of models with variable numbers of parameters. We have used 100,000 sample points, which we have shown to be statistically sufficient in Figure S22 for the faithful calculation of D KL .

Thymocyte Stimulation and Immunoblot Analysis
All mice were housed in The Scripps Research Institute specific pathogen-free vivarium monitored by The Scripps Research Institute Department of Animal Resources. All animal studies were approved by The Scripps Research Institute IACUC and conform to all relevant regulatory standards.
DP cells were prepared as in [2] and rested at 37uC for 3 hours. Then, 10 7 DP cells per sample were incubated on ice for 15 min with 2.4 mM MHCI tetramers pre-loaded with either one of the altered peptide ligands OVA, Q4R7, Q4H7 or G4 [40], stimulated by rapidly adding 37uC warm PBS for the indicated times and quickly lysed in 100 mM Tris, pH 7.5, 600 mM NaCl, 240 mM n-octyl-b-D-glucoside, 4% Triton, 4 mM EDTA and a protease/phosphatase inhibitor cocktail (Roche). Lysates were cleared by centrifugation at 14000 rpm for 10 minutes at 4uC, resolved by SDS-PAGE and analyzed via immunoblot as previously described [2]. Band intensities were quantified via densitometry using NIH ImageJ software, and phosphoY 783 -PLCc1 intensities normalized to total PLCc1 amounts.  (Table S3) using the Gillespie algorithm. The curve in red is the solution of the mass action kinetics given by a set of ODEs. We use the same kinetic rates and initial concentrations for the stochastic simulations and the ODEs. (TIFF) Figure S3 Comparison between the ODE solutions and the stochastic trajectories averaged over a small number of cells. We compared the temporal profiles of concentrations of PIP 3 bound Itk obtained in simulations including stochastic copy number variations due to intrinsic noise fluctuations (red) with the solutions of the deterministic mass action reaction kinetics that ignored such fluctuations (solid black lines). The stochastic simulations were carried out by using Gillespie's method which provided exact numerical solution of the Master equations associated with the models. We used the same rate constants and initial concentrations for the stochastic simulations and ODE solutions. The kinetic trajectories were averaged over 500 realizations (or in silico ''cells'') for the stochastic simulations.  Figure S19 Checkerboard plot of the most robust models as R avg and A avg are varied for different doses of anti-CD3 and anti-CD3/CD4 antibodies. a) Itk 0 = 40 and PIP 3 0 = 130 molecules are used to emulate the 1 mg/mL anti CD3 stimulation. The t avg is held at 1 mins. The checkerboard diagram of the most robust models is shown as R avg and A avg are varied. b) Same as plot a) but Itk 0 = 100 and PIP 3 0 = 370 molecules are used as the initial concentrations. c) Itk 0 = 100 and PIP 3 0 = 370 molecules are used to emulate the 1 mg/mL anti CD3/CD4 stimulation. The t avg is held at 5 mins. The checkerboard diagram of the most robust models is shown as R avg and A avg are varied. d) Itk 0 = 140 and PIP 3 0 = 530 molecules are used to emulate the 5 mg/mL anti CD3/CD4 stimulation. The t avg is held at 1 mins. The checkerboard diagram of the most robust models is shown as R avg and A avg are varied. (TIFF) Figure S20 The plot of D KL for all the 7 models for a specific amplitude and different initial conditions for different doses of anti CD3 or anti CD3/CD4 antibodies. a) Itk 0 = 40 and PIP 3 0 = 130 molecules are used to emulate the 1 mg/mL anti CD3 stimulation. The t avg is held at 1 mins. The D KL is shown for an A avg = 16 molecules. b) Same as plot a) but Itk 0 = 100 and PIP 3 0 = 370 molecules are used as the initial concentrations and A avg = 60 molecules. c) Itk 0 = 100 and PIP 3 0 = 370 molecules are used to emulate the 1 mg/mL anti CD3/CD4 stimulation. The t avg is held at 5 mins. A avg = 60 molecules. d) Itk 0 = 140 and PIP 3 0 = 530 molecules are used to emulate the 5 mg/mL anti CD3/CD4 stimulation. The t avg is held at 1 mins and A avg is set equal to 80 molecules. The vertical orange bar shows the observed experimental values. (TIFF) Figure S21 Addition of parameters which weakly affect the Itk-PIP 3 kinetics, do not lead to any significant difference in the D KL . For Itk 0 = 100 and PIP 3 0 = 370, a) we have looked at the relative difference in the D KL of our old M3 (black) and M3 with the added reactions (magenta) for an amplitude average of 30 molecules and peak time average of 2 mins. b) We have looked at the relative difference in the D KL of our old M3 (black) and M3 with the added reactions (magenta) for an amplitude average of 40 molecules and peak time average of 2 mins. (TIFF)      S16 Values of peak time, peak width, and asymmetry ratio R calculated from the PLCc1 activation kinetics in Figure  S18. (DOCX)