A Model for the Coexistence of Competing Mechanisms for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ca}^{\text {2}+}$$\end{document}Ca2+ Oscillations in T-lymphocytes

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ca}^{\text {2}+}$$\end{document}Ca2+ is a ubiquitous signaling mechanism across different cell types. In T-cells, it is associated with cytokine production and immune function. Benson et al. have shown the coexistence of competing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ca}^{\text {2}+}$$\end{document}Ca2+ oscillations during antigen stimulation of T-cell receptors, depending on the presence of extracellular \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ca}^{\text {2}+}$$\end{document}Ca2+ influx through the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ca}^{\text {2}+}$$\end{document}Ca2+ release-activated \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ca}^{\text {2}+}$$\end{document}Ca2+ channel (Benson in J Biol Chem 29:105310, 2023). In this paper, we construct a mathematical model consisting of five ordinary differential equations and analyze the relationship between the competing oscillatory mechanisms.. We perform bifurcation analysis on two versions of our model, corresponding to the two oscillatory types, to find the defining characteristics of these two families.


Introduction
The role of calcium ions (Ca 2+ ) as an intracellular second messenger in cells has been studied since the late 1800s (Ringer 1883;Carafoli 2002).Regulation of the Ca 2+ concentration in the cytosol is known to control a wide array of cellular functions across a diverse range of cell types.As this regulation is central to cellular function, cells are equipped with a complex set of pumps, channels and exchangers designed to maintain a tight control of the cytosolic Ca 2+ concentration.Ca 2+ extrusion from the cytoplasm is done primarily through plasma membrane Ca 2+ ATPases (PMCA) and sarco/endoplasmic reticulum ATPases (SERCA) (Higgins et al. 2006;Michelangeli and East 2011;Wuytack et al. 2002).These pumps require ATP to transport Ca 2+ across their respective membranes and establish a steep concentration gradient.Ca 2+ entry into the cytoplasm occurs via Ca 2+ channels that allow the ions to flow down the electrochemical gradient.At the endoplasmic reticulum (ER) membrane, the two main Ca 2+ release channels are the inositol 1,4,5-trisphosphate receptor (IP 3 R) and the ryanodine receptor (RyR) (Foskett et al. 2007;Kevin Foskett and Daniel Mak 2010); at the outer plasma membrane (PM), Ca 2+ channels are classified as receptoroperated (ROCC), store-operated (SOCC) or voltage-gated (VGCC), depending on their activation mechanism (Dupont et al. 2016;Thompson and Shuttleworth 2013;Catterall 2011;Yoast et al. 2020).
In T-lymphocytes, agonist stimulation by α-CD3 induces two distinct types of Ca 2+ oscillations (Benson et al. 2023).When stimulated in a Ca 2+ -rich environment, the cytoplasmic Ca 2+ concentration in these cells exhibits transient Ca 2+ spikes before transitioning to a wide, 'sinusoidal-like' oscillation (Fig. 1A).To sustain these oscillations, extracellular Ca 2+ is transported into the cytoplasm through the store-operated Ca 2+ Entry (SOCE) pathway via a SOCC located at the junction between the PM and ER membrane: the Ca 2+ release-activated Ca 2+ (CRAC) channel (Trebak and Kinet 2019;Dupont et al. 2016;Dolmetsch and Lewis 1994).
This Ca 2+ transport through CRAC has physiological relevance: in T-lymphocytes, translocation of the nuclear factor of activated T cells (NFAT) proteins from the cytoplasm to the nucleus, for example, is a pivotal step for cytokine production and correct immune function.Ca 2+ transported through this SOCC binds with the calmodulin protein to trigger the translocation of NFAT (Trebak and Kinet 2019;Srikanth and Gwack 2013).Hence, models of Ca 2+ oscillations in lymphocytes are usually concerned with CRAC activity above that of other channels (Dolmetsch and Lewis 1994;Naik 2020;Xf et al. 2008;Lewis 2011).
Replication of the same experiment in a Ca 2+ -free environment, however, can lead to the appearance of a different class of oscillatory patterns including narrow Ca 2+ spikes (Fig. 1B), wide spikes (Fig. 1C) and wide spikes with intermittent bursting on a raised plateau (Fig. 1D) (Benson et al. 2023).The presence of Ca 2+ spikes in this scenario suggests that, in the absence of SOCE, a secondary oscillatory mechanism arises in this cell type, one related to fast Ca 2+ release from the ER, but that is typically 'silenced' by CRAC activity.Additionally, prior models of Ca 2+ oscillations have linked the appearance of Ca 2+ bursts and wide spikes to IP 3 production and decay and its corresponding effect on IP 3 R activity (Cloete et al. 2021).These experimental results led Benson et al. (2023) to suggest the existence of two fundamentally different kinds of oscillations during Ca 2+ signalling: IP 3 Rmediated and CRAC-mediated.The IP 3 R-mediated oscillations correspond to those observed in the Ca 2+ -free environments, and the CRAC-mediated ones to the Ca 2+rich environment.They supported this conclusion with a mathematical model capable of reproducing these two oscillatory types.
In this paper, we perform a more detailed study of the model presented in Benson et al. (2023) and show that it can reproduce the CRAC-mediated oscillations when extracellular Ca 2+ transport is considered, and the multiple types of IP 3 R-mediated spikes when Ca 2+ influx is greatly decreased or removed.We show that the CRACmediated oscillations are sustained due to delayed Ca 2+ influx via CRAC, that the CRAC channel has a steep activation curve dependent on ER Ca 2+ concentration and that CRAC activity suppresses the IP 3 R-based Ca 2+ spikes when present.After that, we show that the model is capable of reproducing all three types of IP 3 R-mediated Ca 2+ spikes (narrow, wide, wide with an oscillatory plateau) depending on the level of IP 3 R stimulation as well as the level of the ER Ca 2+ stores.receptors (TCR) of this cell type is linked to increased IP 3 R activity (Trebak and Kinet 2019;Benson et al. 2023).Additionally, we assume that cytoplasmic Ca 2+ increases the rate of IP 3 degradation (Dupont et al. 2016).Finally, we include the effect of Ca 2+ extrusion from the cytoplasm via SERCA and PMCA@.
Let c and c e , respectively, denote the Ca 2+ concentration in the cytoplasm and in the ER@; h is an auxiliary variable that controls the rate of Ca 2+ -dependent inactivation of the IP 3 R; p denotes the IP 3 concentration in the cytoplasm; s denotes the rate of Ca 2+ influx through the CRAC channel.The schematic version of the model is shown in Fig. 2. The model equations are ) Here J IP 3 R is the Ca 2+ flux out of the ER through the IP 3 Rs, J SERCA is the flux through a bidirectional SERCA pump, and δ J PM is the flux through a unidirectional PMCA pump.The expressions for J IP 3 R , J SERCA , h ∞ and J PM are taken directly from Sneyd et al. (2017) and are reproduced in Appendix A, where default parameter values are also provided.A detailed construction of the IP 3 R model can be found in Dupont et al. (2016).The details of the J deg term are discussed in Sect.2.1.The parameter γ in (1b) is the volume ratio between the cytoplasm and the ER@.To enable easy simulation of the Ca 2+ -free scenario, we include an additional parameter δ in front of the s and J PM terms.δ is a measure of the rate of extracellular Ca 2+ exchange relative to the rate of the ER fluxes (Sneyd et al. 2004).The case when δ > 0 is referred to as the open-cell model.When δ = 0, Ca 2+ transport through the PM is not allowed and the model exhibits one conserved quantity: the total amount of free Ca 2+ inside the cell, given by C t = c + c e /γ (Dupont et al. 2016;Sneyd et al. 2004).We exploit this observation by rewriting c e as a function of c and reducing the dimension of our system by one.By setting δ = 0, s decouples from the rest of the system, further reducing the dimension of our system by one, giving us a set of three ODEs given by This version of the model is referred to as the closed-cell model.Since c e is not an explicit variable in this version the functional form of the J IP 3 R and J SERCA fluxes is different to their open-cell counterpart; the closed-cell version of the fluxes and corresponding parameter values are given in Appendix A.

