Mathematical Analysis of Membrane Transporters Dynamics: A Calcium Fluxes Case Study

A tight control of intracellular [Ca2+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{2+}$$\end{document}] is essential for the survival and normal function of cells. In this study we investigate key mechanistic steps by which calcium is regulated and calcium oscillations could occur using in silico modeling of membrane transporters. To do so we give a deterministic description of intracellular Ca2+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{2+}$$\end{document} dynamics using nonlinear dynamics in order to understand Ca2+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{2+}$$\end{document} signaling. We first present the ordinary differential equations (ODEs) system for cell calcium kinetics and make a preliminary work on Sobol indices. We then describe and analyze complex transporters action. Besides, we analyze the whole system. We finally perform numerical simulations and compare our results to real data.


Introduction
The ability of a cell to perceive and correctly respond to its microenvironment depends on very complex signaling systems. This includes various membrane and transmembrane proteins, which are at the interface between the extracellular and intracellular media. Because of the wide repertoire of spatio-temporal fluctuations in its intracellular concentrations Berridge et al. 2000 ) signaling is exquisitely poised to link extracellular mechanisms with intracellular modifications and to control cell fate and many cell functions including regulation of metabolism, proliferation, death, gene transcription, cell migration, exocytosis, and contraction (Berridge 2016). The pattern of intracellular calcium signaling such as intracellular calcium oscillations that can be repeated over a longer period (Berridge 1990;Fewtrell 1993) is now considered as a universal mechanism of signal transduction and determines specific cell states; it is also involved in specific stem cells activity (Terrié et al. 2019) or in various cancer cells hallmarks (Déliot and Constantin 1848;Prevarskaya et al. 2018).
Since a tight control of intracellular [Ca 2+ ] is essential for the survival and normal function of cells, resting calcium activity is maintained at a low level (between 50 and 200 nM (Behringer et al. 2017)) in order to keep a large dynamics range for the calcium signal. Membrane ionic channels and transporters, cytosolic calcium buffers and calcium buffering organelles regulate calcium influx, storage and extrusion to maintain [Ca 2+ ] below the activation thresholds and extraphysiological values.
The large concentration gradient of Ca 2+ is created by specialized proteins named Ca 2+ -ATPases, which are considered in this study, namely, (i) Plasma membrane Ca 2+ -ATPase (PMCA) located on the plasma membrane, which extrude Ca 2+ to the extracellular space or (ii) Sarco/endoplasmic reticulum Ca 2+ -ATPase (SERCA) inserted in the membrane of the endoplasmic reticulum (ER) that is an important intracellular Ca 2+ pool. Calcium signal diversity and the generation of receptor-triggered Ca 2+ signals also rely on cell proteins specialized in calcium transport through cell membranes, which ensures the transitory elevation of free cytosolic calcium concentration through Ca 2+ flux into the cytoplasm from different compartments. We consider two main Ca 2+ fluxes, which could contribute to responses to extracellular signals and/or intracellular calcium oscillations: (i) a calcium entry along its concentration gradient from the extracellular compartment across the plasma membrane into the cytoplasm of the cell and (ii) a calcium release into the cytoplasm of the cell from intracellular stocks mainly contained in the ER. In the model we propose, both Ca 2+ fluxes (influx and release) depend on the activity of phospholipase C (PLC), which generates the intracellular second messengers inositol 1,4,5-trisphosphate (IP3) and diacylglycerol (DAG). We also consider in our model the intracellular fluctuations of these two essential second messengers triggering both Ca 2+ fluxes into the cytoplasm of the cells: (i) IP3 diffuses rapidly within the cytoplasm and activates the inositol 1,4,5-trisphosphate receptor (IP3R), which releases Ca 2+ ions from the ER (Mikoshiba 2015;Berridge 2016), and usually lead to rapid and transient increase of cytosolic Ca 2+ ; (ii) DAG, the receptor-operated calcium channel (ROCC), activates plasma membrane which supports entry of Ca 2+ and participates to slow and sustained increase of cytosolic Ca 2+ Vazquez et al. 2001;Dietrich et al. 2005. We also focus on the participation of a second crucial pathway of Ca 2+ influx, the store-operated calcium channel (SOCC), which is highly dependent of the intracellular Ca 2+ pool of ER (Parekh and Putney 2005) and activates after Ca 2+ release.
The intracellular calcium dynamics in various non-excitable cells of the brain such as astrocytes, glial cells, neural stem cells and also glioblastoma stem cells involves the interplay of all these calcium fluxes from and into the ER and across the plasma membrane. This interplay conditions the pattern of the intracellular calcium signal such as a calcium transient with various durations in response to extracellular factors activating receptors or a calcium oscillation for instance. Numerous computational models were previously presented to describe intracellular Ca 2+ oscillations (see the reviews Dupont et al. 2011 andDupont 2014). Other studies assume that Ca 2+ dynamics remains a stochastic process even at the cellular level and that there is poor communication between Ca 2+ channels, for instance between IP3R Ca 2+ -release channels. Most of the computational models for intracellular Ca 2+ dynamics focus on the IP3R Ca 2+ -release channels which play a central role in all the cells (Dupont 2014). Only some calcium studies give a mathematical analysis (including well-posedness) of their models.
In this study we consider a positive feedback of Ca 2+ on IP3R, due to the mechanism of Ca 2+ -induced Ca 2+ release (CICR), and we use the bell-shaped dependency of IP3R activity on cytoplasmic IP3 second messenger (see for example Houart et al. 1999, Lavrentovich andHemkin 2008). Some of these models assume that the IP3 concentration is at a constant stimulatory level, as in Borghans et al. (1997). It has become possible to monitor IP3 concentration variations in intact cells. These experiments have shown that the IP3 concentration can oscillate together with cytoplasmic calcium (Hirose et al. 1999;Thore et al. 2004). Based on other works (see Politi et al. 2016) we assume that concentrations of both IP3 and DAG produced by PLC activity are not constant but oscillating in the cytoplasm. The existence of both positive and negative feedbacks of Ca 2+ on IP3 metabolism could mediate fluctuations in cellular IP3 levels. We also consider a positive and negative feedbacks of Ca 2+ on IP3 metabolism: Ca 2+ activation of PLC producing IP3 and Ca 2+ activation of IP3 3-kinase degrading IP3, respectively. Since the steep calcium gradient across the ER membrane is sustained by active pumping by SERCA, most studies considering the release of Ca 2+ through IP3R calcium release channels also take into account this pumping activity that can be modeled by a Hill expression (Dupont 2014).
Some newly identified players in intracellular Ca 2+ dynamics have been less often incorporated in computational models, mostly because their kinetics are poorly described quantitatively. As already mentioned, release of stored Ca 2+ from the ER activates influx through SOCCs, which depend on activation of STIM1 oligomers, sensing the decrease of Ca 2+ in the ER lumina. SOCCE thus play an important part in the interplay between Ca 2+ of the ER and Ca 2+ fluxes from the external pool through the plasma membrane. Moreover, increasing evidences show that changes in SOCCE or ROCCE may affect in some cells, which are more dependent of external Ca 2+ , the amplitude and duration of Ca 2+ signals and frequency of Ca 2+ oscillations. The study of Liu et al. (2010) proposed a model that reproduces the steady-state dependence of the SOCCE on the level of Ca 2+ in the ER lumina, assuming a cooperative binding of ER Ca 2+ to STIM. Another model, introducing a delay between store emptying and Orai opening (Croisier et al. 2013), has been proposed. To our knowledge no mathematical analysis takes into account at the same time the calcium 14 Page 4 of 32 release and the pumping flux with Ca 2+ entries through SOCCE and ROCCE, both depending on PLC activity, yet.
In this study we investigate key mechanistic steps by which calcium oscillations could occur using in silico modeling. To do so we present here a deterministic description of intracellular Ca 2+ dynamics using nonlinear dynamics in order to understand Ca 2+ signaling. We first present the ordinary differential equations (ODEs) system for cell calcium kinetics and make a preliminary work on Sobol indices. We then describe and analyze complex transporters action. Besides, we analyze the whole system. We finally make simulations and compare our results to real data.
The mathematical analysis focuses first on the action of several transporters separately and study their combined effect on calcium concentrations. We have chosen to study the combined effect of the IP3-dependent flux leading to release of ER Ca 2+ stores, with two different pathways of Ca 2+ entry, ROCCE (depending on DAG) and SOCCE (depending on Ca 2+ ER stores), together with transporters supporting Ca 2+ uptake into the ER and Ca 2+ extrusion through the plasma membrane by two different transporters (PMCA and NCX). We also introduce a continuous Ca 2+ flux from ER to cytoplasma (leak channels) in addition to the IP3-dependent Calcium release, since this leak can be observed in all cell when the pumping activity of SERCA is inhibited. The IP3-dependent calcium release, ROCCE and SOCCE, all depend on the activity of PLC. This enzyme is also under the control of membrane receptors responding to various extracellular signals, whose patterns of change over time are unknown if not controlled experimentally. These external signals stimulating the intracellular activity of PLC are set to constants in some mathematical studies cited previously, but could also be defined as oscillatory signals. Modeling allows us to explore the combined effects of membrane Ca 2+ transporters in presence of extracellular stimulation or in absence, considering only the mechanistic interplay between the membrane transporters.

