Resolving Early Signaling Events in T-Cell Activation Leading to IL-2 and FOXP 3 Transcription

Signal intensity and feedback regulation are known to be major factors in the signaling events stemming from the T-cell receptor (TCR) and its various coreceptors, but the exact nature of these relationships remains in question. We present a mathematical model of the complex signaling network involved in T-cell activation with cross-talk between the Erk, calcium, PKCθ and mTOR signaling pathways. The model parameters are adjusted to fit new and published data on TCR trafficking, Zap70, calcium, Erk and IκBα signaling. The regulation of the early signaling events by phosphatases, CD45 and SHP1, and the TCR dynamics are critical to determining the behavior of the model. Additional model corroboration is provided through quantitative and qualitative agreement with experimental data collected under different stimulating and knockout conditions. The resulting model is analyzed to investigate how signal intensity and feedback regulation affect TCRand coreceptor-mediated signal transduction and their downstream transcriptional profiles to predict the outcome for a variety of stimulatory and knockdown experiments. Analysis of the model shows that: (1) SHP1 negative feedback is necessary for preventing hyperactivity Processes 2014, 2 868 in TCR signaling; (2) CD45 is required for TCR signaling, but also partially suppresses it at high expression levels; and (3) elevated FOXP3 and reduced IL-2 signaling, an expression profile often associated with T regulatory cells (Tregs), is observed when the system is subjected to weak TCR and CD28 costimulation or a severe reduction in CD45 activity.

Abstract: Signal intensity and feedback regulation are known to be major factors in the signaling events stemming from the T-cell receptor (TCR) and its various coreceptors, but the exact nature of these relationships remains in question.We present a mathematical model of the complex signaling network involved in T-cell activation with cross-talk between the Erk, calcium, PKCθ and mTOR signaling pathways.The model parameters are adjusted to fit new and published data on TCR trafficking, Zap70, calcium, Erk and IκBα signaling.The regulation of the early signaling events by phosphatases, CD45 and SHP1, and the TCR dynamics are critical to determining the behavior of the model.Additional model corroboration is provided through quantitative and qualitative agreement with experimental data collected under different stimulating and knockout conditions.The resulting model is analyzed to investigate how signal intensity and feedback regulation affect TCR-and coreceptor-mediated signal transduction and their downstream transcriptional profiles to predict the outcome for a variety of stimulatory and knockdown experiments.Analysis of the model shows that: (1) SHP1 negative feedback is necessary for preventing hyperactivity

Introduction
The actions of CD4+ T-cells are controlled by the stimulatory signals and cytokine milieu in their surrounding tissue environment [1,2].Extracellular signals naturally drive CD4+ T-cells to activate and assume one of many important immunological roles, each characterized by a distinct profile of signaling events and cytokine secretion.These signals are processed by the intracellular signaling pathways to alter the transcription profiles that direct T-cell activation and differentiation.To explore this process from a quantitative perspective, a mathematical model is constructed from the known molecular interactions and regulations.
A critical second signal driving T-cell activation and fate determination is achieved by incorporating CD28-mediated costimulation [27] and its downstream signaling events into the model.In spite of the ability of the TCR to upregulate factors necessary for full activation, TCR-antigen recognition ultimately leads to anergy or a state of hyporesponsiveness in the absence of a second signal [28].Alternatively, TCR signaling with additional costimulation of CD28 on the surface of T-cells greatly enhances IL-2 production, T-cell proliferation and the prevention of anergy.This synergistic effect is believed to be mediated through PKCθ and Akt (also known as protein kinase B), which recruit and regulate several signaling molecules, including NFκB [23,29] and mammalian target for rapamycin (mTOR) [30].The serine/threonine protein kinase mTOR is another important player in CD4+ T-cell activation.mTOR complex 1 (mTORC1) has many known functions, including the ability to control the expression of forkhead box P3 (FOXP3), a master transcriptional regulator in the development of T regulatory cells (Tregs) [30,31].While much less is known about mTORC2, it was recently shown to promote PKC, Akt and NFκB signaling and to regulate T helper type 1 (Th1) and Th2 commitment [32].
A number of mathematical models exist that span the entirety or portions of these T-cell signaling events.The models presented in [33,34] span these events and are able to capture many dominant behaviors.However, these discrete-time logical models do not support a detailed analysis of the interactions in a continuous state and parameter space and the ability to predict responses to analog knockdown or mutation scenarios.Other models are very detailed, but only consider portions of the system, such as the MAPK [35] and NFκB [36,37] pathways.Some of these models consider every hypothesized molecular interaction or complex formation.This level of detail is prohibitive to complete for the entire signaling network, predominantly due to computational tractability, but also due to uncertainties in the existence or rates of these events.A few models set out to address the role of location as a regulatory event, analyzing T-cell signaling events within a spatial context [38,39].These models tend to be limited in scope in order to offset the added computational expense of such an analysis.Our model differs from past approaches in that we consider a broad scope of early TCR and CD28 signaling leading to gene transcription in continuous time, consider only the major known phosphorylation, dephosphorylation, association and dissociation events and address location through compartments and not as a partial differential equation with diffusion.
Here, we present a mathematical model of T-cell activation that focuses on the early steps in TCR trafficking and signal amplification by kinases, phosphatases and other second messengers, as well as CD28 coreceptor signaling and gene transcription of IL-2 and FOXP3.Sections 2 and 3 describe the formulation of the model and its fitting to training data.Model simulations and corroborating experimental evidence are provided in Section 4, demonstrating the model's ability to reproduce many known experimental results.This section also describes predictions made by the model about the roles of phosphatases CD45 and SHP1 in T-cell activation.Finally, conclusions are summarized in Section 5.

Materials and Methods
This section provides an overview of the proposed T-cell activation model and details the methods and materials employed in its development.These include the computational methods and experimental datasets used for model parameter identification.