IP 3 Dynamics
Binding of an agonist to a G protein-coupled receptor (GPCR) leads to the activation of a G protein, which in turn activates phospholipase C (PLC).PLC then acts on phosphatidylinositol bisphosphate (PIP 2 ), cleaving it and producing both diacylglycerol (DAG) and IP 3 .While PLC activity has been shown to depend on c in some cell types, this is not always the case, and for now we assume that the maximal rate of PLC activity (V PLC ) is determined solely by the level of agonist stimulation (Lee et al. 2023).
On the other hand, IP 3 can be metabolized into inositol 1,3,4,5-tetrakisphosphate (Ins (1,3,4,5) P 4 ) by a Ca 2+ -dependent 3-kinase (Dupont et al. 2016).The activity of this 3-kinase can be enhanced by the binding of the Ca 2+ /calmodulin complex.Through this pathway, increases in cytoplasmic Ca 2+ concentration can then lead to faster metabolism (or degradation) of IP 3 .This degradation is modelled by where V deg represents the maximal degradation rate.Finally, we consider that changes in IP 3 concentration occur on a timescale τ p .
Experimental trials show that lymphocytes with genetic mutations that prevented the expression of both the STIM1 and STIM2 proteins, and thus were unable to activate the CRAC channel, did not show signs of ER refilling (or CRAC activity) under conditions of maximal ER depletion induced by blockage of SERCA pumps with thapsigargin (Tg), even in Ca 2+ rich environments.Cells that had been modified to exhibit only one variant of STIM1 or STIM2 did show a CRAC current and ER refilling, even if reduced when compared to the wildtype (WT) cells.Given this evidence, we consider that SOCE through the STIM-Orai channel is the only relevant pathway for extracellular Ca 2+ influx (Benson et al. 2023).
We model Ca 2+ transport through this ER-dependent channel using the phenomenological expression which is a sigmoidal function in c e .Measurements of the CRAC current for T cells as a function of the ER Ca 2+ concentration result in a curve with a similar form (Luik et al. 2008).The parameter s 1 controls the steepness of the activation curve, while K e controls the midpoint.The maximal rate of transport through the channel is given by V SOCE .We do not make a distinction between the contribution of STIM1 and STIM2 at this stage, and address this in Sect. 4. Forced ER depletion via the use of Tg has been shown to induce oscillations (Dolmetsch and Lewis 1994;Benson et al. 2023).At the same time, experimental measurements have shown that increases in CRAC activity lag behind ER depletion, which suggests that the formation of the CRAC channel is slow in comparison to the ER depletion rate (Dolmetsch and Lewis 1994;Zweifach and Lewis 1995;Lewis 2011).This led Dolmetsch and Lewis to hypothesize that Ca 2+ influx via the CRAC channel and Ca 2+ release from the ER oscillate out of phase with each other (Dolmetsch and Lewis 1994).To model this phenomenon, we consider that s evolves on a timescale given by τ s .This timescale must be slow relative to the changes in the Ca 2+ variables, to generate the delayed CRAC response.1.All simulations use initial values corresponding to the equilibrium solution for Parameter values not specified here can be found in Table 1 3 Results