Mathematical Modeling
The present model is proposed in order to study spontaneous calcium oscillations in glial cells or neural stem cells based on in vitro observations and literature results. It is built in vitro. Therefore, it is considered closed with exchanges with the extracellular calcium. Two compartments are considered: the cytosol and the endoplasmic reticulum (ER) as defined in Fig. 1. We also introduce four metabolic concentrations: calcium on the cytosol (C), calcium on the ER (R), inositol trisphosphate or IP3 on the cytosol (I) and diacylglycerols or DAG on the cytosol (D) which are produced by the Phospholipase C (PLC) activity. The external calcium concentration E is highly regulated; we therefore assume that it is constant through time. We do not follow its dynamics.
Cytosolic calcium concentration is subject to two main phenomena: exchanges with the ER and exchanges with the extracellular space. First we focus on calcium exchanges between the cytoplasm and the ER. Three different calcium transporters are included in the model. First an ATP-dependent calcium uptake from the cytosol to the ER done by the SERCA pump at the maximal rate V s and with the half saturation k s . Then the IP3 receptor channel (IP3R) with the maximal transport rate V p depending on the reticulum calcium concentration (with the half saturation k r,p ) and the IP3 concentration (with the half saturation k i,p ). It also varies according to the cytosolic calcium concentration with a maximal affinity of p and a standard deviation of p . Finally there is an additional calcium leak flux at the maximal rate V f . Figure 1 shows the different interactions.
We have the four different kinds of calcium exchanges between the cytosol and the extracellular space. Among them there is the sodium-calcium exchanger NCX with the maximal rate V n and with the half saturation k n for calcium exit. There is also the PMCA pump of maximal rate V m and of half saturation k m allowing calcium outflux. The ROCC channel of maximal rate V T depends on external stimulations f, DAG concentration (with half saturation k D,T ) and on a balance between the cytosolic calcium the concentration and the external calcium concentration. To do so we define the quantity where E is the external calcium concentration and k E,T is the constant of half saturation for external calcium through the ROCC channel. We also denote by k C,T the half saturation for cytosolic calcium through the same channel. Finally the protein STIM1 is able to activate the SOCC in the plasma membrane. When the reticulum endoplasmic calcium concentration is lower than R o , the SOCC channel opens, stimulated by the protein STIM1, allowing a calcium influx of maximal rate V o from the extracellular space to the cytosol.
The molecule IP3 is a product of the hydrolysis of PIP2 by the enzyme PLC that, through the same reaction and in the same proportions, produces DAG. This creation has a maximal rate V L and depends on both external sinusoidal stimulations f and the cytosolic calcium concentration with the half saturation k L . These molecules are then degraded in both a linear way and on a way linearly depending on the cytosolic calcium concentration. We denote by i the proportion of IP3 linearly degraded, V i the maximal rate of degradation related to the cytosolic calcium concentration and k i the constant of half saturation. Degradation of DAG is supposed to follow the same pattern with different rates (results not shown). Therefore we denote by d the proportion of DAG linearly degraded, V d the maximal rate of degradation related to the cytosolic calcium concentration and k d the constant of half saturation.
Finally, we have the following ordinary differential equations (ODEs), ∀ t ∈ ℝ + , The first equation stands for cytosolic calcium concentration kinetics. The three terms on the first line describe exchanges between the cytosol and the ER. The five terms on the second and third lines stand for exchanges between the cytosol and the extracellular space. Accordingly, equation 2 describes ER calcium concentration kinetics. Equation 3 relates IP3 influx and outflux while equation 4 describes DAG dynamics. All the parameters meaning are detailed above. The initial condition is given by: Assumption 1 If not precised the involved parameters are assumed to be positive.