The Mathematical Model: Overview and Scope
The mathematical model reflects the current understanding of the intracellular signaling events leading to T-cell activation.The reaction scheme, as shown in Figure 1, encompasses the primary events that occur upon TCR engagement and CD28 costimulation.Appendix 1 describes each signaling element and assigns a symbol and initial value for use in constructing the mathematical model.Appendix 1 also lists the rate equations, parameter values and summarizes the reactions primarily in terms of protein binding, phosphate group transfers and kinase and phosphatase activation.The part of the model representing the TCR-mediated Erk pathway is based on the ordinary differential equation (ODE) model developed by Zheng [40] with updates in Perley et al. [41].The Zheng model takes into account the major events of early T-cell signaling leading to the activation of the protein tyrosine kinases Zap70 and SFK, the downstream Erk pathway, as well as negative feedback regulation by phosphatase SHP1.
The Zheng model was extended to specifically incorporate TCR trafficking, tuning of CD45 and SHP1 phosphatase activity, refinement of Erk signaling with AP1 activation and the addition of the calcium, NFκB and mTOR pathways leading to IL-2 and FOXP3 transcription.These changes are further detailed in Section 3. Molecular association and complex dissociation are designed to obey the law of mass action, i.e., assuming that the rate of a reaction is proportional to the product of reactant concentrations.Most enzyme-catalyzed reactions utilize second-order kinetics derived from a simplification of the Michaelis-Menten enzyme kinetics by assuming that the substrate concentration is limiting.Such a simplification contributes to reducing the complexity of the nonlinearities while retaining the structure of the pathway connectivity.For reactions catalyzed by enzymes not explicitly included in this model, the expression further reduces to first-order kinetics, where the enzyme concentration becomes a constant.In addition, single molecular transitions are also governed by first-order kinetics.Alternatively, enzyme-catalyzed reactions in which cooperative binding between substrates plays a significant role utilize the Hill equation.Scenarios characterized by these equations include the cooperative transcriptional regulation of IL-2 and FOXP3 and the activation of certain enzymes.The model is also designed to mimic the cell response to stimulation by αCD3, αCD28, phorbol 12-myristate 13-acetate (PMA) and the calcium ionophore ionomycin, as well as mimic perturbations in kinase and phosphatase activity by various small molecule inhibitors.The complete model consists of 48 ODEs with 154 reaction parameters and simulated in MATLAB using the ode15s solver.

Global Sensitivity Analysis
Due to the size and complexity of the model, it is necessary to reduce the dimension of the uncertain parameter space to facilitate an effective and computationally tractable analysis of the model.Parameter sensitivity analysis provides a powerful tool for analyzing mathematical models of complex biological systems and can be used to facilitate parameter estimation by identifying and ranking the various contributors of parametric uncertainty in the model.Parameters with the lowest sensitivity ranks, i.e., parameters whose uncertainty causes the least amount of variability in the output space, can be neglected during the model fitting process and fixed at their initial guesses.For the purpose of this study, a variance-based sensitivity analysis is used to determine how the model output changes over the global parameter space: where S i,j is the sensitivity index of the i − th output to the j − th parameter.The sensitivity indices are numerically computed using Sobol's method by the efficient sparse-grid-based algorithm proposed by Buzzard [42,43].Parameters having little effect on model outputs upon perturbation are determined to be insignificant and fixed at their nominal values.Significant parameters are retained for parameter identification.

Parameter Identification
The ability of the model simulated with a given parameter vector to reproduce available experimental data is quantified by the objective function: where e t i is the error at each sample time point, W t i is the inverse of the measurement error covariance matrix at time t i , y(t i ) and ŷ(t i ) are the simulated model outputs and the mean values of the available data at time t i , respectively, C(t i ) is a binary matrix that indicates which outputs are measured at time t i and N s is the total number of sampled time points.This cost function is computed in log space in order to compress the range of and smooth the function, making its approximation with sparse-grid interpolation more effective.The "1" in the cost function is added to avoid the singularity of log at zero.Specific parameter values that are capable of producing data-consistent model simulations are extracted directly from the literature whenever possible.For those not found in the literature, they are either derived using known conditions or constraints of the system or estimated from experimental data using a combination of manual and automated calibration.First, the starting guess values are roughly estimated by back-of-the-envelope calculations and manual tuning to within a reasonable level of accuracy.This step is used to establish our parameter search space and its bounds, defined to be a hypercube spanning the nominal parameter vector by an order of magnitude both above and below.Next, Equation ( 2) is evaluated at a number of samples from the parameter hypercube, generated using Latin hypercube sampling in log space.Finally, starting from the parameter vector with the lowest sampled cost, MATLAB's nonlinear constrained optimization solver, fmincon, is employed to identify the parameter vector that minimizes Equation (2) constrained within the hypercube.To increase the chance of identifying the global minimizer, multiple starting points were considered that span the region of low cost.This procedure of automated parameter identification is aided by the use of sparse-grid interpolation, a powerful tool for generating fast-evaluating approximations of mathematical models [44].

Experimental Datasets for Model Development
Experimental datasets for the purposes of model development are compiled from the published literature or generated from experiments performed in our lab specifically for this manuscript.For datasets compiled from the published literature, only those that are commensurate with our model in terms of cell type, stimuli, output species and time scale are considered.Table 1 summarizes the attributes of the compiled datasets used to train the model.These datasets consist of dynamics (i.e., time courses), dose response experiments (i.e., multiple doses of an input or stimulus) and input response experiments (i.e., single doses of different input or stimuli combinations).In order to compare the model and observed behaviors directly, simulations are generated to mimic the experimental conditions used in our experimental set-ups (see Appendices 3 and 4) or as indicated in the corresponding source articles.Model simulations are shown with the corresponding data and error bars whenever possible; however, for cases in which data could not be visualized on the same plot, the reader is referred to the original source.

Model Development
The following subsections describe the process of formulating the proposed T-cell activation model as shown in Figure 1 from the original ODE model presented by Zheng [40].These changes include the incorporation of TCR trafficking mechanisms, tuning of CD45 and SHP1 phosphatase activity, refinement of Erk signaling with AP1 activation and the addition of the calcium, NFκB and mTOR pathways leading to IL-2 and FOXP3 transcription.These changes are further detailed in the following subsections.