The Model is Capable of Reproducing the Two Families of Oscillatory Patterns
We implement the open and closed-cell versions of (1) using XPPAUT (Ermentrout 2012) and the AUTO package (Doedel et al. 2021).Data processing and graph generation are done using the numpy 1.19.1 and matplotlib 3.4.0Python packages (Harris et al. 2020;Hunter 2007).The code used to generate the figures can be found in github (0Huitzil 2023).Using the parameter values found in Table 1, model (1

Bifurcation Analysis of the Open-cell Model
Figure 4 shows the bifurcation diagram and representative oscillatory traces obtained from simulation of the open-cell model (1), using V PLC as a bifurcation  parameter, for a fixed value δ = 2.The CRAC-mediated oscillations appear on the diagram as a branch of stable periodic orbits (green) situated between two Hopf bifurcations (HB).The distance between the Hopf bifurcations can be increased or decreased by modifying the K p parameter accordingly; this change does not affect the c values during oscillation, but does increase or decrease the mean IP 3 concentration, respectively.
For V PLC = 0.1, the oscillation has a period of approximately 40 seconds and has a sinusoidal shape, with a period of slow increase in c followed by a period of slow decrease, both lasting approximately 20 seconds each.The maximal loss of ER Ca 2+ during a single oscillation is approximately 20 μM.This amount of ER depletion causes a constant, but moderate, formation of the CRAC channel, which causes s to oscillate between 25% to 50% of its maximal rate.The value of c oscillates from 0.1 to 0.15 approximately; the amplitude of this oscillation can be increased by rescaling c, but due to its involvement in most of the fluxes this would necessitate a substantial rescaling of the values in Table 1.

CRAC-mediated Oscillations Require a Steep and Delayed SOCE Mechanism
The oscillation shown in Fig. 4B (V PLC = 0.1) uses a value of s 1 = 0.2.This choice of s 1 was made to ensure that J SOCE was steep enough as a function of c e for the CRAC-mediated oscillations to exist.This value can be decreased to ensure a shallower activation.To investigate the required steepness of this curve, we performed parameter continuation over s 1 on both the steady state and the stable periodic orbit observed with V PLC = 0.1 (Fig. 5A).The branch of CRAC-mediated oscillations exists for s 1 ≥ 0.18.Similar diagrams can be obtained for close values of V PLC , but the exact position of the HB seen in Panel A will change.We performed two parameter continuation on this HB point and plot the corresponding curve in s 1 − V PLC parameter space (Fig. 5B).This HB curve divides the plane into two regions, one where CRAC-mediated oscillations are observed (shaded green) and one where no oscillations are observed.Continuation reveals a minimum value of s 1 needed to sustain the CRAC-mediated oscillations, given by s 1 ≈ 0.15.This minimum value suggests that the formation of the CRAC channel in these cells follows a very steep 'trigger-like' activation curve.For comparison, our SOCE activation curve has a similar steepness as a reverse Hill function with a coefficient n = 100 (Fig. 5E).
A similar analysis is performed on the parameter τ s (Fig. 5C, D).The time series shown in Fig. 3A uses a value of τ s = 15.CRAC-mediated oscillations are observed for a wide range of parameter values but, crucially, the diagram in Panel C shows the branch of stable periodic orbits ending at a HB at both high and low τ s values.This observation suggests that, although a delayed SOCE activation is necessary to exhibit CRAC-mediated oscillations, this delay must be bounded for the oscillations exist.Panel D shows the curve obtained by following the two HBs in the τ s − V PLC parameter space.This curve once again divides the parameter plane into a region with CRAC-mediated oscillations and one with no response.This region is bounded, in contrast to the region obtained in the s 1 − V PLC scenario, which further emphasizes the upper and lower bounds on the delay needed to observe the CRAC-mediated oscillations.
These two observations are not independent of each other; it could be argued that the same result could be generated by using a less steep curve alongside a faster formation timescale.To investigate this, we followed the HB point in the s 1 − τ s parameter space for a fixed value of V PLC = 0.1 (Fig. 5F).While CRAC-mediated oscillations are present for a region with a less steep activation curve and a faster formation of the CRAC channel, the HB curve defines a positive boundary on both parameters, and no oscillations are observed for low values of s 1 and τ s , indicating that both features of our model are necessary for the existence of the CRAC-mediated oscillations.1) for V PLC = 0.1 and δ = 2. Panel C shows a similar diagram using τ s as a bifurcation parameter.In both panels, the black curve represents the equilibrium solution, while the green solid lines show the maximum and minimum c values of a branch of stable periodic orbits.Panels B, D and F show the locations in the s 1 − V PLC , τ s − V PLC and s 1 − τ s (respectively) parameter planes of the Hopf bifurcations (black curves) and regions of CRAC-mediated oscillations (green-shaded regions).Panel E shows the SOCE activation curve with a steepness value s 1 = 0.2 in green, and a reverse Hill curve with a corresponding coefficient n = 100.Parameter values not specified here can be found in table 1 (Color figure online)