Assumption 2
The function f is supposed to be periodic, nonnegative and bounded. Therefore ∃F ∈ ℝ * such that ∀ t ∈ ℝ + : (1) In cells the IP3-dependent calcium release, ROCCE and SOCCE all depend on the activity of PLC. This enzyme is also under the control of membrane receptors responding to various extracellular signals that is described in our model by using the function f. This function f representing the external signals stimulating the intracellular activity of PLC is set constant in our study, or as oscillatory signals.

Preliminary Work: Sensitivity Analysis
We perform a global sensitivity analysis to assess the impact of model parameters on the main model output. In this study, we analyze the effect of each parameter on the cytosolic calcium concentration for each time through first order Sobol indices (Iooss and Lemaître 2015) given by the following expression : where C is the cytosolic calcium concentration for a given time, X i is a model parameter, V(C) represents the total variance of C, C[Y|X i ] is the conditional mean of C given X i , and N is the number of model parameters. Therefore, the value S i measures the part of Y variance that is explained by parameter X i . In other terms, first order Sobol indices determine how much the cytosolic calcium concentration varies when a chosen parameter value varies. Sobol indices are always between 0 and 1. The closest to 1 a Sobol indice is, the higher the model output is sensitive to the related parameter variations. Using Sobol analysis helps to determine parameters contribution to the output throughout time. We limit our study to first order Sobol indices. The impact of the interaction of several parameters on model output can be assessed through other Sobol indices. However they are more difficult to interpret. The sum of all Sobol indices is equal to 1. Each sub-figure in Fig. 2 displays Sobol indices variations through time for a given parameter. We display the higher Sobol indices. First order indices for other parameters are lower than 10 −3 . The first order Sobol indices sum is approximately 30%.
This preliminary work underlines the importance of maximal transporter rates and in particular those of the SERCA ( V s ), SOCC ( V o ), NCX ( V n ) and PMCA ( V m ). In this study we will therefore focus on these transporters impact. Moreover the impact of R o has to be put into perspective because it is only high on the early stages. Besides, sums of first order Sobol indices in long times are close to 0.3 (id est not close to 1). This means that most of the time, parameters interact with each other to impact the variability of the cytosolic calcium concentration. We will investigate the nature of these interactions in Sect. 5.