Modeling TCR Trafficking
The TCR expression level at the cell surface is the result of a dynamic equilibrium maintained by the membrane expression of newly synthesized TCR, receptor internalization, recycling to the cell surface and degradation [46,47].As depicted in Figure 1, the TCR complex is represented in the model by five states: free TCR (T CR f ), ligand-bound TCR (T CR b ), phosphorylated ligand-bound TCR (T CR p ), internalized TCR (T CR i ) and degraded TCR (T CR deg ).TCR trafficking is modeled by allowing TCR complexes from all surface states (T CR f , T CR b and T CR p ) to be internalized [38].Internalized TCR subunits can be degraded or recycled back to the unbound state.Newly synthesized TCR contributes to the free receptor state with a rate proportional to degraded TCR.
Early studies demonstrated that the TCR is a constitutively cycling receptor [46,47].Thus, at steady state, a certain amount of TCR is endocytosed, while at the same time, an equal amount of TCR is exocytosed.The associated model rate parameters are set by the following observations.Several studies seem to agree that the constitutive endocytic rate constant for the TCR in resting T-cells is ~0.01 min −1 , meaning that ~1% of the cell surface-expressed TCR is internalized each minute [46].
TCR ligation induces downregulation and degradation of the TCR in a dose-dependent manner [45,46].In theory, TCR downregulation can be accomplished by an increase in the endocytic rate constant, a decrease in the exocytic rate constant or a combination of both.However, most studies found that TCR downregulation is caused by an increase in the endocytic rate constant to ~0.038 min −1 after TCR triggering [46], whereas the exocytic rate was constant.In addition, TCR degradation was also found to accelerate after TCR ligation.The TCR subunits in non-stimulated Jurkat cells were observed to degrade with rate constants of ~0.0011 min −1 , resulting in a half-life of ~10.5 h.Triggering of the TCR by anti-TCR Abs resulted in a three-fold increase in the degradation rate constants to ~0.0033 min −1 , resulting in a half-life of ~3.5 h [45].This process was modeled by making the receptor internalization and degradation rate parameters functions related to T-cell activation.These mechanisms were mathematically encoded in the model using: Equation ( 3) ensures a continuous transition between the constitutive endocytosis or degradation rates (k min ) and the ligand-induced rates (k max ) using a Hill equation, where X is the induction substrate concentration, K is the substrate concentration of the transition midpoint and n is the Hill coefficient denoting positive (n > 1) or negative (n < 1) cooperativity.Since the TCR cycling pathway is dependent on the CD3γ di-leucine-based motif and is activated by PKC [46,47], for receptor cycling, we set k = k int , X = [P KCθ], n = 2, and K is set to 5% of total PKCθ with k min = 0.01 min −1 and k max = 0.038 min −1 .Since the activation of Lck and Zap70 leads to recruitment of Cbl and ubiquitination (Ub) of the CD3 and ζ chains to induce degradation of the TCR in the lysosomes [46], for receptor degradation, we set k and K is set to 5% of total SFK with k min = 0.0011 min −1 and k max = 0.0033 min −1 .
At equilibrium without any stimulus, most studies also seem to agree that the pool of recycling TCR was distributed with approximately 75%-85% at the cell surface and 15%-25% inside the cells.Assuming the minimum constitutive rates for internalization and degradation with 80% surface TCR and 20% inside the cell, of which half are being degraded, we can compute the constitutive exocytosis rate as k exo = 0.0789 min −1 and the constitutive synthesis rate as k synth = k deg min to achieve the observed dynamic equilibrium in resting T-cells.
The resulting simulated surface and degraded TCR dynamics are consistent with expectations at equilibrium and after receptor engagement with a ligand.The dose response curves in Figure 2A,B show good agreement with the behaviors reported by Menne et al. [47], Geisler et al. [46] and von Essen et al. [45].