IP 3 R based Ca 2+ Spikes are Suppressed by CRAC Activity
For the CRAC channel to form and allow the influx of extracellular Ca 2+ the ER must itself be depleted enough to activate the STIM proteins and create the STIM-Orai complex.In our model, this Ca 2+ release from the ER must necessarily be done through the IP 3 R.The literature on the activity of the IP 3 R and its regulation by both Ca 2+ and IP 3 is extensive and a thorough review on the subject can be found in Foskett et al. (2007).
In our model, the release of Ca 2+ from the ER is controlled by the open probability of the channel (P 0 ).Once open, the Ca 2+ released from the channel activates the IP 3 R in a process known as Ca 2+ -induced Ca 2+ release (CICR), which is the process responsible for the fast upswing in concentration observed in Ca 2+ spikes (Fig. 3B-D) (Dupont et al. 2016;Sneyd et al. 2017;Jelbart et al. 2022).After opening, the receptors are inactivated and remain closed for a time.The variable h represents the rate of activation of the IP 3 R in response to the binding of Ca 2+ ; the rate at which this activation occurs is regulated by c via modulation of the timescale τ h (1c).Additionally, it is worth remarking that increases in c also promote degradation of IP 3 (1d).
Inspection of the time series for the P 0 term during the CRAC-mediated oscillation reveals that no fast CICR is observed during this scenario, in contrast with the similar time series for a narrow spike (Fig. 6).The mathematical expression for the P 0 term can be found in Appendix A. Instead of a fast opening and closing of the channel, the P 0 curve indicates that the IP 3 R channel remains open for a longer time which allows for a gradual Ca 2+ leak from the ER on the timescale of the CRAC dynamics.This gradual leak allows for greater depletion of the ER stores during CRAC-mediated oscillations when compared to the narrow spikes (Fig. 10B-D).
The effect of the suppression of the fast CICR by CRAC activity can be captured by simulating the open-cell model (1) under two different sets of initial conditions, a 'low' ER load (c e = 809) such that CRAC is partially activated at the start of the simulation (Fig. 7A-B), and a 'high' ER load (c e = 850).The initial conditions for c and s are obtained by solving equations 1b and 1e assuming an unstimulated cell (V PLC = 0).Under these conditions, CRAC is completely inactivated at the start of the simulation in the high ER load scenario (Fig. 7C-D).The time series of the cytoplasmic Ca 2+ in the low ER load does not exhibit transient spiking as it converges to the limit cycle.The relatively large value of s positions this initial state close to the limit cycle (Fig. 7B).The figure also shows the J SOCE curve as a function of c e .As the c e load decreases and crosses the midpoint of the J SOCE curve (K e = 800 μM) it triggers the activation of CRAC which increases the value of s.Superimposing this trajectory on the J SOCE activation curve shows c e oscillating between 790 and 820 μM, which correspond to the complete activation (J SOCE = 3) and complete deactivation (J SOCE = 0) values of the curve, respectively.However, due to the slow evolution of s, the actual flux though the channel never reaches either of these values.
In the high ER load scenario, the time series for c shows transient Ca 2+ spikes before converging to the limit cycle.These spikes are less notable when looking at the c e − s plane, due to the relatively low amount of ER Ca 2+ being released during this type of oscillation.The spikes persist for approximately 30 s, a time frame similar to the formation of the CRAC channel, before being replaced by the CRAC-mediated oscillations once the channel has been formed.These narrow spikes are only observed at values of c e close to the midpoint of the J SOCE curve (K e ).Using the expression for total Ca 2+ , C t = c+c e /γ , we can notice that these transient spikes appear at C t values between 140 and 150 μM.Narrow spikes are the only type of oscillation observed in the closed-cell model for this range of C t values (Sect.3.3).Thus, the model suggests that restricting extracellular exchange during a CRAC-mediated oscillation could result in the appearance of narrow spikes, but not wide spikes.

