A Mathematical Model of T Lymphocyte Calcium Dynamics Derived from Single Transmembrane Protein Properties

Fate decision processes of T lymphocytes are crucial for health and disease. Whether a T lymphocyte is activated, divides, gets anergic, or initiates apoptosis depends on extracellular triggers and intracellular signaling. Free cytosolic calcium dynamics plays an important role in this context. The relative contributions of store-derived calcium entry and calcium entry from extracellular space to T lymphocyte activation are still a matter of debate. Here we develop a quantitative mathematical model of T lymphocyte calcium dynamics in order to establish a tool which allows to disentangle cause-effect relationships between ion fluxes and observed calcium time courses. The model is based on single transmembrane protein characteristics which have been determined in independent experiments. This reduces the number of unknown parameters in the model to a minimum and ensures the predictive power of the model. Simulation results are subsequently used for an analysis of whole cell calcium dynamics measured under various experimental conditions. The model accounts for a variety of these conditions, which supports the suitability of the modeling approach. The simulation results suggest a model in which calcium dynamics dominantly relies on the opening of channels in calcium stores while calcium entry through calcium-release activated channels (CRAC) is more associated with the maintenance of the T lymphocyte calcium levels and prevents the cell from calcium depletion. Our findings indicate that CRAC guarantees a long-term stable calcium level which is required for cell survival and sustained calcium enhancement.