Tuning the Roles of CD45 and SHP1
The protein tyrosine phosphatases, CD45 and SHP1, are central players in signal amplification following antigen recognition by the TCR and the activation of many downstream second messenger and signaling molecules.CD45 is a leukocyte-specific transmembrane glycoprotein and a receptor-like protein tyrosine phosphatase (PTP) [8].The SFK member Lck is the best characterized CD45 substrate in T-cells [7].Lck exists in dynamic equilibrium with three main sub-populations: (1) open and activated (SF Kact in the model); (2) open and not activated ("primed") (SF Kdp); (3) closed and not activated (SF K) [9].Phosphorylation of SFKs at the negative regulatory site (Y505 in Lck) by the kinase Csk results in an intramolecular interaction with the SH2 domain, creating a folded inactive conformation [48].Dephosphorylation at this site by CD45 opens up the molecule, creating a "primed" molecule.Clustering of these primed SFKs results in the transphosphorylation of the activation loop (Y394 in Lck), which displaces it from the catalytic site and creates an active kinase by allowing substrate access.Dephosphorylation at this site by CD45 or other PTPs, such as SHP1, downregulates SFK activity and returns them to the primed state [9].Thus, CD45 functions as both a positive and negative regulator of the T-cell antigen receptor and in setting the threshold of activation.
The dual role of CD45 as both positive and negative regulator of T-cell activation is modulated in part by the distribution and movement of CD45 and its substrates proximal to the receptor complex inside and outside lipid rafts.Lck inside lipid rafts has been reported to be hyperphosphorylated and less active when CD45 is excluded [9].The engineered inclusion of CD45 to the lipid domains also decreases TCR signaling, consistent with its ability to downregulate signaling by dephosphorylating the positive regulatory site at later time points [9].These observations show that the localization of both CD45 and the SFKs can affect the phosphorylation state of SFKs.It is believed that upon receptor-mediated clustering, the SFKs are activated and often relocalize to lipid rafts, where the concentration of CD45 is much lower and the kinases can experience sustained signaling.The later recruitment of CD45 to these domains will then primarily dephosphorylate the activation site and downregulate activity [9].
The model is augmented to include both the positive and negative regulatory roles of CD45, as well as spatial localization effects, on T-cell activation with two explicit states: CD45 p and CD45 n (see Figure 3).The first, CD45 p , reflects its ability to activate substrates immediately following TCR triggering and promote signaling.Initially, this state exists at a high level with its major role being to dephosphorylate negative regulatory sites and "prime" SFKs (R 2 ).Following TCR engagement, the receptor cluster formation (represented in simulation by ligand-bound TCR) separates CD45 and its substrate.This causes CD45 p to enter an ineffective state, an implicitly modeled state called CD45 (R 23 ), and decreases its ability to activate its substrate.The second state, CD45 n , represents the negative regulatory function of CD45 on SFK.Initially, this state is minimally active with most existing in the inactive state, CD45.As TCR signaling progresses (represented in simulation by activated SFK), CD45 is recruited back to the lipid rafts (R 23a ), where its ability to dephosphorylate activated SFKs is the dominant role (R 3 , R 4a , R 22 and R T CRp ).Since CD45 can possess both positive and negative roles simultaneously, the species CD45 p and CD45 n are not mutually exclusive; however, the intersection between these two sets, a CD45-related state possessing both roles (CD45 p n ), is not explicitly modeled or tracked.CD45 tr , representing completely inactive CD45 that is translocated away from its substrates, is also not explicitly modeled.n , which is the intersection between the two modeled states outlined in purple, has both positive and negative roles, but is not explicitly modeled.CD45 tr represents completely inactive CD45 that is translocated away from its substrates and is also not explicitly modeled.In contrast to the dual role of CD45, Src homology region 2 domain-containing phosphatase-1 (SHP1) functions primarily as a negative regulator of TCR signal transduction.SHP1 associates with and negatively regulates the Syk family kinase Zap70 upon T-cell activation (R 8 and R 9 ), thereby suppressing TCR signaling [49].SHP1 also forms a negative feedback loop that is composed of (1) SHP1 phosphorylation by activated Lck, (2) binding of phospho-SHP1 to Lck and (3) Lck inactivation by SHP1-mediated dephosphorylation (R 3 ) [50].Erk, on the other hand, antagonizes SHP1 activity by modifying Lck (S59 phosphorylated by MAPK/Erk), which interferes with SHP1 recruitment and Lck inactivation [50,51].These events are captured by partitioning the SFK states with S59 phosphorylation (i.e., SF KdpS59p and SF KactS59p) from those without (R 20 and R 21 ).
Parameters for this updated model structure, which include reactions R 1 -R 10 and R 20 -R 23a , are trained to recapitulate the behaviors and roles of CD45 and SHP1 in early TCR signaling.To effectively isolate the early TCR signaling module from downstream feedback, the Erk positive feedback loop is suppressed with the presence of the Mek1/2 inhibitor U0126.Figure 4 shows the changes in the Zap70-Y319 phosphorylation level with different doses of αCD3 stimulation in the presence of U0126.The results show that with increasing doses of stimulation, the response of Zap70 increases slowly initially, then rapidly between 2 and 5 µg/mL of αCD3 and eventually saturates and even decreases slightly at the highest stimulus level (100 µg/mL).The dose response curve demonstrates the existence of a threshold αCD3 stimulation concentration.Below the threshold level, CD45 is unable to activate enough SFK and indirectly Zap70 to overcome the negative regulation by SHP1 and CD45 itself to sustain signaling.Once the threshold is crossed, the negative feedback barrier caused by SHP1 is overcome, and Zap70 signaling increases rapidly and eventually saturates at a high stimulus level.The Zap70 dose response simulations show good agreement with the observations originally reported by Zheng [40].Simulation Data (Zheng, 2005) αCD3 (µg/ml)

Tuning Erk Signaling
The Erk signaling cascade plays an important role in IL-2 activation via the transcription factor, AP1.Activated Zap70 mobilizes a number of adapter proteins to the receptor complex, including linker of activated T-cells (LAT) and growth factor receptor-bound protein 2 (Grb2) [52].Son of Sevenless (SOS), a Ras guanine exchange factor (GEF), forms a complex with Grb2 and LAT to facilitate the activation of Ras in T-cells.LAT also facilitates the activation of Ras through the RAS guanyl nucleotide-releasing protein-1 (RasGRP1), which is itself activated by DAG [19].Ras triggers the extracellular signal-regulated kinase (Erk) cascade that results in the activation of AP1 [20].As mentioned in Section 3.2, Erk also promotes its own activation by phosphorylating Lck at S59 and antagonizing SHP1 activity [50,51].
These reactions (R 11 -R 19 ) are modeled using the first-and second-order mass action kinetic equations of the base model presented by Zheng [40].The kinetic parameters in this module are fitted to Erk phosphorylation data presented by Perley et al. [41].Figure 5A,B shows the changes in the Erk phosphorylation with different doses of the MKP inhibitor sanguinarine and Mek1/2 inhibitor U0126.For each experiment, Jurkat cells were stimulated with 10 µg/mL αCD3, followed by a dose of inhibitor once the indicated amount of time had past.Samples were taken before and after the inhibitor dose administrations to show how the system changes over time with varying degrees of inhibition.At the lowest concentrations (where inhibition is negligible), the transient nature of Erk phosphorylation can be seen by the decreasing level over time.As the inhibitor concentrations increase, a threshold is crossed, resulting in more forceful inhibition: 10 µg/mL for sanguinarine and 1 µg/mL for U0126.The phospho-Erk simulations show good agreement with the time course and dose response observations originally reported by Perley et al. [41].
Figure 5. Erk activation in response to various doses of (A) MKP inhibitor sanguinarine and (B) Mek1/2 inhibitor U0126.Doses of sanguinarine and U0126 were administered at 15 and 6 min post-stimulation (10 µg/mL αCD3), respectively.Samples of phospho-Erk were taken at indicated times relative to the inhibitor doses and measured via western blot.The data shown are the means and standard errors from at least three independent experiments (data source: [41]).
The model is modified to recapitulate calcium signaling with these reactions (R 13 , R 27 -R 30 ), modeled using first-and second-order mass action kinetic expressions similar to that of the Erk signaling module.The kinetic parameters in this module are fitted to Ca 2+ flux data for which the experimental methodology is presented in Appendix 3. Figure 6 shows the experimental observations of intracellular calcium in Jurkat cells stimulated with 10 µg/mL αCD3 and the corresponding model simulation.Upon TCR stimulation, calcium ions are rapidly released from intracellular stores, peaking after 2 min before returning to a slightly elevated level.It is evident from the plot that the model readily captures the rapid transience of intracellular calcium flux.