Bifurcation Analysis of the Closed-Cell Model
To study the narrow-spike oscillations generated by the closed-cell model (2), we compute bifurcation diagrams with V PLC as the bifurcation parameter and various choices of fixed C t .The case for C t = 140 is shown in Fig. 8A.Here, there are two supercritical Hopf bifurcations with a branch of stable periodic solutions between these bifurcations.Corresponding time series for c and p for a representative value of V PLC are shown in Fig. 8B.
When C t = 95 the two HBs are subcritical, and four new bifurcation points are detected: two saddle-nodes of periodic orbits (SNPO) and two periodic doubling (PD) bifurcations.The branch of periodic solutions is now only stable for parameter values between the SNPO1 and the PD1 bifurcations (Fig. 8C) A representative time series corresponding to the stable branch of periodic solutions is shown in the upper part of Panel D. Secondary branches of periodic solutions are created from the PDs, but they are not plotted here.The region between the PD2 and the HB2 points has no attractors on either the equilibrium line or the branch of periodic orbits.Instead, the corresponding attractor occurs on a branch of periodic solutions that is isolated from the rest of the diagram and is not plotted here.A representative time series for this region can be seen in the lower part of Panel D, which corresponds to the wide spike with an oscillatory plateau shown in Fig. 3.
The diagram obtained at C t = 75 shares all the qualitative features of the C t = 95 case.However, the stable region of the main periodic orbit branch has shrunk.The top part of Panel F shows a representative trace from the stable region of the main branch of periodic orbits, between the SNPO1 and PD1 points of Panel E@; the lower part shows an attractor belonging to a different branch of stable periodic orbits that emanates from the main branch at the PD2 point.
The information about the different types of oscillations produced by the closed-cell model ( 2) can be summarized by looking at the two parameter bifurcation set shown in Fig. 9. Narrow spike oscillations can be found in the region bounded by the HB, PD and TR curves (Region I).In reality the left boundary is set by a curve of SNPO that lies very close to the HB curve, but it is not drawn here for simplicity.
Delimiting the regions in parameter space where the wide spikes and wide spikes with plateau exist is a harder task.The PD curve in Fig. 9 can be used as proxy for the boundary of the region containing wide spikes (Region II) due to it separating the region where the narrow spikes are stable, but it must be remarked that this curve is not the actual boundary.We have found wide spikes with oscillatory plateaus in the area bounded by the PD, TR and SNPO curves (Region III).It is unlikely that these curves are the exact boundaries between Regions I, II and III, but they are close to the boundaries and provide useful proxies for distinguishing the regions in which the two types of wide spikes are observed.

Qualitative Differences Between the Two Families of Oscillations
For δ ≈ 2, the open-cell model (1) exhibits CRAC-mediated oscillations.In contrast, for δ ≈ 0 the model exhibits IP 3 R mediated Ca 2+ spikes.These two behaviours are shown in Panels B and D of Fig. 10, where we compare the attractor obtained for δ = 2 against the corresponding attractor for δ = 0.01.The CRAC-mediated oscillation has a longer period than the IP 3 R mediated spike, as well as a wide variation in ER concentration values, needed to trigger the formation of the CRAC channel.  1) using V PLC for a fixed value of δ = 2. Panel B shows the time series obtained by stimulating the model with V PLC = 0.1, the cytoplasmic Ca 2+ concentration is plotted as a solid curve, and the ER concentration as a dashed curve.Panels C and D show analogous figures, but with a fixed parameter value of δ = 0.01.Parameter values not specified here can be found in Table 1 Bifurcation diagrams for examples of these two cases are shown in Panels A and C of the same figure.While the overall structure of the two bifurcation diagrams is the same, the time series show there are significant differences in the behaviour of the attracting periodic solutions.
To explain the transition between the CRAC-mediated and the IP 3 R-mediated oscillations in the open-cell model (1) we treat δ as a bifurcation parameter and perform two parameter continuation on the HB points found on both the δ = 2 and the δ = 0.01 diagrams (Panels A and C of Fig. 10).The resulting diagram shows two non-intersecting curves of HB (Fig. 11A).These curves define two regions of space, one where the CRAC-mediated oscillations are observed (shaded green) and These two regions can be made to intersect by decreasing the parameter K e from 800 μM to 400 μM (Fig. 12A).This change causes both the green and red regions to grow in parameter space and to intersect at δ ≈ 0.6.Physiologically, lowering K e reduces the activation threshold for the CRAC channel, with in turn causes the resting ER Ca 2+ load to decrease.At δ = 2, the V PLC bifurcation diagram exhibits the same features as the analogous diagram for K e = 800 μM.
At the other end of the diagram (δ = 0.01) the structure of the 1-parameter V PLC bifurcation diagram differs from that for the case K e = 800 μM.It is worth stressing that at small values of δ the behaviour of the open-cell model can be approximated by that for the closed-cell model for an appropiate value of C t .The K e = 800 case shows a V PLC diagram (Fig. 10C) qualitatively similar to the one generated by the closed-cell model with a value of C t = 140 (Fig. 8A).Meanwhile, the V PLC diagram for the K e = 400 μM case is qualitatively closer to the closed-cell bifurcation diagram corresponding to C t = 75 (Fig. 8E), with both of them containing branches of wide and narrow spikes.A detailed investigation of any interaction between the red and green regions as δ, among other system parameters, varies is left for future work.