Preliminary Work on Transporters
Different kinds of transporters are involved in the studied dynamics. We therefore make a step by step analysis by focusing first on each transporter dynamics separately. In this section, we analyze three kinds of possible reactions from a transporter to a substrate, • Symport like transports. Symport like transporters work in the membrane and several molecules are transported across the cell membrane at the same time, and is, therefore, a type of cotransporter. It possess a saturation value. The symport like transport is defined by a maximal rate V, a half saturation value k and a degree of reaction i, i ∈ {1, 2} , such as ∀ x ∈ ℝ , • The bell-shaped transporter. It describes a transporter possessing a maximal affinity with a substrate for a given value of the latter. The bell-shaped transport involved in IP3R dynamics has a maximal affinity p and a standard deviation p with cytosolic calcium concentration such as ∀ x ∈ ℝ ,

Lemmas on Symport Like Transporters
We Lemma 1 For the symport transporter of degree 1 (id est i = 1 ) and for x ∈ ℝ + , then T i is V k -Lipschitz continuous, nonnegative and bounded by V.
Proof We set, ∀x ∈ ℝ + , We therefore have ∀x ∈ ℝ + , Let x and y be in ℝ + . Then, so that T 1 is V k -Lipschitz continuous. ◻ Lemma 2 For the symport transporter of degree 2 (id est i = 2 ) and for x ∈ [0, A] , A ∈ ℝ + then T i is 2VA k -Lipschitz continuous, nonnegative and bounded by V.
Proof We set, ∀x ∈ ℝ + , We therefore have ∀x ∈ ℝ + , Let x and y be in [0, A], A ∈ ℝ + . Then, so that T 2 is 2VAk-Lipschitz continuous. ◻ In Fig. 3 we give examples of symport like functions.

Study of the Bell-Shaped Transporter
We define for p ∈ ℝ + , p ∈ ℝ + ∀ x ∈ ℝ, Lemma 3 For x ∈ ℝ then H p is Lipschitz continuous, nonnegative and bounded by 1.
Proof We set ∀ x ∈ ℝ, Then its derivative is ∀ x ∈ ℝ, and the second derivative is ∀ x ∈ ℝ, We therefore have the following variation table.
The derivative H ′ p is bounded by 1 p exp(− 1 2 ) and therefore H p is Lipschitz continuous. ◻ In Fig. 4 we give an example of a bell-shaped H p function.