Modeling CD28 Costimulation
In spite of the ability of the TCR to upregulate many factors necessary for full activation, TCR-antigen recognition ultimately leads to anergy, or a state of hyporesponsiveness, in the absence of a second signal [28].By contrast, signaling via CD28 on the surface of T-cells, in addition to TCR signaling, greatly enhances IL-2 production, T-cell proliferation and the prevention of anergy.As a result, CD28-mediated costimulation provides a critical second signal in T-cell activation and fate determination [27].CD28 signaling is believed to be mediated by phosphoinositide 3-kinase (PI 3 K) (R 33 ) [54].PI 3 K phosphorylates PIP 2 to become phosphatidylinositol 3,4,5-trisphosphate (PIP 3 ), an activity that is directly antagonized by phosphatase and tensin homolog (PTEN) (R 34 ).PIP 3 serves as pleckstrin homology (PH) domain membrane anchors for 3-phosphoinositide-dependent kinase-1 (PDK1) (R 35 ).PDK1 activation leads to the membrane recruitment of PKCθ (R 36 ) and protein kinase B (PKB, or Akt) (R 40 ) to regulate several signaling pathways, including NFκB [23,29] and mTOR [30].
The model is augmented to simulate the effects of CD28 costimulation on T-cell activation.As with the other signaling modules, these reactions (R 32 -R 35 ) are modeled primarily using first-and second-order mass action kinetic expressions.The corresponding kinetic parameters are tuned after accounting for downstream signaling events, such as NFκB and mTOR signaling, which are described in the following sections.

Modeling NFκB Signaling
Several lines of evidence indicate that the NFκB pathway is perhaps the most relevant biochemical or transcriptional target for the costimulatory activity of CD28 [29].TCR-and CD28-mediated induction of the NFκB signaling pathway intersect at the central regulator, PKCθ (via mechanisms described in Section 3.5).PKCθ induces NFκB signaling by activating IκB kinase (IKK) (R 37 ) [21], which phosphorylates the inhibitor IκB (at positions S32 and S36).This triggers the rapid polyubiquitination (at positions K21 and K22) and proteolysis of IκB in the 26S proteasome complex (R 38 ).IκB degradation exposes the nuclear localization signal of NFκB, allowing its rapid translocation into the nucleus (R 39 ) [23,24].
Physiologically, the defining characteristic of IκBα is its ability to regulate rapid, but transient, induction of NFκB activity, owing to the participation of IκBα in an autoregulatory feedback loop.That is, the activation of NFκB causes the upregulation of transcription of IκBα, which, in turn, serves to shut off its own nuclear localization signal [36,37].This upregulation occurs due to the presence of κB sites in the IκBα promoter.Thus, IκBα is thought to maintain the transient effect of inducing agents on the transcription of NFκB responsive genes [23].It was demonstrated that in vivo degradation of IκBα is required for the appearance of NFκB in the nucleus.In addition, some investigators were able to demonstrate by co-immunoprecipitating rel proteins with IκBα, from stimulated cells treated with proteasome inhibitors to block the degradation of phosphorylated IκBα, that IκBα undergoes degradation mediated by the 26S proteasome after phosphorylation, but before dissociation [23].
We model the autoregulatory relationship between IκBα and NFκB (R 38 and R 39 ) as a two-component negative feedback system: The dynamic behavior of the negative feedback depends on the relative efficiency of the feedback regulation (α and β) regulating oscillation persistence versus self-regulation (γ and δ), causing oscillation damping.The output, N F κB, can range from persistent oscillations (high feedback efficiency and no damping, α = δ = 0) to gradual rising to a plateau level (low feedback efficiency and high damping).Figure 7 demonstrates the phosphorylation of IκBα by IKK, subsequent rapid degradation by the 26S proteasome and synthesis of IκBα promoted by nuclear NFκB.The corresponding NFκB signal is also depicted.Model parameters are fitted to phospho-IκBα data for which the experimental methodology is specified in Appendix 4. In our case, the fitted parameters correspond to damped oscillations (intermediate feedback efficiency and intermediate damping).Although we note that the sample size is small, the model simulation does show good agreement with the trend seen in the data.Mammalian target of rapamycin (mTOR) is an evolutionary conserved serine/threonine protein kinase that is well known for its ability to control T-cell activation and differentiation.mTOR is necessary for TCR-induced signaling to drive differentiation into Th1, Th17 or Th2 effector types in vitro or in vivo under polarizing conditions.On the other hand, mTOR-deficiency drives CD4+ T-cells to the generation of FOXP3+ T-cells, even under normal activating conditions [30].The ability of mTOR to regulate the expression of FOXP3 makes it a key player in the development of Tregs [30,31].
The model is augmented with the described reactions (R 40 -R 47 ), which consist primarily of first-and second-order mass action kinetic expressions.Hill equations are employed in cases, such as PTEN deactivation by TCR triggering (R 45 ), TSC2 phosphorylation (R 41 ) and IL-2 (R 46 ) and FOXP3 (R 47 ) transcription, to better describe apparent cooperativity between transcriptions factors [26,31].Kinetic parameters are chosen to reproduce qualitative observations of Akt and PKC phosphorylation in response to αCD3 and αCD28 costimulation originally reported by Lee et al. [32], the results of which are shown in Figure 8.In the model, as with the data, while αCD3 and αCD28 are each capable of inducing Akt and PKCθ phosphorylation, both are required to induce full activation of these pathways.
Figure 8. T-cell signaling in response to doses of αCD3 and αCD28 as measured by (A) P-Akt and (B) P-PKCθ.CD4+ T-cells were stimulated (40 min) with 0.5 mg/mL plate-bound αCD3, 2.5 mg/mL of soluble αCD28 or both.Bar graphs quantify phosphorylation of Akt and PKC, with each sample normalized to the level of unphosphorylated protein in one experiment representative of three replicates (data source: [32]).