Discussion
Literature on Ca 2+ dynamics inside lymphocytes is extensive, and experimental insight about the particular type of Ca 2+ oscillations observed in Fig. 1 dates back forty years and has highlighted the importance of SOCE through the CRAC channel for the oscillations to persist with time (Zweifach and Lewis 1995;Dolmetsch and Lewis 1994).Benson et al. have provided experimental evidence to corroborate this hypothesis, as well as evidence for a second type of oscillation that emerges in the absence of SOCE@.With this in mind, we have created a mathematical model based on physiological principles that can reproduce both the CRAC-mediated and IP 3 R-mediated Ca 2+ oscillations.
In cells without extracellular Ca 2+ exchange, the closed-cell model ( 2) is able to generate narrow spikes, simple wide spikes and wide spikes with an oscillatory plateau, depending on the rate of IP 3 production due to PLC activity and a degradation rate dependent on the cytoplasmic Ca 2+ concentration.Positive feedback mechanisms relying on phospholipase C (PLC) activity on c were also explored but were unable to generate oscillations with oscillatory plateaus (Lee et al. 2023).
Regarding the broad spikes with oscillatory plateaus, these solutions were observed near the right HB in the case C t = 95, and below belong to an isolated branch of periodic solutions located near the HB2 point on the diagram (Fig. 8D).While it is possible to numerically observe this branch by performing continuation on the periodic solution, we do not have a working hypothesis for the onset of this branch or the connection it might have to the main branch emanating from the HB@.We believe that the onset of this branch would require the use of slow-fast analysis on the model, a technique that is outside the scope of this paper.
The role of the parameter C t in the closed-cell model is also worth emphasizing.The model predicts that lower values of total Ca 2+ in the cell give rise to wide spike oscillations and reduces the parameter range where narrow spike oscillations are observed.This prediction could help to explain the changes in the behaviour of the response seen in the Ca 2+ -free experiments (Fig. 1); with the absence of a refilling mechanism, cells will empty their Ca 2+ stores and this is accompanied by a transition in the Ca 2+ response, from a narrow spike to a wide spike with (or without) an oscillatory plateau.
In cells with extracellular Ca 2+ exchange, the open-cell model ( 1) is capable of reproducing the CRAC-mediated oscillations observed experimentally.The time series of the open probability of the IP 3 R (Fig. 6) during such an oscillation highlights the gradual IP 3 R activity happening in the background of the strong SOCE signal, where fast-CICR has been suppressed.This finding is supported by prior experiments that were able to induce CRAC-mediated oscillations in cells by the blockage of SERCA pumps with Tg, causing a constant, but small, Ca 2+ leak from the ER (Dolmetsch and Lewis 1994).The time series presented in Fig. 6 suggests that, under TCR stimulation, the IP 3 R of these cells act in a similar fashion, causing a slow but constant leak of Ca 2+ just strong enough to induce the formation of the CRAC channel and induce SOCE@.
Our model introduces a new variable s which represents the Ca 2+ influx through the CRAC channel.The activation curve for the CRAC channel in our model needs to be rather steep for these oscillations to occur.The formation of this channel also needs to be done over a sufficiently long period to allow the ER to refill and deplete enough so that the oscillations can regenerate.This time delay on the formation of the channel can be attributed to the time needed to activate the STIM isoforms and for these to move away from the ER membrane and closer to the Orai isoforms on the PM (Trebak and Kinet 2019).We have shown that the steepness and the time delay are both necessary in our model to exhibit the CRAC-mediated oscillations.Our model currently does not make a distinction between STIM1-mediated and STIM2mediated SOCE@.Experimental evidence suggests that the two STIM molecules have a similar level of importance in this cell type, with STIM1 being slightly more important (Benson et al. 2023).Benson et al. have also performed experiments in cells without STIM1 or STIM2, and found that no CRAC-mediated response was obtained in cases where at least one STIM protein was missing.Our J SOCE activation curve could be modified to account for different types of STIM and investigate the lack of a CRAC-mediated oscillation in those cells.
The CRAC-mediated oscillations generated by the open-cell model (1) also replicate the behaviour initially hypothesized by Dolmetsch and Lewis (Dolmetsch and Lewis 1994), who predicted that the ER Ca 2+ concentration and Ca 2+ influx oscillate out of phase with each other, as seen in Fig. 4. Exploration of the corresponding time series of the IP 3 R and SERCA fluxes suggest that Ca 2+ transport through the ER membrane also oscillates with the same period, overriding the fast CICR from the IP 3 R.We conjecture then that the appearance of CRAC-mediated oscillations is dependent on the timescales of SOCE activation (τ s ) and activation of the IP 3 R (τ max ) being on the same order of magnitude.Experimentally, this can be seen in the fact that increasing either τ s from 15 to 150 or τ max from 7.5 to 75 both show no oscillations for any V PLC value, but increasing both at the same time does produce a CRAC-mediated oscillation (not shown here).This conjecture is not explored further in this paper and is left for future work.
The presence of the CRAC-mediated oscillation is also dependent on δ, as seen in Figs.11 and 12, which can be interpreted as a dependence on the PM fluxes being 'fast enough' in comparison to the ER fluxes.A proper discussion on the relative speed of the PM fluxes, however, necessitates a proper non-dimensionalization of the model, which is left for future work.
Finally, we have shown that the oscillations observed in the presence and absence of extracellular Ca 2+ have fundamentally different dynamical properties.One of the primary distinctions between the two families of oscillations is the amount of Ca 2+ transported through the ER during oscillation; in narrow spikes the ER stores are only slightly depleted, with a total loss in Ca 2+ concentration no greater that 5 μM@; in CRAC-mediated oscillations, the ER stores lose up to 30 μM of Ca 2+ .This difference in the magnitude of Ca 2+ loss from the ER is not visible from the cytoplasmic Ca 2+ concentration, which has the same order of magnitude in both narrow spikes and CRAC-mediated oscillations.
We showed that the IP 3 R-based narrow spikes and the CRAC-mediated oscillations belong to different families of stable periodic orbits that can both be found in the opencell (1) under the correct parameter conditions.We have found it possible to intersect the two curves of HB seen in Fig. 11 by varying the parameter K e .In this intersecting region, two branches of periodic orbits can coexist, which further emphasized their differing nature.We have not been able to 'pull down' the HB curve associated to the CRAC-mediated oscillations towards the boundary δ = 0, however, which emphasizes the dependence that this oscillation type has on the relative speed of the PM fluxes.It is worth remarking that these two families of oscillations can be found in the open-cell model simply by varying the parameter δ, with CRAC-mediated oscillations present at 'high' values of the parameter, and narrow spike oscillations present at 'low' values (Fig. 11).This switch in the type of oscillation implies that changes in δ are able to modulate the rate of Ca 2+ exchange through the ER membrane during an oscillation.The rate of transport through the PM is known to affect the frequency of oscillations in Ca 2+ models (Dupont et al. 2016;Sneyd et al. 2004).
There is no intrinsic reason why the two oscillation types could not 'coexist' on the same periodic solution.In fact, the open-cell model can exhibit a hybrid oscillation where a narrow spike appears on top of a CRAC-mediated oscillation using the parameter values K e = 400 μM, δ = 0.62 and V PLC = 0.1 μM alongside the values shown in Table 1.We did not include this type of orbit in the paper as it did not correspond to any of the experimental traces observed.The hypothesis that this narrow spike occurs as one of our variables crosses an activation threshold (the main candidates being c and h) is something have explored, but confirming that hypothesis would require the use of mathematical techniques that exceed the scope of this work.