Study of the Switch Transporter
We define for .
Then its derivative is ∀ x ∈ ℝ, and we have ∀ x ∈ ℝ, Therefore G 0 is V o -Lipschitz continuous and we have the following variation table.

Lemma 5 The product of Lipschitz continuous and bounded functions is Lipschitz continuous.
Proof We assume that the function f is k f -Lipschitz continuous such that ∃F ∈ ℝ + ∕ |f | ⩽ F . We also assume that the function g k g -Lipschitz such that ∃G ∈ ℝ + ∕ |g| ⩽ G. Let x, y ∈ ℝ . We have

Mathematical Analysis
In force of Sect. 4, we will prove the existence and nonnegativity of the solutions. We will also exhibit bounds on the solutions and a steady-state.

Upper Bounds
From 1-4, we have the following inequalities : Inequalities 7 and 8 imply, using Gronwall's lemma, and Lemma 6 We can exhibit a sufficient, but not necessary, condition to ensure a bound on C and R. Let V n , V m , F, V T , E T and V o be such that : In that case, ∃B C , B R ∈ ℝ + * Proof Setting ∀ x ∈ ℝ + and using 1-4 we have We want to study u variations for x ∈ ℝ + . Its derivative reads so that u is a strictly decreasing function and Therefore, assuming that V n + V m ⩾ FV T E T + V o , we have, using the intermediate value theorem, ∃!y > 0 such that u(y) = 0 . Moreover for all x ∈ ℝ + , x > y , we have u(x) < u(y) = 0 . If this condition is not satisfied, u is nonnegative. Therefore using Eq. 12 and assuming V n + V m ⩾ FV T E T + V o , there exists an unique nonnegative constant N c such that u(N c ) = 0 and for C > N c u is strictly negative meaning that R � + C � is strictly negative.
Let C be such that C > N c . Then R � + C � is negative and at least C ′ or R ′ has to be negative. Let us show that, assuming R ′ is negative, C ′ has to be negative too. Using 2, we have ∀ t ∈ ℝ + so that, using 1 Therefore C � (t) is negative, so that R � + C � < 0 implies C ′ < 0.
We will now demonstrate that under this condition R et C are bounded. For C > N c , we proved that C ′ < 0 , so that Using 5, this implies that so that, using Gronwall lemma, and therefore R and C are bounded. ◻

Remark 1 Lemma 6 ensures bounds on C and R under the condition
This result has to be compared with Sobol indices results in Sect. 3. It is assumed that V n , V m , V o and V T have a strong impact on the dynamics in both studies.

Existence of the Solution
Since we have nonnegativity and bounds of C, R, D and I then, based on Sect. 4, Lemma 1 to Lemma 5, we can ensure that terms related to symport like, bell-shaped or swtich transporters are Lipschitz continuous. Moreover products and sums of Lipschitz continuous and bounded functions are Lipschitz continuous. Then we can rewrite (1)-(4), setting to have ∀t ∈ ℝ + : where H is locally Lipschitz continuous with respect to the second variable. We finally conclude, thanks to the Cauchy-Lipschitz theorem, that we have existence and uniqueness of the solution to the system ∀t ∈ ℝ + .

Steady-State Study
Assuming f = 0 , the steady-state of system 1-4 is given by We define the function S such that, ∀x ∈ ℝ 2 ,

X(t) ∶= (C(t);R(t);I(t);D(t)),
X � (t) = H(t, X(t)), X(0) = X 0 , We are looking for a solution at E(C) = 0 . The related derivative reads We therefore have Using the intermediate value theorem, ∃ ! x c > 0 such that S(x c ) = 0 . The related steady-state is therefore given by We deduce that there is only one steady-state. Moreover at equilibrium, the ER calcium concentration is higher than the cytosolic concentration. This difference is directly linked to the maximal capacity of the SERCA pump ( V s ) and the one of the leak channel ( V f ). This underlines the importance of these transporters capacity.
We will show that the stead-state is always included on the viability domain defined by the bounds. We know that C and R are nonnegative. We have to verify they are lower than the related upper bounds. In other words we want to show that x c ⩽ B c = max(C(0), N c ) . To do so it is sufficient to prove that x c ⩽ N c . We also want to prove that R < B r Remark 2 This verification is only possible when upper bounds are exhibited, that is to say when V n + V m ⩾ FV T E T + V o . If this inequality does not hold there is no more analysis to do and the steady-state is in the viability domain because concentrations are greater than 0. Having a steady-state does not mean that variables are bounded; in fact the steady-state can be unstable.