Results and Discussion
The following results and discussion subsections are divided into two main groups.First, the complete and calibrated model, as described in Section 3, is evaluated against additional data not used in the model tuning process.This experimental evidence and corroboration of the model's predictive accuracy is described in Section 4.1.The sections that follow then describe predictions made using the model.As mentioned in the introduction, these will mainly focus on the roles of signal strength and feedback loops on the regulation of T-cell activation.

Experimental Datasets for Model Corroboration
Experimental datasets for the purposes of model corroboration are compiled from the published literature.Only those that are commensurate with our model in terms of cell type, stimuli, output species and time scale are considered.Table 2 summarizes the attributes of the compiled dataset.Each dataset is used to analyze the model's ability to recapitulate the dynamics of the biological system as indicated.These datasets consist of dynamics (i.e., time courses), dose responses (i.e., multiple doses of an input or stimulus), input responses (i.e., single doses of different input or stimuli combinations) and knockdown or knockout (i.e., activity of a particular species is reduced or eliminated) experiments.In order to compare the model and observed behaviors directly, simulations are generated to mimic the experimental conditions used in the published experimental set-ups.Model simulations are shown with the corresponding data and error bars whenever possible.In order to corroborate the model, we perform several in silico experiments and compare the simulations to published experimental results (described in Table 2).Figure 9A shows the changes in the intracellular calcium release in response to different combinations of stimuli.As shown in the figure, only αCD3 and ionomycin are capable of inducing calcium signaling, with ionomycin (administered with PMA) being the stronger of the two inducers.For the most part, the model simulations are in good qualitative agreement with the observations reported by Smeets et al. [27]; however, there are two noticeable discrepancies.First, these data show that the model slightly overestimates the system's calcium sensitivity to ionomycin stimulation relative to that of αCD3.Since the model is shown to fit calcium signaling quite well (see Figure 6), more information on the effects of ionomycin would likely resolve this issue.Second, the model's predicted calcium response demonstrates a slight sensitivity to PMA, a trend that is not observed in the data.This model behavior is likely caused by the presence of the positive feedback loop mediated by Erk.PMA is considered a DAG substitute, activating the substrates of DAG, such as RasGRP, which leads to the activation of Erk.Erk, in turn, promotes signaling by preventing SFK from activating the SHP1 negative feedback loop.In the model, the effect of elevated SFK facilitates further signaling through LAT, PLCγ, IP 3 and eventually calcium.This may be evidence that the model, particularly the positive feedback loop, may be more sensitive than these experimental data suggest.
Figure 9B shows changes in nuclear translocation of NFAT, a downstream target of calcium signaling, with varying doses of the calcineurin inhibitor cyclosporin A (CsA).Model simulations show excellent agreement with the observations reported by Clipstone et al. [53] in terms of half-maximal inhibitory concentration (IC 50 ) at ~3.5 ng/mL CsA.However, the model does appear to be slightly less sensitive to changes in inhibitor concentrations, as the slopes of the two curves are slightly different.
Figure 10 shows the activity of transcription factors NFAT, NFκB and AP1, as well as IL-2 synthesis in response to different combinations of stimuli.While there are a few obvious quantitative discrepancies between the model simulations and the data presented by Smeets et al., there is good qualitative agreement among them in the input/output relationships.Nuclear localization of NFAT is only achieved when stimulating with combinations involving αCD3 (Figure 10A).This is because NFAT translocation requires strong calcium signaling, which αCD28 and PMA stimulation alone could not induce, consistent with [17] and Figure 9A.On the other hand, both αCD3 and PMA are sufficient to induce nuclear localization of NFκB and AP1 (Figure 10B,C).This is promoted by αCD28, but a combination of αCD3 and PMA have the greatest influence on translocation.IL-2 transcription in response to these same stimuli is shown in Figure 10D.The results indicate that TCR signaling is necessary for IL-2 transcription, but not sufficient.Coupling with CD28 coreceptor signaling causes a moderate increase in IL-2 transcription, but coupling with PMA results in the greatest observed increase, consistent with [17,27].These experimental results corroborate the model and suggest that it is capable of accurately recapitulating the transcription profiles and, indirectly, the upstream signaling pathways, for a wide variety of stimuli.

Corroboration of CD45 Activity
In this model, CD45 exists initially as fully active (proximal to the SFKs), then is deactivated (translocates away from the SFKs) by the formation of the TCR complex.At later time points, CD45 is reintroduced to the receptor cluster with an increased negative regulatory role.This results in dephosphorylating active SFK and phosphorylated ligand-bound TCR, thus terminating the TCR signal.Figure 11A depicts Lck-Y505 phosphorylation as a function of CD45 activity after 15 min of stimulation.At normal expression levels (i.e., 100% of WT), Y505 phosphorylation is reduced by approximately 75%.Phosphorylation at Y505 is inversely related to CD45 activity, which is corroborated by the observations of McNeill et al. [55].The model is able to approximate the data from 100% WT activity down through 5% WT activity very well.Only when CD45 is completely suppressed do the model and data diverge, although the trend is still present.[55]).(B) Model simulations of Erk phosphorylation 3 min after stimulation with 0 or 50 µg/mL αCD3 as a function of CD45 activity (data sampled from [55]).(C) Model simulations of intracellular calcium release over time in wild-type (WT) and CD45 knockdown mutant (5% of WT activity) stimulated with 10 µg/mL αCD3 at 90 s as a function of CD45 activity (data sampled from [56]).In CD45-null thymocytes, the p56Lck and p59Fyn tyrosine kinases are hyperphosphorylated, and p56Lck is found in its inactive conformation [56].Both basal and TCR-stimulated tyrosine phosphorylation of TCRζ and CD3 are also much reduced.These defects are associated with the failure of Zap70 kinase recruitment to the TCRζ chain; however, TCR-induced signaling is not entirely ablated.Figure 11B shows Erk phosphorylation as a function of αCD3 stimulation and CD45 activity.It is clearly evident that pErk increases with CD45 activity; however, above ~50%, there is qualitative switch, and pErk begins to decrease as CD45 expression approaches wild-type levels, which is also corroborated by observations by McNeill et al. [55].
Furthermore, significant inositol phosphate and calcium signals are observed in CD45-null thymocytes.In our simulation, the CD45 defect is approximated with a knockdown to 5% of WT activity as complete CD45 knockout suppressed all signaling.Although greatly reduced from the nominal system, calcium signaling does not appear to be insignificant in CD45-defective CD4+ T-cells (Figure 11C).The molecular analysis presented by Stone et al. suggests that the threshold for TCR signal transduction is greatly increased in CD45-null T-cells, thus explaining the profound defects in thymic development [56].