Fig. 1
Fig. 1 Experimental observations in Jurkat T-cells.Representative traces of Fura-2 imaging of Ca 2+ oscillations caused by TCR stimulation of Jurkat T-cells with 125 ng/mL α-CD3.Cells in column A are stimulated in the presence of 2 mM extracellular Ca 2+ ; cells in columns B, C and D are stimulated in a Ca 2+ -free environment

Fig. 2
Fig. 2 Schematic representation of the model.The bottom panel shows a schematic representation of the open-cell model (1).Top left panel shows a closeup of the model at the ER membrane.The red arrow represents the J IP 3 R flux, transporting Ca 2+ from the ER into the cytoplasm; Ca 2+ is transported back into the ER via the J SERCA flux.Increases in cytoplasmic Ca 2+ lead to faster degradation of IP 3 (J deg ), as well as a decrease in the rate of inactivation of the IP 3 R, given by h (h ∞ ).Top right panel shows a similar closeup of the model at the PM@.The green arrow represents Ca 2+ influx into the cytoplasm given by the variable s, and Ca 2+ removal from the cytoplasm is done via the J PM flux.A depleted ER (represented here with a faded color) sends a signal to the PM to initiate influx and refill the ER (J SOCE ) (Color figure online)

Fig. 3
Fig. 3 Types of Ca 2+ oscillations produced by the model.Cytoplasmic Ca 2+ traces obtained by simulation of the open-cell (1) (Panel A) and closed-cell models (2) (Panels B, C, D) using the values in Table 1.All simulations use initial values corresponding to the equilibrium solution for V PLC = 0. Panel A uses parameter values δ = 2, V PLC = 0.1; Panel B uses δ = 0, V PLC = 0.1, C t = 140; Panel C uses δ = 0, V PLC = 0.1, C t = 75; Panel D uses δ = 0, V PLC = 0.195, C t = 95.Panels A through D use different horizontal axes to increase visual clarity.Parameter values not specified here can be found in Table 1 ) is capable of producing the following types of oscillations:1.CRAC-mediated: Parameter values δ = 2, V PLC = 0.1 μM/s.This choice simulates an open-cell environment and generates an attractor with an amplitude of c ≈ 0.2 μM and oscillatory period of T ≈ 40 seconds (Fig.3A).2. Narrow spike: Parameter values δ = 0, V PLC = 0.1 and C t = 140 μM.This choice simulates a closed-cell environment with an ER Ca 2+ concentration c e = 770 μM.The attractor has an amplitude of c ≈ 0.2 μM and period T ≈ 4 seconds (Fig.3B).3. Broad spike: Parameter values δ = 0, V PLC = 0.1 μM/s and C t = 75 μM.This choice simulates a closed-cell environment with an ER Ca 2+ concentration c e = 412 μM.The attractor has an amplitude of c ≈ 0.4 μM and period T ≈ 10 seconds (Fig.3C).4. Broad spike with an oscillatory plateau: Parameter values δ = 0, V PLC = 0.195 μM/s and C t = 95 μM.This choice simulates a closed-cell environment with an ER Ca 2+ concentration c e = 522 μM.The attractor has an amplitude of c ≈ 0.4 μM and period T ≈ 20 seconds (Fig.3D).