Hence, recalling that S(x c ) = 0,
Recalling that u is strictly decreasing and that u(N c ) = 0 , then u(x c ) > u(N c ) and we have Hence C is always on the viability domain of C.
We will now prove that R < B r . We know that so that we have and R is always on the viability domain of R.

Simulations
In this section, we first comment on data from the relevant literature and from our in vitro experiences. We then present several numerical simulations, assuming a constant function f = F . We also compare the numerical simulations with different values of V n , V m , V o , V T and V s . Then we present simulations with a periodic function f. We finally compare our results to experimental data. These simulations have been done with the Matlab software. We used the functions ode45 and ode23s to simulate the whole dynamics. These functions adjust integration time steps. Because of the important number of the involved parameters, we will not give exact parameters values for each simulation in the article. They are available upon request.

Literature and Experimental Data
Spontaneous oscillations were measured in human glioma stem cells. Cells were incubated with a permeable ratiometric calcium probe Fura2. This probe binds free cytosolic Ca 2+ and allows the measurement of cytoplasmic calcium variations in a solution containing 1.8 mM of calcium. The intensity of fluorescence of each cell were measured and next converted into calcium concentration by calibration. Figure 6 shows 3 representative measurements of spontaneous oscillations of cytoplasmatic calcium in 3 different resting human glioma stem cells. Only few cells in different recording area produced oscillation in basal conditions. The figure shows the oscillation recorded during 60 to 90 s. The amplitude and the frequency of these oscillations are different in the 3 cells. For example, the amplitude of the first graph is between 0,6 and 0,3 M and around 0,01 M in the 2 last graph.
Other calcium oscillation cell shapes can be found in the work of Borghans et al. (1997) and Keener and Sneyd (2009). They display different frequencies, amortizations and period shapes. Therefore it is well-known that cytosolic calcium may present oscillations. What is questionable is the origin of these oscillations (Parekh 2011). There are two major hypotheses: • External stimulus variations lead to cytosolic calcium oscillations through pump fluxes. While the external calcium concentration can be directly impacted by protein variations, other variations may affect pumps and channels behavior and therefore change the calcium dynamics. Fig. 6 Representative cytosolic calcium oscillations. Experiments were done in human glioma stem cells with the calcium sensitive probe Fura2 • Since cell calcium is self-regulated, oscillations arise from the mechanical properties of the dynamical system. In other words, if there is not enough cytosolic calcium, pumps and channels sum of efforts will make it increase. A contrario, if it is too high according to a reference value pumps and channels sum of efforts will make it decrease.
We use in silico modeling to address this question. For that purpose, we base our work on the study of (Lavrentovich and Hemkin 2008) adding the results of These values help to find accurate model dynamics.