Corroboration of SHP1 Activity
SHP1 functions as a negative regulator by deactivating Zap70 and SFK upon T-cell activation [49].SHP1 also forms a negative feedback loop that is composed of SHP1 phosphorylation by activated Lck, binding of phospho-SHP1 to Lck and Lck inactivation by SHP1-mediated dephosphorylation [50].To evaluate the ability of the model to capture the role of SHP1 feedback on T-cell activation, we simulate the SHP1 knockdown experiment presented by [49]. Figure 12 shows the stimulation of the IL-2 reporter with different combinations of stimuli and two different levels of SHP1 activity.Model-simulated results and data are depicted by bars and dots, respectively.In the wild-type system, neither PMA nor αCD3 alone are sufficient to upregulate IL-2 transcription.However, the combination of the two is able to generate a substantial signal.Pairing PMA with ionomycin produces the maximally-observed signal.In the SHP1-knockdown system, IL-2 synthesis is drastically increased as a result of αCD3 and PMA + αCD3 stimulation, thus demonstrating the lack of inhibition.PMA + ionomycin stimulation is not affected by SHP1 knockdown because SHP1 functions upstream of their substrates: RasGRP1, PKCθ and calcium.The experimental results show mostly good agreement with the model, suggesting that the model is able to adequately capture the reported behavior of SHP1.

Weak CD28 Costimulation Predicted to Elevate FOXP3 Transcription
TCR signal intensity and costimulation is known to differentially affect the activation of T-cells.We study this by performing a two-way dose response experiment, measuring the effects of αCD3 and αCD28 on IL-2 and FOXP3 transcription.Figure 13A,B shows IL-2 and FOXP3 transcription, respectively, following 30 min of stimulation by various combinations of αCD3 and αCD28.At trivial doses of both stimuli, neither species is active, indicating that the cell is at rest.At high doses of both stimuli, IL-2 is the dominant species, indicating that the cell is fully active.However, when costimulation through CD28 is relatively weak, the balance of dominance shifts toward FOXP3, while leaving IL-2 at a much reduced level.Indeed, this behavior is corroborated by the observations by Kretschmer et al. [57] that weak TCR signals and limited costimulation have been linked to FOXP3 induction and the development of the regulatory phenotype in CD4+ T-cells.
Figure 12.IL-2 reporter stimulation in response to SHP1 knockdown.The model simulates the wild-type and SHP1 knockdown mutant (C453S mutation, resulting in catalytically inactive SHP1) stimulated with 10 µg/mL αCD3.IL-2 reporter stimulation was measured 2 h post-stimulation (data sampled from [49]).We investigate the roles of CD45 and SHP1 on T-cell activation by simulating a two-way knockdown response experiment and measuring the effects on a small number of important outputs.Figure 14A-F show the model response to various CD45-and SHP1-knockdown scenarios.SHP1 clearly has the effect of downregulating most species playing a role in T-cell activation.As SHP1 activity is reduced (i.e., moving top to bottom on the right axis in each plot), calcium, Erk, PKCθ and IL-2 are all upregulated.This is due to the loss of their transient behaviors and constitutive activation at elevated levels.For CD45, on the other hand, the analysis demonstrates a dual role as both activator and inhibitor of T-cell activation.As CD45 activity is reduced (i.e., moving right to left on the top axis in each plot), TCR signaling increases until CD45 activity reaches ~1%-10% of the wild-type, then decreases to complete inactivation as CD45 is completely suppressed.This result indicates that CD45 is required for TCR signaling, without which the pool of SFKs would remain in the closed and inactive state; however, the highest level of activation is not actually achieved for wild-type CD45 expression, but rather at a reduced rate between approximately 1%-10% of the wild-type with intact SHP1 activity.Stronger TCR signaling can be achieved at full CD45 activity by downregulating SHP1.It is evident that the system is quite sensitive to CD45 activity, particularly at the lower end of expression.By contrast, SHP1 is much more effective closer to wild-type levels; however, its activity is heavily dependent on the presence or absence of CD45.The results suggest that there exists a distinct threshold of CD45 activity (~1% of WT) that may be critical to determining the outcome of TCR-mediated activation.Below this threshold, the system is presumably unable to activate with every state remaining at baseline.However, at the threshold and above, the signaling molecules relating to Th development become active and the system is rescued.Recently, McNeill et al. reported that only 3% of normal CD45 activity is sufficient to reconstitute CD45-deficient mice with normal numbers of mature T-cells, biased toward CD4+ T-cell lineage commitment [55].
CD45 also appears to have an effect on the activation of FOXP3.At total CD45 suppression, the cell remains unresponsive, even at full stimulation.For wild-type activity, TCR signaling promotes mTORC1 activation, thereby also suppressing FOXP3.However, a reduction in CD45 activity to around the aforementioned threshold actually induces a substantial upregulation in FOXP3 expression while IL-2 remains low.This may be due to the transcriptional regulators of FOXP3 becoming sufficiently induced, while signaling remains too low to promote mTOR-mediated inhibition.This result suggests that CD45 downregulation may lead to weak IL-2 induction and increased FOXP3 expression, which is often observed in Treg lineage commitment.