Fig. 4
Fig. 4 Open-cell bifurcation diagram.Panel A shows the bifurcation diagram for the open-cell model (1) using V PLC as a bifurcation parameter; the black curve represents the equilibrium solution, while the green solid lines represent the maximum and minimum c values of a branch of stable periodic orbits.The HB labels mark Hopf bifurcations of the equilibrium point.Panel B shows the s and c e traces in solid and dashed green lines, respectively, obtained from simulation of the open-cell model (1) with V PLC = 0.1 and δ = 2, as well as the parameters fromTable 1 (Color figure online)

Fig. 5
Fig. 5 Open-cell bifurcation diagrams for steepness and delay parameters.Panel A shows the s 1 -bifurcation diagram for the open-cell model (1) for V PLC = 0.1 and δ = 2. Panel C shows a similar diagram using τ s as a bifurcation parameter.In both panels, the black curve represents the equilibrium solution, while the green solid lines show the maximum and minimum c values of a branch of stable periodic orbits.Panels B, D and F show the locations in the s 1 − V PLC , τ s − V PLC and s 1 − τ s (respectively) parameter planes of the Hopf bifurcations (black curves) and regions of CRAC-mediated oscillations (green-shaded regions).Panel E shows the SOCE activation curve with a steepness value s 1 = 0.2 in green, and a reverse Hill curve with a corresponding coefficient n = 100.Parameter values not specified here can be found in table 1 (Color figure online)

Fig. 6
Fig. 6 Modulation of IP 3 R activity by CRAC.Time series of the open probability of the IP 3 Rs during a CRAC-mediated oscillation (δ = 2) in green, and a similar series during a narrow spike (δ = 0) in red.In both cases the time series were obtained by stimulating the open-cell model (1) with V PLC = 0.1.Parameter values not specified here can be found in table 1 (Color figure online)

Fig. 7
Fig. 7 Initial high loads of ER Ca 2+ lead to transient narrow spikes.Panels show stimulation of the opencell model (1) (represented by an increase in V PLC from 0 to 0.1) under different initial c e values.The corresponding initial c and s values are obtained by solving equations 1b and 1e assuming the cell is unstimulated (V PLC = 0).The top panels use initial values c = 0.809, c e = 809, s = 0.421, while the bottom panels use c = 0.850, c e = 850, s = 0. Panels A and C show the time series for the cytoplasmic Ca 2+ concentration (c), panels B and D the J SOCE curve (4) as a function of c e alongside a projection of the corresponding trajectory on the s − c e plane.The initial conditions are marked with a red dot for visibility.All panels use a value of δ = 2. Parameter values not specified here can be found in table 1 (Color figure online)

Fig. 8
Fig. 8 Closed-cell model bifurcation diagram.Panel A shows the complete bifurcation diagram of the closed-cell model (2) using V PLC as a bifurcation parameter, for C t = 140.Panels C and E show similar (partial) bifurcation diagrams for C t = 95 and C t = 75, respectively.The black curve represents the equilibrium solution, while the red lines show the maximum and minimum values of c on the branches of periodic orbits; stable orbits are represented by solid lines and unstable orbits by dashed lines.Panels B, D and F show the time series for the cytoplasmic Ca 2+ concentration in solid lines, and the cytoplasmic IP 3 concentration in dashed lines.Panel B shows the attractor obtained for parameter values V PLC = 0.1, C t = 140.Panel D shows two types of attractors obtained for C t = 95, for V PLC = 0.09 (top) and V PLC = 0.195 (bottom).Panel F shows two types of attractors obtained for C t = 75, for V PLC = 0.06 (top) and V PLC = 0.1 (bottom).Parameter values not specified here can be found in table 1 (Color figure online)

Fig. 9
Fig. 9 Two-parameter bifurcation set for the closed-cell model.The black line represents a curve of Hopf bifurcations (HB); the red line a curve of period doubling bifurcations (PD); the orange line represents a curve of saddle nodes of periodic orbits (SNPO) emanating from the right branch of the HB curve; the dotted green line represents a curve of Neimark-Sacker bifurcations (TR).Labels I, II, III and IV indicate the approximate region in parameter space where narrow spikes, wide spikes, wide spikes with plateau, and no oscillations are observed, respectively.Parameter values not specified be found in table 1 (Color figure online)

Fig. 11
Fig. 11 Two-parameter bifurcation set for the open-cell model.Panel A shows the two-parameter bifurcation set for varying δ and V PLC .The border of the upper green region represents a curve of HB observed for high values of δ, while the green region itself represents the parameter region where the CRAC-mediated attractor is observed.The border of the red region is a curve of HB observed at δ values close to zero.Panels B, C and D show the bifurcation diagrams of the open-cell model (1) using V PLC as a bifurcation parameter, for fixed values δ = 2, 0.62, 0.01, respectively.The black curve represents the equilibrium solution, while the red and green lines represent branches of periodic orbits.Stable attractors are represented by solid lines and unstable attractors by dashed lines.Parameter values not specified here can be found in table 1 (Color figure online)

Table 1 (
Color figure online)