Mechanical Periodicity
To test our model, we provide simulations setting f = F . In Fig. 7 we give cytosolic calcium, ER calcium, IP3 and DAG trajectories for a chosen set of parameters values. This set of parameters is kept as reference in this subsection All the concentrations remain nonnegative through time. On the one hand and since they have the same production rate and comparable degradation rates, IP3 and DAG concentrations are almost the same, showing only a higher degradation rate for DAG. On the other hand, cytosolic calcium and ER calcium have complementary dynamics leading to a calcium transfilling (Fig. 8). Fig. 7 Trajectories of cytosolic calcium, ER calcium, IP3 and DAG with no external stimulation (F constant). We can observe a periodic behavior for all the concentrations with a period of little more than an hour This model supports the idea that spontaneous Ca 2+ oscillations can be generated without the aid of external stimulation (Parri and Crunelli 2003). However related oscillations occur in long time and look different that recordings obtained in living cells, reported by Fig. 6.
To figure out how these mechanical periods may vary, we test several values of some maximal fluxes: V n for the NCX (Fig. 9), V m for the PMCA (Fig. 10), V T for the ROCC (Fig. 11) and V f for the leak flux (Fig. 12). Figure 13 investigates R o impact which is the value of switch for SOCC. Moreover since the SERCA is supposed to have a strong impact on the dynamics, we give results for crossed k s and V s variations on Figs. 14 and 15. Trajectories of cytosolic calcium, ER calcium, IP3 and DAG with no external stimulation (F constant) and different values of V n . Three kinds of dynamics set apart: the one for V n < 0.25 having a "high" stable state, the one for 0.25 < V n < 0.35 with a periodic behavior and the one for V n > 0.35 having a "low" stable state. Figure fig:meca1 gives reference parameters Fig. 10 Trajectories of cytosolic calcium, ER calcium, IP3 and DAG with no external stimulation (F constant) and different values of V m . Two kinds of dynamics set appart: the one for V m > 0.07 having a stable state and the one for V m < 0.07 with a periodic behavior. Figure fig:meca1 gives reference parameters Fig. 11 Trajectories of cytosolic calcium, ER calcium, IP3 and DAG with no external stimulation (F constant) and different values of V o . Again three kinds of dynamics set appart: the one for V o < 0.03 having a "low" stable state, the one for 0.03 < V o < 0.07 with a periodic behavior and the one for V o > 0.07 having a "high" stable state. Figure fig:meca1 gives reference parameters Fig. 12 Trajectories of cytosolic calcium, ER calcium, IP3 and DAG with no external stimulation (F constant) and different values of V f . The value of V f has an impact on the period observed, the higher V f is, the faster the dynamics is. Figure fig:meca1 gives reference parameters Fig. 13 Trajectories of cytosolic calcium, ER calcium, IP3 and DAG with no external stimulation (F constant) and different values of R o . Two kinds of dynamics set appart: the one for R o < 300 having a stable state and the one for R o > 300 with a periodic behavior. When it exists, the steady-state for R is given by the value of R o . Figure fig:meca1 gives reference parameters  (F constant) and different values of V s and k s . Not only is the value of V s by itself important but also the union of all given transport parameters ( V s but also k s for the SERCA). Figure fig:meca1 gives reference parameters Fig. 15 Trajectories of ER calcium with no external stimulation (F constant) and different values of V s and k s . Not only is the value of V s by itself important but also the union of all given transport parameters ( V s but also k s for the SERCA). Figure fig:meca1 gives reference parameters Results from Figs. 9,10,11,12,13,14,and 15 prove that there exists a set of values for which the mechanical system exhibits a periodic behavior. Moving a parameter value outside the related domain makes the system reach a steady-state. Moreover, when all the parameters but one are fixed, letting the last parameter vary makes the system reach at most two steady-states. However all these values do not impact the whole dynamics separately. Letting two parameters vary at the same time on 14 and 15 we prove the existence of crossed dynamics. De facto it seems impossible to explain the dynamics based on only one or two parameters. We thus prove similar results to those with Sobol indices in Sect. 3.

Adding External Proteic Variations
The goal of this subsection is to investigate in more details the conditions in which complex oscillatory phenomena occur in the model and to characterize more thoroughly the various modes of dynamical behavior that can be obtained. To test our model, we provide simulations for F > 0 , freq > 0 and ∀ t ∈ ℝ + , This function f is periodic, nonnegative and bounded by F. It describes external proteic variations. It is displayed in Fig. 16. In Fig. 17 we give cytosolic calcium, ER calcium, IP3 and DAG trajectories for a chosen set of parameters values. This set of parameters is kept as reference in this subsection. Figure 18 gives portrait phases of the related concentrations. Cytosolic calcium displays two kinds of periods. These results are in accordance with literature data (see Sect. 7). It is difficult to access ER calcium using in vivo experiences. Figure 16 describes ER calcium periods using in silico modeling. Each period is made of two phases: one growing with oscillations and one roughly decreasing. IP3 and DAG concentrations oscillate in phase. Figure 19 displays cytosolic calcium variations in time according to different values of freq.
This figure displays a particular phenomenon for a given interval of frequencies (here [0.35, 0.45]). At this frequency, two kinds of oscillations appear with different periods. For the majority of tested frequencies, the dynamics has only one kind of Fig. 17 Trajectories of cytosolic calcium, ER calcium, IP3 and DAG with external stimulation (periodic f). We can observe a periodic behavior for all the concentrations with a period of ∼ 60s = 1m  Figure fig:osc2 gives reference parameters oscillation. It is possible to test several frequencies of f to see how the cytosolic calcium concentration varies. Figure 20 displays cytosolic calcium variations in time according to different values of freq.
These Figures display particular oscillations behavior with two frequencies. The most present oscillations are huge spikes at low frequencies and small spikes at huge frequencies.
Finally Fig. 21 displays what happens when f = 0 . In that case one can observe a steady-state as predicted by the mathematical analysis.