INTRODUCTION
The dynamics of the free cytosolic calcium concentration upon stimulation of T lymphocytes (TCs) is crucial for TC activation and fate decision processes. While it is clear that calcium rises upon stimulation of the TC receptor (TCR), the calcium pattern associated with different fates of TCs has not been deciphered (1)(2)(3)(4). However, it is likely that the calcium signal is correlated with the later fate of the activated TC, i.e., anergy, division, acquisition of the regulatory phenotype or apoptosis (3,5,6). Dysregulation in TC calcium signaling has been linked to inflammatory and autoimmune diseases as well as to allograft rejection (2). A relevant player in calcium dynamics is the calcium-release activated channel (CRAC) which is located in the plasma-membrane and activated by calcium depletion in the intracellular calcium stores like the endoplasmic reticulum (ER). As ion-gating transmembrane proteins in the plasma-membrane (PM) are possible targets of drug applications in the context of various clinical settings (2,7,8), insight into the specific calcium dynamics is essential for an efficient control of TC behavior.
The relative contribution of ER-derived calcium entry versus CRAC to the calcium signal following TC stimulation is a matter of ongoing debate (9)(10)(11)(12). While a considerable number of scientists argue for CRAC being the major player of TC calcium dynamics (12,13), others argue for a dominant role of calciuminduced calcium-release (CICR) (10,14,15) or for a dominant role of second messenger-induced calcium-release from ER (16,17). All three contributions are required for a functional TC calcium signal, however, the sequence of the contributions might be essential. Second messenger-induced activity appears as a very early signal (18), which might act as triggering event for CICR and subsequent CRAC activation (2,(19)(20)(21). Quantitative analysis of the components of calcium signaling during TC activation is essential for the development of strategies for an efficient control of TC responses. A mathematical analysis of the calcium dynamics in a model including calcium stores and CRAC may shed light on the relation and relevance of both calcium sources (12,22). This is the major motivation for the present work.
T lymphocytes are non-excitable cells in the sense that they do not exhibit bursts like pancreatic beta-cells, spikes like neurons, or comparable fast dynamics (23)(24)(25) even though single cell measurements detected calcium oscillations (10,19,26). The non-excitability of TCs might have prevented a larger interest of mathematical modelers in lymphocyte calcium dynamics. The few existing models (22,(27)(28)(29) mostly focused on modeling of CRACchannel dynamics or a special part of the store-operated-calciumentry signaling pathway (27) and its contribution to intracellular calcium dynamics. Also the dependence of ORAI1 assembly to a tetrameric CRAC on calcium oscillations was considered (28). In one approach a spatial resolution of CRAC currents and of the calcium dynamics in TCs was considered in the context of immunological synapse formation (22). Two different models for inositol 1,4,5-tris-phosphate-receptor (IP3R) activity were compared and it was shown that they differ in their impact on TC calcium dynamics (27), a result that will be used in the present model as well. The plasma-membrane calcium-ATPase (PMCA) was modeled in Jurkat TCs and a reversible modulation of PMCA activity was postulated (12), which is a further topic addressed in this investigation.
The aforementioned theoretical studies on TC calcium dynamics were all based on whole cell current models. In the context of excitable cells we have shown that it is possible to derive the whole cell currents from the single transmembrane properties (30). To achieve this goal, specific quantitative measurements of protein activation, inactivation, dependencies and conductance were used. The big advantage of this approach is that most parameters of the model are determined by independent experiments, which increases the predictive power of the model. The present paper applies this strategy to calcium dynamics of TCs. The model also includes the dynamics of transmembrane channel expression kinetics in order to represent CRAC recruitment upon calcium store depletion. The mathematical model was validated using dynamic calcium data measured under specific experimental conditions. The validated model allowed to reassess the relative role of store-derived and CRAC-mediated calcium entry on a quantitative basis. A new role of the CRAC-channel is postulated, which is associated with maintenance of TC calcium levels rather than TC activation.

MATERIALS AND METHODS
The modeling framework is presented in this section. Three compartments, extracellular space (ES), cytosol, and ER are distinguished, each being represented by ordinary differential equations. The nucleus is only included as an object which reduces intracellular space (Section 2.4). The compartments are encased by the PM and the membrane of the ER. Both membranes contain transmembrane proteins (Figure 1), which allow for a flow of ions from one compartment to the other. The surface densities and the properties of these proteins in terms of conductance and control parameters determine the resting and activated states of the TC. While the protein properties are derived from measured single protein data, which are independent of the whole cell calcium experiments used for validation, the protein densities are in parts derived from steady state conditions (Section 2.9) and in parts determined by fitting to whole cell calcium dynamics (Section 2.10).
In the following, the equations for calcium dynamics, calciumbuffering, and for the second messenger d-myo-inositol-1,4,5tris-phosphate (IP3) are introduced. The details of the compartment sizes and the surface between the compartments are explained in Section 2.4. Particular emphasis is put on the geometry of the Jurkat TC, the specific cell line, which is used in the subsequently described experiments. The single transmembrane protein characteristics will be described and the corresponding FIGURE 1 | Scheme of the transmembrane proteins included in the mathematical model. The elements of a TC considered in the mathematical model are shown with particular emphasis on calcium-conducting transmembrane proteins. The outer border of the cell represents the PM. In the PM CRAC and PMCA are located which induce calcium influx and extrusion, respectively. The intracellular organelle around the nucleus (yellow) depicts the ER. The membrane of the ER contains IP3R and SERCA (sarco/endoplasmic reticulum calcium-ATPase) which control the exchange of calcium between the ER and the cytosol. mathematical models introduced. Wherever possible we implemented data from experiments performed with Jurkat TCs. Finally, all model pieces are merged to the proposed model of TC calcium dynamics in Section 2.10. This includes the determination of the remaining unknown parameters.

CALCIUM DYNAMICS AND BUFFERING
Two major players determine exchange of calcium through the PM (Figure 1): PMCAs actively transport calcium from the cytosol into ES (12,31), while CRAC-channels allow for a passive electrochemical current of calcium ions into the cell (21,32,33). The density of active CRAC-channels ρ CRAC in the PM is increased in dependence on ER calcium (C ER ) depletion (34). The free cytosolic calcium concentration (C) is further affected by modulations of the calcium flow between cytosol and ER (Figure 1). Sarco/endoplasmic reticulum calcium-ATPase (SERCA) transports calcium ions from the cytosol into the ER (10,35) and by this maintains a chemical gradient of calcium from the ER to the cytosol. Conversely, calcium can passively leave the ER into the cytosol when IP3R channels open in dependence on calcium and the second messenger IP3 (1, 36, 37) associated with calciuminduced-calcium-release (CICR) (38). Furthermore other second messengers like cyclic ADP ribose (cADPR) (16) and nicotinic acid adenine dinucleotide phosphate (NAADP) (39) were found to influence calcium dynamics. Their effect is mediated by activation of the ryanodine receptor (RyR) which leads to calcium-release from intracellular stores (17,(40)(41)(42)(43). In order to avoid an overparametrization of whole cell calcium curves, which the present model focusses on, we restrict ourselves to the dynamics of IP3. The inclusion of cADPR and NAADP and their effect on RyR Frontiers in Immunology | T Cell Biology requires a more detailed data basis and should be addressed with a more complex model in the future.

Calcium in the cytosol
The four sources and sinks of free cytosolic calcium C are described as where z Ca = 2 is the valence of calcium ions, F the Faraday constant, ξ and ξ ERC geometrical surface to volume ratios for PM [equation (16)] and ER membrane [equation (17)], respectively. ρ X is the surface density and I X the single transmembrane protein current, which are defined in the subsequent sections. By convention, positive ions that enter the cytosol are represented by negative electrical currents. B C represents the cytosolic calcium-buffer in the rapid buffer approximation where b 0 is the total buffer concentration, and K b the calciumbuffer dissociation constant. The fraction of free calcium in the cytosol then reads ( The main buffer within the cytosol is calmodulin (CaM) with 4 calcium binding sites per CaM. There is a diversity of measured CaM concentrations depending on cell type and organ (44)(45)(46). A realistic average value is 25 µM of CaM, corresponding to b 0 = 100 µM of calcium binding sites. The dissociation constant K b was determined by the required fraction of free calcium of f C ≈ 0.1% in non-excitable cells (47), leading to K b = 0.1 µM.

Calcium in the ER
The dynamics of the calcium concentration in the ER C ER is described by an equation analogous to equation (1) with a different geometrical factor ξ ER [equation (18)], and a different calcium-buffer B C,ER defined by The fraction of free calcium in the ER reads The resting ER calcium level is C ER,0 = 400 µM (2) which holds true for Jurkat TCs considered here (34).
In the ER calcium is buffered by calsequestrin, with three high and three low-affinity binding sites (48), and by calreticulin, with two distinct domains, one with high-affinity (K = 0.01 mM) but low capacity (0.6-1 mol Ca 2+ /mol protein), and one with lowaffinity (K = 2 mM) but high capacity (18 mol Ca 2+ /mol protein) (49). As calreticulin binds more than 50% of the luminal ER calcium (50) only this buffer is considered here. In pancreatic acinar cells it was estimated that 20-times more calcium would be free in the ER compared to the cytosol (47), suggesting f C,ER ≈ 0.02. This is achieved by using b ER,0 = 30 mM and K ER,b = 0.1 mM, which corresponds to an intermediate dissociation constant of both calcium binding domains.

IP3 DYNAMICS
IP3 (P) is generated in a TCR-and calcium-dependent way described by where β P is the production and γ P the degradation rate. The Hill-function is defined as where K is the concentration of X at which the Hill-function reaches its half value, and n the Hill-coefficient which determines the steepness of the Hill-function. The degradation rate in equation (7) is determined by steady state conditions for IP3 in equation (35). The production rate is the tonic production rate and is modulated by increased calcium with the Hill-function in equation (7), leading to a positive feedback loop between calcium and IP3. β P is fitted as described in Section 2.10 and mainly influences the speed of the early calcium response after TC stimulation.
The production is further modified by the time-dependent input function T (t ) representing the degree of TCR stimulation of the cell. T = 1 is assumed in the resting state.
The resting concentration of IP3 P 0 is identified as critical parameter. It strongly determines the responsivity of the cell via activation of IP3R (see Section 2.6). It was fitted as described in Section 2.10.

MEMBRANE AND REVERSAL POTENTIALS
The resting membrane potential is set to V = −60 mV (51)(52)(53). It is assumed that the membrane potential is not changed by the calcium currents (V = V 0 ) and that the electrical current corresponding to calcium fluxes in or out of the cell is equilibrated by other ions.
Further it is assumed that V ER,0 = V 0 = −60 mV, thus, the ER and the cytosol are electrically equilibrated (54). ER calcium efflux may lead to small fluctuations (55) which are neglected. Thus, V ER = V is assumed at all times.
In this approximation, the reversal potentials depend on the chemical gradient only. The Nernst-equation is used to calculate the reversal potential during dynamical changes of calcium www.frontiersin.org concentrations: where R = 8.315 J/(K mol) the Rydberg (molar) gas constant, T = 310 K, and F the Faraday constant. For many channels, the real I-V-relationship is not linear as assumed for the currents I X in equations (23) and (28). Therefore, the reversal potential is corrected for the CRAC-channel by a shift V C = 78 mV in order to achieve the correct linear extrapolation of the I-V-relation of CRAC-channels with a zero around V C ≈ 50 mv [(34) Figure 1, (33) Figure 2]. This approximation is only valid for depolarization below V = 50 mV.
The reversal potential for ER calcium V C,ER is treated in complete analogy to the cytosolic case which leads to a correction of V C,ER . As the value is not known it is derived using the fitting routine in Section 2.10.

TC GEOMETRY
Most measurements on TC calcium dynamics are performed in Jurkat TCs which are small but still larger than normal human blood derived TCs. In an approach based on ordinary differential equations the effect of an ion current onto the concentration of the ion in the cytosol or ER is not spatially resolved. While local calcium entry can induce transient high concentrations of calcium (22) the comparably small cytosolic volume of TCs justifies this approach for the description of whole cell calcium dynamics because local inhomogeneities will quickly equilibrate. In the model this is reflected as change of the average concentration. How an ion current changes the average calcium concentration depends on the geometry of the cell. In the dynamic equations for the ion concentrations the current I X through an individual ionconducting protein X is multiplied by the surface density ρ X . Thus the concentration change is derived from a current surface density. The latter has to be translated to the actual change in concentration by a surface to volume ratio, which is considered here.
Given a cell radius R cell , the cell volume V cell and cell surface A cell are known as well. However, the volume relevant to changes in concentration is not the cell volume V cell but the cytosolic volume V cyt which can be approximated as using the volumes of ER and nucleus. This is important because the nucleus, with a radius of is substantially reducing the resulting cytosolic volume. f R ≈ 0.8 is assumed for human TCs (56), and f R ≈ 0.25 for Jurkat TCs. The volume of the ER is expressed as a fraction of the total cell volumẽ with f V ≈ 0.1 (57). However, electron micrographs of TCs suggest that f V ≈ 0.01 (25) which is used here. Taking this together, the cytosolic volume becomes The surface of the ER is also needed in order to translate the current surface densities calculated on the ER surface into concentration changes in cytosol and ER. While the TC itself is approximated as a sphere, the ER is absolutely non-spherical. The exact surface of the ER is difficult to be measured and accordingly approximated as whereÃ ER is the surface of a spherical ER with volumeṼ ER [determined in equation (12)], and f A is the fold increase of the ER surface with respect to the surface of a spherical ER. f A = 30 was roughly estimated from the folding degree of the ER in electron micrographs of TCs (25). Note that only the product of f A with f 2/3 V enters the model, such that both parameters are redundant and were only separated because of their physiological meaning.
The size of human blood TCs can be estimated starting from the capacity of C m = 0.028 pF/µm 2 (58,59) and using the whole cell capacitance of C cell = 2 pF [(60), p. 606]. C cell = 1.7 pF was found in Fomina et al. (14). Using C cell /A cell = C m this yields a radius R cell = C cell 4πC m (15) and the resulting R cell ≈ 2.4 µm corresponds to A cell = 72.4 µm 2 .
The experiments described below were performed with Jurkat TCs and the same authors determined the cell volume to V cell = 2 pl (12). This determines the values R cell = 8 µm and A cell = 804.2 µm 2 used in the present simulations. Given the cell radius R cell , the fractions f V and f R , as well as the factor f A , the surface to volume ratios required in equations (1) and (4) can be calculated by with

CRAC-CHANNEL
The open CRAC-channel current is determined by the electrochemical gradient This approach closely follows the model of Martin et al. (22). The validity of the Ohm's law approximation is only guaranteed within limited ranges of membrane potentials.

Single channel conductance
The single channel CRAC conductance was found to be extremely small in the order of g CRAC = 2 fs (61).

CRAC recruitment
The density of active CRAC-channels, estimated by the steady state CRAC-channel current, is a dynamic function of the ER-calcium concentration [(34) Figure 1C] described by where with C CRAC = 169 µM and n CRAC = 4.2. To our knowledge, this is the first time that the surface density of active CRAC-channels in the PM is modeled as a dynamic quantity. When estimating the same quantity from the degree of STIM1-redistribution, a rather similar relation is found with C CRAC = 187 µM and n CRAC = 3.8 [(34), Figure 2]. The uniformity of both curves supports the view that CRAC-channels are recruited and open in response to ER calcium depletion (34). It can be deduced from the similarity of both curves that the opening dynamics is substantially faster than the redistribution of STIM1. As no opening dynamics of the CRAC-channel is included in the model, the dynamics of the current itself and not of STIM1-redistribution is used.

CRAC density
In equation (25) ρ ± CRAC are the upper and lower limits of possible active CRAC densities. The resting value ρ CRAC,0 is not known from experiment and is determined by parameter fitting to calcium dynamics upon TCR stimulation (Section 2.10).
The density of CRAC-channels upon activation with PHA increased about 9-fold (33) which constraints ρ + CRAC . A 10-fold increase has been reported for the whole cell CRAC current in response to stepwise reduced C ER [ Figure 1C in Luik et al. (34)]. These findings translate into the condition where ρ + CRAC was determined by parameter fitting within the experimental boundaries in Section 2.10. A value for ρ − CRAC is not known and is determined by the steady state condition equation (36).

CRAC time scales
The time scale of CRAC recruitment can be estimated from the rising time of calcium curves which provides an upper bound of τ CRAC < 100 s for the activation time. It is likely that this time is associated with CRAC recruitment rather than with CRAC opening because opening time scales are typically much shorter. The time scale of CRAC recruitment is set to τ CRAC = 5 s. Larger values could also be used as the fit was insensitive to τ CRAC . Inactivation of CRAC-channels is difficult to be assessed (62). As the time scale of inactivation is in the order of 1000 s (14) and thus longer than the typical experimental durations used here, the present model ignores CRAC inactivation and assumes that the reduction of active CRAC-channels is a secondary effect of C ER recovery.

IP3R IN THE ER
The ER releases its calcium content if activated by IP3 and cytosolic calcium (63,64). The release of calcium from intracellular stores is based on the opening dynamics of RyR and IP3R in the membrane of the ER. TCs express both, RyR and IP3R and even different subtypes of them.
IP3Rs have binding sites for IP3 and calcium and exhibit complex forms of cooperativity (65). For the present purpose the heuristic description of IP3R activation and inhibition is sufficient. The characteristic feature of the IP3R conductance is a calciumdependent log-bell-shaped opening probability curve (63) which has been measured for ER vesicles from canine cerebellum and further reviewed in Foskett et al. (38). The open probability curve was best fitted by the product of an activating and an inhibiting Hill-function, both with Hill-coefficient 2 (63).

Open probability
TCR signaling leads to the generation of IP3, the ligand of the IP3R, which modulates the open probability of IP3R and is described as the product of an activation term g IP3R and an inactivation term h IP3R . The properties of single channel openings were quantitatively determined in Xenopus laevis oocytes (37, 38) and we assume that the single channel properties are transferable to TCs. The measured dynamics are well described by the previously published Mak-McBride-Foskett model (37, 38).
The dependencies of the model equation (27) on calcium and IP3 are depicted in Figure 2 and correctly reproduce the measurements in Mak et al. (37) suggesting that the modulating effect of IP3 is mediated by IP3R inactivation (38). At low calcium this effect is hardly visible and IP3R activation remains unaffected by changes in IP3 for resting concentrations beyond 50 nM. The dynamic IP3 range is between 100 nM and 1 µM (see (66), Table 1), a regime of IP3 at which the IP3R-type1 exhibits saturation (27,67). We do not aim at resolving whether the resting IP3 is lower in TCs or whether the IP3R characteristics are different in TCs, such that the DeYoung-Keizer model should be employed instead (27,65). It is assumed that the resting concentration of IP3 is in the range of 5-10 nM which ensures that an increased IP3-concentration has an impact on the IP3R opening probability as depicted in Figure 2.

IP3R calcium current
The steady state activation function (Figure 2) can be used to define the calcium current through the IP3R which follows the electrochemical gradient between the cytosol and the ER with g IP3R = 0.064pS. Note that the conductance differs between tissues (38). V − V ER is the potential difference between ER and cytosol, and V C,ER is the ER-reversal potential for calcium, calculated from the Nernst-equation equation (9). We assume an electrical equilibrated relation of cytosol and ER such that V = V ER holds true.

(In)activation dynamics
The activation and inactivation factors g and h are treated dynamically and approach equation (27) in steady state, while the adaptation of C IP3R,inh in equation (27) is treated in quasi steady state:

Activation time
Activation time scales can be determined from Mak and Foskett (68), Figure 5, and are in the range of less than 5 and 20 ms for depolarizations to 20 and 40 mV. Two exponentials were needed to fit the opening frequency. Using rat hepatocytes the activation and inactivation time scales were found to depend on the IP3-levels (69): activation varied between 100 and 500 ms for 10 µM and 400 nM of IP3, respectively (see Figure 1 therein). The time delay reported in Marchant and Taylor (69) is consistent with the IP3dependent time delay of channel opening of 1 s > τ IP3R > 100 ms using basophilic leukemia cells from rats (70). As the model focusses on calcium dynamics on the scale of minutes, a constant τ IP3R = 100 ms is assumed.

Inactivation time
Onset of inactivation happens in less than 2 min (68). A slow and a fast current were distinguished (69). The fast current inactivates on a time scale of 200-450 ms, the slow one between 1 and 6 s (see Figure 2 in the same publication). The authors attribute the fast time scale to inactivation of IP3R and the slow one to the depletion of the calcium content in the ER. Accordingly, only the fast time scales are relevant for the single IP3R, and θ IP3R = 300 ms is assumed.

Calmodulin dependence
It was found that the calcium-release from ER is reduced for high concentrations of the calcium-buffer calmodulin (71). Such a dependence is neglected in the present model.

IP3R density
The IP3R density on Jurkat TCs is not known and is determined using steady state condition equation (34).

PLASMA-MEMBRANE CALCIUM-ATPase
Plasma-membrane calcium-ATPase is an ATP-driven calcium pump which extrudes calcium from the cell to the ES. It was characterized in TCs (12). In a first attempt the dependence on the ATP concentration is ignored and assumed to be large enough in order to make the pump work optimally. In this case the activity is mainly dependent on the calcium concentration in the cell. A suitable modeling approach is with The current I PMCA is positive as it carries calcium out of the cell. The Hill-coefficient was determined to be n PMCA = 2 (72).

Turn-over rate
The turn-over rate of PMCA is approximately 30 Hz (73) which corresponds to an activity rate of k a = 0.03/ms which is also used in Sherman et al. (74). This turn-over rate can be translated into an electrical current by using that every pumping event corresponds to the flow of 2 electrical charges which yields I PMCA = z Ca ek a = 60 · 1.6 · 10 −19 C/s ≈ 10 −17 A = 10 −5 pA.

Calmodulin-dependent activation
Note that the values of half-activation also depend on the calmodulin concentration (75)(76)(77). For calmodulin concentrations above 1 µM (which is exclusively the case in all present simulations) full activation of all isoforms is ensured (75,78). Hence, the dependence on calmodulin is weak in this regime and is neglected.

Delay of activation
Binding of calcium to PMCA is a comparably fast process with a rate constant >3 per second (79). However, the activity of PMCA is delayed in some isoforms including the isoform 4b (76) which is relevant for TCs. The rate constant of PMCA activation upon stimulation with 500 nM of free calcium was in the range of 0.02 per second (76), which leads to a delay of PMCA activation in the range of τ PMCA = 50 s in equation (31). This delay was associated with a calmodulin and calcium-dependent activation (76). However, as the exact mechanism is not known this delay is modeled in equation (31) on a phenomenological basis.

PMCA density
For TCs no precise value of the protein density is known and the value is determined by the steady state condition equation (33).

SARCO/ENDOPLASMIC RETICULUM CALCIUM-ATPase
The calcium level in ER is kept high with the help of SERCA calcium pumps. The activity of SERCA is assumed to rapidly adapt to the present calcium concentration in the cytosol and can be described by a Hill-function.
In every turn-over cycle two calcium ions are transported per ATP (80). There are different subtypes of SERCA whose properties differ. In Jurkat TCs as well as in human tonsil lymphocytes the dominant isoform is SERCA2b (81).

Turn-over rate
The turn-over rates of most isoforms are in the range of k = 10 Hz (i.e., I SERCA = α SERCA z Ca ek = 6 · 10 −6 pA). For SERCA2b k 2b = 5 Hz was reported (82) which implies the value I SERCA = 3 · 10 −6 pA used in the model.

Calcium-dependent activation
For the SERCA isoforms 1, 2a, and 2b the half-activation C SERCA = 0.4 µM and the Hill-coefficient n SERCA = 2 are a good approximation [(82) Figure 4]. The half-activation of SERCA3 is around 1µM with the same Hill-coefficient (82). SERCA2b is an isoform active at relatively low calcium concentrations (82). Specifically, in Jurkat TCs as well as in human tonsil lymphocytes, the dominant isoform SERCA2b was characterized with n SERCA2b = 2.0 and C SERCA2b = 0.25 µM (81), which is used here.

SERCA density
Even though the expression of SERCA was shown to be modulated upon activation (83), the expression density of SERCA within ER is not known and is determined by parameter fitting in Section 2.10.

STEADY STATE DETERMINES PROTEIN DENSITIES
The resting state of the TC is determined by setting the dynamics in equations 1, 4, 7, and 24 to zero. Accordingly, the equations for C, C ER , P, and ρ CRAC , read: where I X,0 denotes the currents with all quantities X in the resting configuration X 0 . The parameters of the model are summarized in Table 1. www.frontiersin.org
As not all parameters could be determined by steady state conditions or by experimental constraints, Figure 1A in Bautista et al. (12) was used to determine the remaining free parameters. We used a two-step fitting procedure: at first, the differential evolution algorithm defined in Storn and Price (84) was incorporated into the C ++ -code of the model on the basis of all parameters in Table 1 that were not determined by steady state conditions. The parameters were varied within hard-coded boundaries dictated by experimental constraints (when available). The quality of the fit to the calcium data in Bautista et al. (12) was measured as the mean square deviation with X i and E i representing the simulation and experimental values, respectively. In a second step, the first approximative fit was subject to a sensitivity analysis in which each parameter was varied by 10% while monitoring the effect on QI. The three unknown protein densities ρ SERCA , ρ CRAC,0 , and ρ + CRAC , two sensitive IP3 related parameters P 0 and β P , as well as the very sensitive parameter V C,ER were used for fine-tuning the initial parameter fit with the same differential evolution Frontiers in Immunology | T Cell Biology algorithm. These fit parameters are marked variable in Table 1. The final fit reached QI = 0.100197 with N = 22. The sensitivity analysis was repeated for the final fit and the impact of each parameter on QI in equation (37) is provided in Table 1.
All subsequently described simulations are started with the cell in steady state as defined by the hard-coded equations (33)(34)(35)(36). Starting from these initial conditions, the respective stimulation protocols are applied as described in the results section.

RESULTS
In the methods section single protein characteristics were summarized and specific mathematical models capturing their main properties were proposed or cited. The models for the single proteins were combined to a whole cell model and the unknown parameters were determined using steady state conditions or by data fitting as described in Section 2.10. In this section, we replicate specific experimental setups described in the literature in silico and analyze the calcium dynamics from the perspective of the model.

TCR STIMULATION
The introduced TC-model is used to investigate the experiments of Bautista et al. (12) in Jurkat TCs. TC activation by stimulation of TCR induces an intracellular rise of second messengers like cADPR, NAADP, and IP3. In the present model this rise is collectively reflected in equation (7) for IP3. The IP3 signal activates IP3R and by this induces a calcium-release from the ER. The positive feedback loop of CICR leads to even more calcium-release from the ER, which in turn reduces the ER calcium concentration C ER . CRAC is activated in a C ER -dependent manner, as represented by equation (24). The rising cytosolic calcium is cleared by PMCA and the ER is refilled by SERCA, both being ATP-dependent processes.
Using 2 mM of external calcium a cell in resting state was activated with OKT3 via TCR [see Figure 1A in Bautista et al. (12)]. This induced a calcium peak to more than 1 µM within 50-100 s, which subsequently relaxed to a plateau level of ≈0.7 µM over the following 200 s. This behavior is well reproduced by the model (Figure 3A) and was used to determine the unknown parameters ( Table 1). The peak height relies on both, on a proper activation of IP3R and CRAC. The choice of P 0 turned out to be rather important in order to guarantee a proper activation of IP3R. The maximum activated CRAC density was essential for the height of the plateau. A value of f CRAC = 6.5 was found to correctly reproduce the plateau height as measured in Bautista et al. (12). Bautista et al. reported that the delay of PMCA activation is responsible for the calcium overshoot (12). Accordingly, we investigated whether this conclusion is supported by the model. For that purpose we reduced τ PMCA from 50 s (76) to 1 millisecond in the otherwise unaltered simulation. This modification does not touch the steady state configuration, such that all other parameters remained unchanged. The model still generated an overshoot of calcium, however, with a reduced amplitude (Figure 3A, blue dotted line). The height of the plateau as well as the relaxation time from peak to plateau remained unchanged. Thus, the simulation supports an influence of the PMCA delay onto the overshoot, but it turned out to not be essential for its existence. This surprising result led to the question whether a different delay in the model could explain the calcium overshoot. All delays were tested and the only delay with impact on the overshoot was IP3R inactivation, i.e., the parameter θ IP3R (Figure 3A, green dashed line). In conclusion, the model suggests that intrinsic properties of the IP3R and not the delay of PMCA activation are responsible for the calcium overshoot upon TCR stimulation. Next, the model is used for an analysis of the currents which are associated with the different phases of the calcium dynamics ( Figure 3B). The IP3-current (red full line) is clearly the largest and also persists beyond the overshoot of calcium due to SERCA activity (red dotted line). The calcium peak induces a PMCA current (black dotted line) that drives the calcium out of the cell (black dashed line). Thus, the CRAC current (black full line), induced by the loss of calcium in the ER, does not even induce a net flow of calcium into the cell (black dashed line) but just prevents the cell from running out of calcium. This model result suggests that the role of CRAC is the stabilization of the cell rather than its activation, which is mostly mediated by calcium from the ER.
The single protein currents ( Figure 3C) show that in the model IP3R currents react much more dynamic than CRAC currents. The increased cytosolic calcium even reduces the CRAC currents on the single channel level. However, the reduced ER calcium level induces a strong increase in the active CRAC-channel density ( Figure 3D). Thus, the increased whole cell CRAC current seen in Figure 3B is not a result of single channel responses but of a changed channel density.

TCR STIMULATION WITH ZERO EXTERNAL CALCIUM
Zero external calcium experiments aim at suppressing CRAC currents in order to investigate ER calcium currents in response to TCR stimulation. TCR stimulation of Jurkat TCs hold at zero external calcium results in an intracellular calcium peak after 50-100 s with reduced amplitude in comparison to stimulation with normal external calcium [ Figure 2A in Bautista et al. (12)]. The increased calcium is cleared below baseline level within 100-200 s.
In the model, the Nernst-equation prohibits the use of zero external calcium conditions. At very low calcium concentrations the Nernst-equation loses its validity and the reversal potential diverges. Therefore, external calcium is set to the concentration C * ext at which the CRAC current vanishes in resting state: This mimicks the suppression of CRAC currents as intended in the experiment. The in silico result is shown in Figure 4. The CRAC current vanishes at resting state ( Figure 4B, black full line). This is approximately true during the whole experiment, such that the method for mimicking the zero external calcium experiment appears appropriate.
The calcium peak ( Figure 4A) is lower than the one found in Figure 3 which is qualitatively consistent with the experimental result (12). Furthermore, the time scales of calcium rise and clearance are perfectly matched between simulation and experiment. Even the overshoot of clearance below the baseline calcium level is fully reproduced. However, quantitatively, the peak is higher in www.frontiersin.org

A B
FIGURE 4 | Calcium dynamics in response to TCR stimulation at zero external calcium. Simulation of the experiment in Figure 2A in Bautista et al. (12). (A) The cell is in steady state until t = 10 s, when stimulation is started by setting T (t ) = 1.6 in equation (7) and external calcium is set to equation (38).
Stimulation and external calcium are kept constant throughout the simulation.
(B) Whole cell currents show the almost complete inhibition of CRAC currents. The calcium peak (left) is generated by ER calcium only. Negative currents are calcium currents into the cytosol. theory than in experiment. As the in silico peak is exclusively generated by the ER, one might hypothesize that the ER is too big or its calcium content is too high. As these parameters were already chosen comparably low (Table 1), it is more likely that the lack of external calcium influences the stimulation of the cell. The peak size can be reduced to the measured amplitude by reducing the stimulation T from 1.6 to 1.25 (not shown).

BLOCK OF SERCA AT ZERO EXTERNAL CALCIUM
Thapsigargin (TG) is frequently used to block SERCA activity as it prevents calcium uptake by the ER and leads to a continuous reduction of ER calcium. As low ER calcium recruits and activates CRAC-channels, this would lead to a strong influx of extracellular calcium. In order to prevent this according ER depletion and CRAC activation experiments are performed in zero calcium medium. This strategy was used in Jurkat TCs to generate a cell state in which the ER is mostly void of calcium and CRAC-channels are recruited and activated to a maximum (12,13,25,86). It was reported that this procedure leads to a transient calcium peak of about 0.5 µM after more than 100 s which is slowly cleared and reaches calcium levels below the resting level [ Figure 1A in Quintana et al. (86)].
Having established a strategy [equation (38)] for mimicking a medium with zero calcium, a SERCA block is performed in silico by setting I SERCA = 0 at t = 10 s. No TCR stimulation was applied. The measured dynamics are well reproduced without any further parameter fitting (Figure 5A, black full line). Furthermore, the intended depletion of ER is achieved (Figure 5A, red dashed line): a continuous reduction of ER calcium is observed. Note that also IP3 exhibits some dynamics (Figure 5A, blue dotted line) which further accelerates ER calcium loss by activation of IP3R. As expected, the reduced ER calcium leads to the recruitment of active CRAC-channels ( Figure 5B).

BLOCK OF PMCA IN A TG TREATED TC
As the TG-mediated SERCA block works in silico (Figure 5), the role of PMCA in the clearance of cytosolic calcium is investigated.
Following Figure 6 in Bautista et al. (12) TG was applied in a zero calcium medium as in Figure 5. Then a pulse of 2 mM external calcium was applied for 50 s. This induces a steep rise in calcium which is also steeply cleared again. Together with Lanthan (La 3+ ), a PMCA, and CRAC inhibitor (2), the clearance of such a calcium peak was substantially slower (12). A block of PMCA alone by carboxyeosin led to an only weakly modified clearance time, without a return to baseline levels within 300 s.
The same protocol is applied in the model (Figure 6). As before, ER calcium is depleted (Figure 6A, red dashed line) and the active CRAC density is increased correspondingly (as in Figure 5B). The initial free cytosolic calcium peak (Figure 6A, black full line) is the same as in Figure 5A. Upon increasing external calcium to 2 mM for 50 s at t = 300 s, a strong CRAC current is induced (Figure 6B, black full line), which steeply increases cytosolic calcium ( Figure 6A, black full line). This calcium is also quickly cleared upon return to the mimicked zero calcium medium. Calcium clearance is dominated by the PMCA current ( Figure 6B, black dotted line). However, in the model a small CRAC current is observed supporting extrusion of calcium out of the cell in the case of zero external calcium concentration. Upon permanent block of PMCA and repetition of the transient stimulation by external calcium, it is this backward CRAC current that clears calcium from the cytosol. The time course of the clearance is slower than without PMCA block and the calcium baseline is not reached after 350 s ( Figure 6B, black full line). This is a similar behavior as in the corresponding experiment with carboxyeosin [ Figure 6D in Bautista et al. (12)]. However, it is not known whether the real CRAC allows for such inverse current under zero calcium conditions. The return of calcium to the baseline might also be supported by uptake of calcium by other organelles like mitochondria, which is not covered by the present model.
In the case of PMCA-and CRAC-block with La 3+ a rather slow calcium clearance is observed in experiments which apply the same stimulation protocol [ Figure 6C in Bautista et al. (12)]. In silico no clearance is observed at A B The cell is in steady state until t = 10 s, when SERCA block is applied and external calcium is set to equation (38). At t = 300 s 2 mM external calcium are applied for 50 s followed by a switch back to equation (38). The procedure is repeated at t = 600 s. all. SERCA is blocked and positive I IP3R -currents are not allowed in the model, such that an uptake of calcium into the ER is excluded. A full block of PMCA and CRAC also excludes any expulsion of calcium out of the cell. Thus, the model suggests that the slow clearance of cytosolic calcium is either due to incomplete block of PMCA or to leakage currents.

BLOCK OF PMCA IN AN UNTREATED TC
The model suggests that the role of CRAC for TC activation is mainly the maintenance of the integrity of the TCs during stimulation in the sense that it prevents the activated TC from running out of calcium. If we block PMCA in silico at the time of TCR stimulation (protocol as in Figure 3) in an otherwise untreated TC, the TC would be prevented of losing calcium. According to our interpretation of the role of CRAC we would expect that CRAC activity is strongly reduced in comparison to Figure 3. T lymphocytes receptor stimulation of PMCA-blocked TCs is predicted to induce a strong increase of cytosolic calcium (Figure 7, full black line). Thereby, the steady state CRAC current is not increased (as in Figure 3B, full black line) but reduced (not shown). However, the CRAC current is not reduced to zero such that the total block of PMCA infers a persistent net influx of calcium into the cell and, thus, to a persistent increase of cytosolic calcium. This would ultimately destroy a real TC. As SERCA activity is normal in this in silico experiment, calcium in the ER is only transiently reduced (not shown), leading to a transient and weak increase in the CRAC density (Figure 7, dotted red line). The amplitude is less than twofold instead of sixfold in Figure 3D.  Figure 3 is repeated. At the time of stimulation (t = 10 s), PMCA activity is blocked preventing calcium efflux from the cell. Cytosolic calcium (full black line, left axis) and the response of the CRAC density (dotted red line, right axis) are shown. The axis for the CRAC density was scaled as in Figure 3D for better comparison.
This result supports the interpretation of the role of CRACchannels in silico and may be tested in experiment. It further shows, that the calcium level in the ER may be used as an indicator for the overall calcium status of TCs.

DISCUSSION
Within the presented investigation whole cell calcium dynamics were derived from models of single transmembrane ion-conducting proteins. The model successfully described a number of experimental settings and captured important characteristics of calcium dynamics upon TCR stimulation. In particular, the role of store-derived calciumrelease versus CRAC activation was well represented. We conclude, that we have generated a modeling framework suitable for the quantitative analysis of calcium dynamics of TCs.
The value of using measured single transmembrane protein characteristics is twofold: at first, it substantially reduces the number of free parameters in an otherwise very complex model. All measured single protein properties were just implemented and not altered in the fine-tuning of the model. The reduced number of free parameters increases the predictive power of the model. Secondly, it can be considered as a multi-scale approach which links single proteins to whole cell behavior. This allows the analysis of the dynamics on the single protein level and its implications onto the cellular properties. The drawback of this approach is that not all parameters could be derived from TCs such that we had to assume that single transmembrane properties are universal, which is a correct assumption in many cases (30). Within this framework, cell-specific properties are controlled by the protein densities and the cell-specific properties, both determining the activity range of the respective proteins. However, the predictive power of the model would benefit from corresponding TC related data.
The model of Jurkat TC calcium dynamics as supported by the mathematical model starts from a TCR-derived increase in second messengers like IP3. This ultimately activates IP3R currents and induces an initial calcium current from the ER into the cytosol which triggers CICR. However, calcium uptake via SERCA activity equilibrates the calcium loss in the ER, which leads to a zero net flux on long-term -despite ongoing TCR stimulation. The ER calcium loss, according to the model, was the major contribution to the free cytosolic calcium peak. However, as PMCA activity leads to an overall loss of calcium in the cell, a compensation mechanism is required for a sustained elevation of free cytosolic calcium. The model suggests that the CRAC activity essentially contributes to this compensation mechanism. The initiation of CRAC currents by the depletion of ER calcium levels is in line with this interpretation. The model predicts, that a standard stimulation of TCs together with a block of PMCA activity would lead to a strong calcium rise together with a minor and transient increase of the CRAC density and on long-term to a reduced CRAC current (see Figure 7). This prediction may be tested in experiment in order to validate this interpretation of the role of CRAC-channels.
The emerging hypothesis that the role of CRAC is the stabilization of the TC calcium level, needs to be further strengthened by more detailed modeling work. In particular, early events after stimulation like NAADP generation (87) and subsequent RyR calcium currents from the ER are essential for the early calcium rise (11,18) and will have to be included in the mathematical model for a proper time-resolved coverage of calcium dynamics. The need for additional mechanisms is also underpinned by the strong sensitivity of the model behavior to changes in the IP3 resting concentration P 0 ( Table 1). In the present simulation the long-term calcium plateau height mostly relies on the maximum possible CRAC activation by store-operated calcium depletion. For the long-lasting calcium rise, cADPR was proven relevant (16,40) and has to be considered in the context of the present hypothesis on the role of CRAC. Also the transferability to human blood derived TCs, which exhibit a different geometry, has to be assessed.
As the PMCA block experiment at zero external calcium has shown, it might be important to include leakage currents into the model. However, it should be noted that the postulated role of delayed PMCA activity for the calcium overshoot after TCR stimulation (12) could only partially be confirmed by the mathematical model. With fast activation of PMCA, an overshoot was still observed in our model, and the overshoot could only be suppressed by lack of IP3R inactivation. It should be noted that this result might rely on the Mak-McBride-Foskett model (37), which was used for IP3R dynamics and which exhibits inactivation at rather low IP3. The result might differ if the TC calcium dynamics would be based on the DeYoung-Keizer model for IP3R activity (65).
The proposed model has limitations in its range of applicability. For example, the usage of the Nernst-equation for the chemical gradient and Ohm's law for the current-voltage relationship of ion-conducting pores is justified only in narrow limits. Experiments with zero external calcium drive the model to the very limits of this range of applicability which was circumvented here using a phenomenological approximation. The model may be reformulated in terms of the Fokker-Planck-equation in order to describe ion transport through the pores in more detail. Furthermore, it is known that calcium entrance points lead to spatially inhomogeneous calcium dynamics (22) which are not covered by the present space-averaged model. The value of the present approach lies in the surprising result that quantitative characteristics of single transmembrane proteins are sufficient to determine the cell behavior in the framework of an ordinary-differentialequation based model. The model has proven its predictive power, as it was fitted to data of one experiment in Figure 3 and could be used to predict and explain further data of calcium dynamics generated under other experimental conditions in Figures 4  and 6. It is planned to elaborate the potential and the limits of the model predictions by application to further experimental settings.

ACKNOWLEDGMENTS
We thank Dr. Harald Kempf for revising the manuscript. Christine Schmeitz was supported by the Helmholtz International Graduate School for Infection Research. This study was partially supported by the Deutsche Forschungsgemeinschaft (grants GU 360/13-1 and GU 360/15-1 to Andreas H. Guse). Esteban Vargas and Michael Meyer-Hermann were supported by the BMBF within the GerontoSys initiative. Michael Meyer-Hermann was supported by HFSP. www.frontiersin.org