Conclusions
Signal intensity and feedback regulation are known to be major determinants of the signaling profile for the TCR and its various coreceptors.While the exact nature of these relationships remains under investigation, it is believed to involve a complex signaling network with cross-talk between the calcium, Erk, PKCθ and mTOR signaling pathways.
In this manuscript, we present a mathematical model encompassing the signal transduction events relating to the process of TCR-mediated cell activation and gene expression.The model is able to reproduce key behaviors in: (1) ligand-induced TCR trafficking, synthesis and degradation; (2) early kinase and phosphatase interactions between SFK, Zap70, CD45 and SHP1; (3) CD28 costimulation; and (4) downstream signal transduction pathways leading to IL-2 and FOXP3 synthesis, including calcium, Erk, PKCθ and mTOR.
In addition to corroborating many experimentally-observed behaviors, the model is able to provide insight into the positive and negative regulatory roles of the phosphatases CD45 and SHP1 during T-cell activation.Analysis of the model demonstrates that: (1) SHP1 negative feedback is necessary for preventing hyperactivity in TCR signaling; (2) CD45 is required for TCR signaling, but also partially suppresses it at high activity levels; and (3) elevated FOXP3 and reduced IL-2 signaling, an expression profile often observed in developing Tregs, can be achieved either by weak TCR and CD28 stimulation or a severe reduction in CD45 activity.However, we do note that further investigation and experimental evidence is required in order to corroborate these predictions.
While the proposed model is a demonstrably powerful tool for predicting many events involved in CD4+ T-cell activation, there are certain caveats that need to be considered to accurately define the model's scope and capability.First, as a finite mathematical approximation, this model is inherently an abstraction of the biological reality.As such, the model does not attempt to explain every possible mechanism involved in the processes in question.For example, the Foxo family protein, Foxo1, is thought to bind to the FOXP3 locus and induce FOXP3 gene transcription.Foxo1 activity is subject to modulation by Akt kinase signaling, and Tregs have dampened Akt signaling in response to TCR stimulation compared with conventional T-cells.However, it is not within the scope of this work to test and corroborate every one of these mechanisms.Second, the purpose of the model is to simulate the dynamics of a small number of key species involved in T-cell activation or, more specifically, and clinically relevant, activation of human primary CD4+ T lymphocytes.Due to practical limitations on the availability of experimental data, however, the model is partially calibrated with datasets from alternate sources, including Jurkat cells.While Jurkat cells are an immortalized cell line of human T lymphocytes, they share many similarities with their primary counterparts and are very useful in studying T-cell signaling and IL-2 production.However, we do recognize that this constitutes an amalgamation of data from a variety of sources and that further investigation using primary T-cells is highly desirable.
As demonstrated in this manuscript, the mathematical model captures the key events in TCR and CD28 co-mediated signal transduction events leading to IL-2 and FOXP3 activation.As such, the model enables researchers to study these processes with a combination of broad scope, using a large-scale highly-connected network and with quantitative detail and accuracy not allowed by comparable models.In the future, we plan to use this mathematical model to design informative and hypothesis-driven experiments to refine our understanding of the dynamical nature of CD4+ T-cell activation.Our goal is to use such a model as a foundation for the design of strategies to drive and control T-cell activation and differentiation for therapeutic gain.Table A2 presents a comprehensive list of all reaction equations used to model the biochemical reactions of this system (illustrated in Figure 1).These expressions form the basis for the rate equations of the ODE-based model presented in Table A1.The list includes the reaction identifier, mathematical equation and biological meaning for each reaction.Table A3 presents a comprehensive list of all reaction parameters included in the model.The list also includes the biological meaning, value, 95% confidence interval (if applicable), units and source (if applicable) for each parameter.For parameters not provided by external sources, they are either estimated from data or explicitly derived to satisfy a condition of the model, for example, to ensure equilibrium when the system should be at rest.Confidence intervals are estimated for the model calibration data (Table 1) using the likelihood ratio method [59].

Figure 1 .
Figure 1.Reaction network of the T-cell activation model.Solid arrows denote reactions for which the forward and reverse directions are indicated; dashed arrows connect reactions (either forward or reverse) with their catalysts.Colored arrows denote reactions that are catalyzed by specific species as indicated.Symbols and reactions are described in Appendix 1.

Figure 2 .
Figure 2. Resting and ligand-induced (A) surface TCR expression and internalization and (B) degradation activity in response to varying doses of αCD3 stimulation (data shown in [45]).

Figure 3 .
Figure 3. Representation of the phosphatase CD45 in the model.The model states CD45 p and CD45 n are represented by the regions outlined in blue and red, respectively.CD45 p n , which is the intersection between the two modeled states outlined in purple, has both positive and negative roles, but is not explicitly modeled.CD45 tr represents completely inactive CD45 that is translocated away from its substrates and is also not explicitly modeled.

Figure 4 .
Figure 4. Zap70-Y319 phosphorylation in response to various doses of αCD3 stimulation in the absence of Erk feedback.Jurkat cells were incubated in the presence of a Mek1/2 inhibitor (2 µg/mL U0126) and stimulated with αCD3 at the indicated concentrations.Samples were taken 5 min post-stimulation and analyzed by western blot.The data shown are the means and standard errors from at least three independent experiments (data source:[40]).

Figure 6 .
Figure 6.Intracellular calcium release in response to stimulation with 10 µg/mL αCD3.The data shown are the means and standard errors from 12 independent experiments.

Figure 7 .
Figure 7. IκBα phosphorylation following 2 µg/mL αCD3 and αCD28 costimulation.The data shown are quantified western blot data from one experiment.

Figure 11 .
Figure 11.Downstream TCR signaling in response to CD45 knockdown.(A) Model simulations of Lck phosphorylation at the negative regulatory residue Y505 as a function of CD45 activity (data sampled from[55]).(B) Model simulations of Erk phosphorylation 3 min after stimulation with 0 or 50 µg/mL αCD3 as a function of CD45 activity (data sampled from[55]).(C) Model simulations of intracellular calcium release over time in wild-type (WT) and CD45 knockdown mutant (5% of WT activity) stimulated with 10 µg/mL αCD3 at 90 s as a function of CD45 activity (data sampled from[56]).

Figure 13 .
Figure 13.(A) IL-2 and (B) FOXP3 transcription in response to stimulation by combinations of αCD3 and αCD28.Results show outputs 30 min after stimulation.

Table 1 .
Description of experimental datasets used for model development.

Table 2 .
Description of experimental datasets used for model corroboration.

Table A1 .
Summary of model states.

Table A2 .
Summary of model equations.

Table A3 .
Summary of model parameters.