Conclusion
In this study we analyze a model for calcium dynamics at a cellular level using nonlinear dynamics. This model is a first step in view of a better understanding of calcium dynamics between the cytosol and the ER. Indeed, even if calcium oscillations are often described and discussed (Keener and Sneyd 2009), this phenomenon is complicated by the wide diversity of the nature of implicated dynamics. To the best of our knowledge, no mathematical analysis focuses on transporters actions. In this paper, we study the action of several transporters separately. We then study their combined effect on calcium concentrations. We obtain existence, uniqueness and bounds on the solutions for the related system. We also underline the particular impact of given maximum transporter rates. We finally give several numerical simulations, for different values of the main parameters, and we compare the model with in vitro data.

Fig. 21
Trajectories of cytosolic calcium, ER calcium, IP3 and DAG with f equal to zero. Note that the time scale is not the same. IP3 and DAG concentrations quickly reach a steady-state at 0 (time scale 500 s) while calcium concentrations reach more slowly a non-zero steady-state (time scale 100,000 s ∼ 28 h). Figure 17 gives reference parameters We show that the mechanism through which cytosolic calcium oscillates is complex and involves several transporters action with different dynamics described previously. In particular, this study suggests that 1. Not just the mechanical action of transporters nor the external proteic variations are responsible for observed oscillations, but a combined action of both is. In fact, while taken apart, these actions may make the system oscillate but the resulting dynamics are not the one observed in in vivo nor in vitro experiences. Moreover it displays only one kind of period when the experimental dynamics have two periods. 2. There is a special parameters domain for which the system oscillates. This domain is defined by the combined action of involved parameters. In other words, we cannot define the precise value of one parameter to make the system oscillate focusing only on the given parameter. This result may explain why oscillations are not always observed in in vitro experiences and may have different shapes (Borghans et al. 1997;Keener and Sneyd 2009). 3. The periodic behavior is only true for a narrow set of parameters including V s and k s for the activity of the SERCA but also maximal rates of calcium extrusion through NCX exchanger and the calcium pump PMCA or calcium entry through TRPC channels. Moving a parameter value outside this window of values makes the system reach a steady-state. 4. Spontaneous periodic variations of calcium can be obtained in the cytoplasm and the ER lumina that display complementary dynamics. The resulting calcium transfilling between the two compartments maintains the necessary calcium stores in the ER lumina but also sustains a periodic variation of the cytoplasmic calcium concentration that follows the trajectories of IP3 and DAG, which are product by the basal activity of PLC without external stimulation. These mechanical periods have however nothing to do in frequency and shape with the cytoplasmic calcium oscillations that can be recorded in cells, but highlight that the mechanical properties of the system per se and the interplay between the transporters and the second messenger can support a basic calcium oscillation. 5. Some transporters stand out in both the Sobol indices and the deterministic analysis. It appears that the SERCA, and in particular its maximal transport rate, but also the leak channels are the most impactful on the definition of the steadystate. In other words, this modeling shows that the maximal transporters rates of SERCA ( V s ), Orai1 ( V o ), NCX ( V n ) and PMCA ( V m ) have a high impact on the cytosolic Ca 2+ flux and variation, for small changes of their values. It suggests that the activity of these transporters is critical for calcium homeostasis. 6. The fluxes maximal rate between the cytosol and the extracellular space (both inside and outside through the SERCA pump and the leak channel) maintain the difference in calcium between the two compartments and is directly linked with the maximal capacities V s and V f of these transporters. 7. The Store-operated channels (SOCC) seem to have a role of safeguard. The switch value permits a quick regulation of all the concentrations. Indeed with different value of maximal SOCC transport rate R o , the maximal calcium, IP3 and DAG concentrations are not the same and possibly limit oscillations (Fig. 13).
This in silico study aims at giving some clues on calcium oscillating behaviors between ER, cytosol and extracellular space. It could be completed with in vivo or in vitro experiments to support or reject results on calcium dynamics detailed previously and improve calcium modeling. Moreover, more precise quantitative data on calcium, IP3 or DAG concentration or transporters rates would significantly improve model calibration and de facto model prediction capacities.