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

Angélique Perrillat-Mercerot, Nadine Déliot, Alain Miranville, Rémy Guillevin, Bruno Constantin

► To cite this version:

Angélique Perrillat-Mercerot, Nadine Déliot, Alain Miranville, Rémy Guillevin, Bruno Constantin. Mathematical Analysis of Membrane Transporters Dynamics: A Calcium Fluxes Case Study. Acta Biotheoretica, 2022, 70 (2), pp.44. 10.1007/s10441-022-09437-3. hal-03856308

HAL Id: hal-03856308

https://hal.science/hal-03856308

Submitted on 17 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Mathematical analysis of membrane transporters dynamics, A calcium fluxes case study

A. Perrillat-Mercerot<sup>a</sup>, N. Deliot<sup>b</sup>, A. Miranville<sup>c</sup>, R. Guillevin<sup>c,d</sup>, B. Constantin<sup>b</sup>

<sup>a</sup> R&D Scientist at Novadiscovery, 1 place Giovanni da Verrazzano,, Lyon, 69009, France

<sup>b</sup> Université de Poitiers, Laboratoire Signalisation et Transports Ioniques Membranaires, ERL CNRS 7003, Equipe 4CS,, Bâtiment B36 - TSA 51106,1 rue Georges Bonnet, Poitiers Cedex 9, 86073, France

<sup>c</sup> Université de Poitiers, Laboratoire commun I3M (CNRS-UP-CHU-SIEMENS), Laboratoire de Mathématiques et Applications, UMR CNRS 7348, Equipe DACTIM-MIS, Site du Futuroscope - Téléport 2, 11 Boulevard Marie et Pierre Curie, Bâtiment H3 - TSA 61125, Poitiers Cedex 9, 86073, France

<sup>d</sup> CHU de Poitiers, 2 Rue de la Milétrie, Poitiers, 86021, France

## Abstract

A tight control of intracellular [Ca²⁺] 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 Ca²⁺ dynamics using nonlinear dynamics in order to understand Ca²⁺ 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.

**Keywords:** Brain calcium kinetics, ODE non linear system, substrate modeling regularity, well-posedness

## 1. 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 [5, 4], Calcium
(Ca²⁺) signaling is exquisitely poised to link extracellular mechanisms with
intracellular modifications and to control cell fate and many cell functions in-
cluding regulation of metabolism, proliferation, death, gene transcription, cell
migration, exocytosis, and contraction [3]. The pattern of intracellular cal-
cium signaling such as intracellular calcium oscillations that can be repeated
over a longer period [2, 12] is now considered as a universal mechanism of
signal transduction and determines specific cell states; it is also involved in
specific stem cells activity [27] or in various cancer cells hallmarks [8, 26].
Since a tight control of intracellular [Ca²⁺] is essential for the survival and
normal function of cells, resting calcium activity is maintained at a low level
(between 50 nM and 200 nM [1]) in order to keep a large dynamics range
for the calcium signal. Membrane ionic channels and transporters, cytoso-
lic calcium buffers and calcium buffering organelles regulate calcium influx,
storage and extrusion to maintain [Ca²⁺] below the activation thresholds and
extraphysiological values.

The large concentration gradient of Ca²⁺ is created by specialized proteins
named Ca²⁺-ATPases, which are considered in this study, namely, (i) Plasma
membrane Ca²⁺-ATPase (PMCA) located on the plasma membrane, which
extrude Ca²⁺ to the extracellular space or (ii) Sarco/endoplasmic reticulum
Ca²⁺-ATPase (SERCA) inserted in the membrane of the endoplasmic retic-
ulum (ER) that is an important intracellular Ca²⁺ pool. Calcium signal di-
versity and the generation of receptor-triggered Ca²⁺ signals also rely on cell
proteins specialized in calcium transport through cell membranes, which en-
sures the transitory elevation of free cytosolic calcium concentration through
Ca²⁺ flux into the cytoplasm from different compartments. We consider two
main Ca²⁺ fluxes, which could contribute to responses to extracellular signals
and/or intracellular calcium oscillations: (i) a calcium entry along its con-
centration 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²⁺ 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 essen-
tial second messengers triggering both Ca²⁺ 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²⁺ ions from the ER [21, 3], and usually lead to rapid and transient increase of cytosolic Ca²⁺; (ii) DAG, the receptor-operated calcium channel (ROCC), activates plasma membrane which supports entry of Ca²⁺ and participates to slow and sustained increase of cytosolic Ca²⁺ [29, 9]. We also focus on the participation of a second crucial pathway of Ca²⁺ influx, the store-operated calcium channel (SOCC), which is highly dependent of the intracellular Ca²⁺ pool of ER [23] and activates after Ca²⁺ 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²⁺ oscillations (see the reviews Dupond et al., 2011 [11] and Dupond, 2014 [10]). Other studies assume that Ca²⁺ dynamics remains a stochastic process even at the cellular level and that there is poor communication between Ca²⁺ channels, for instance between IP3R Ca²⁺-release channels. Most of the computational models for intracellular Ca²⁺ dynamics focus on the IP3R Ca²⁺-release channels which play a central role in all the cells (Dupond, 2014 [10]). Only some calcium studies give a mathematical analysis (including well-posedness) of their models.

In this study we consider a positive feedback of Ca²⁺ on IP3R, due to the mechanism of Ca²⁺-induced Ca²⁺ release (CICR), and we use the bell-shaped dependency of IP3R activity on cytoplasmic IP3 second messenger (see for example Houart et al. 1999 [16], Lavrentovich and Hemkin 2008 [19]). Some of these models assume that the IP3 concentration is at a constant stimulatory level, as in Borghans et al. (1997) [6]. 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 [14, 28]. Based on other works (see Politi et al., 2006 [25]) 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²⁺ on IP3 metabolism could mediate fluctuations in cellular IP3 levels. We also consider a positive and negative feedbacks of Ca²⁺ on IP3 metabolism: Ca²⁺ activation of PLC producing IP3 and Ca²⁺ 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²⁺ through IP3R calcium release channels also take into account this pumping activity that can be modeled by a Hill expression [10].

Some newly identified players in intracellular Ca²⁺ dynamics have been less often incorporated in computational models, mostly because their kinetics are poorly described quantitatively. As already mentioned, release of stored Ca²⁺ from the ER activates influx through SOCCs, which depend on activation of STIM1 oligomers, sensing the decrease of Ca²⁺ in the ER lumina. SOCCE thus play an important part in the interplay between Ca²⁺ of the ER and Ca²⁺ 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²⁺, the amplitude and duration of Ca²⁺ signals and frequency of Ca²⁺ oscillations. The study of Liu et al. [20] proposed a model that reproduces the steady-state dependence of the SOCCE on the level of Ca²⁺ in the ER lumina, assuming a cooperative binding of ER Ca²⁺ to STIM. Another model, introducing a delay between store emptying and Orai opening [7], has been proposed. To our knowledge no mathematical analysis takes into account at the same time the calcium release and the pumping flux with Ca²⁺ 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²⁺ dynamics using nonlinear dynamics in order to understand Ca²⁺ 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 concentra-
tions. We have chosen to study the combined effect of the IP3-dependent
flux leading to release of ER Ca²⁺ stores, with two different pathways of
Ca²⁺ entry, ROCCE (depending on DAG) and SOCCE (depending on Ca²⁺
ER stores), together with transporters supporting Ca²⁺ uptake into the ER
and Ca²⁺ extrusion through the plasma membrane by two different trans-
porters (PMCA and NCX). We also introduce a continuous Ca²⁺ 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 ac-
tivity of SERCA is inhibited. The IP3-dependent calcium release, ROCCE
and SOCCE, all depend on the activity of PLC. This enzyme is also un-
der 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²⁺ transporters in presence of extracel-
lular stimulation or in absence, considering only the mechanistic interplay
between the membrane transporters.

## 2. Mathematical modeling

The present model is proposed in order to study spontaneous calcium os-
cillations 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 Figure 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.
<figure><img src="image_2.png" /><figcaption>Figure 1: Schematic representation of calcium exchanges in a cell. There are two compartments: the cytosol and the endoplasmic reticulum.</figcaption></figure>

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 $\mu_p$ and a standard deviation of $\sigma_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 concentration and the external calcium concentration. To do so we define the quantity $E_T = \frac{E}{E+k_{E,T}}$ 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 $\alpha_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 $\alpha_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),

$\forall t \in \mathbb{R}^+$,
$$
\begin{align}
C'(t) = & -V_s \frac{C(t)^2}{C(t)^2 + k_s} + V_f(R(t) - C(t)) + V_p \frac{R(t)}{R(t) + k_{r,p}} \frac{I(t)^2}{I(t)^2 + k_{i,p}} \exp\left(-\frac{(C(t) - \mu_p)^2}{2\sigma_p^2}\right) \nonumber \\
& + \frac{V_o}{\pi} \arctan(R_o - R(t)) + \frac{V_o}{2} - V_n \frac{C(t)}{C(t) + k_n} - V_m \frac{C(t)^2}{C(t)^2 + k_m} \tag{1} \\
& + f(t) V_T \frac{D(t)}{D(t) + k_{D,T}} \left(E_T - \frac{C(t)}{C(t) + k_{c,T}}\right), \nonumber
\end{align}
$$

$$
R'(t) = V_s \frac{C(t)^2}{C(t)^2 + k_s} + V_f(C(t) - R(t)) - V_p \frac{R(t)}{R(t) + k_{r,p}} \frac{I(t)^2}{I(t)^2 + k_{i,p}} \exp\left(-\frac{(C(t) - \mu_p)^2}{2\sigma_p^2}\right), \quad (2)
$$

$$
I'(t) = f(t)V_L \frac{C(t)}{C(t) + k_L} - (\alpha_i + V_i \frac{C(t)}{C(t) + k_i})I(t), \quad (3)
$$

$$
D'(t) = f(t)V_L \frac{C(t)}{C(t) + k_L} - (\alpha_d + V_d \frac{C(t)}{C(t) + k_d})D(t), \quad (4)
$$

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:

$$
(C(0), R(0), I(0)) = (\bar{C}, \bar{R}, \bar{I}, \bar{D}) \in (\mathbb{R}^{+} \times \mathbb{R}^{+} \times \mathbb{R}^{+} \times \mathbb{R}^{+}),
$$

**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* ∈ ℝ⁺:

$0 \le f(t) \le F.$

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.
### 3. 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 [17] given by the following expression :

$$S_i = \frac{V(E[C|X_i])}{V(C)}, \quad i = 1 \dots N,$$

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 index 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 Figure 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 Section 5.
<figure><img src="image_3.png" /><figcaption>Figure 2: Sensitivity analysis with first order Sobol indices in time for each model parameter. Displayed parameters are $V_s, V_o, V_m, V_n, V_T$ and $R_o$.</figcaption></figure>

## 4. 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 \in \{1, 2\}$, such as $\forall x \in \mathbb{R}$,

$$T(x) = V \frac{x^i}{x^i + k}$$

* 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 $\mu_p$ and a standard deviation $\sigma_p$ with cytosolic calcium concentration such as $\forall x \in \mathbb{R}$,

$$H_p(x) = \exp(-\frac{(x - \mu_p)^2}{2\sigma_p^2}).$$
* The switch transporter. It describes a transporter with an On/Off behavior towards a substrate. The switch transporter involved in SOCC dynamics has a maximal rate $V_o$ and a value of switch for $R(t) = R_o$ such as $\forall x \in \mathbb{R}$,

$$G_o(x) = \frac{V_o}{\pi} \arctan(R_o - x) + \frac{V_o}{2}.$$

### 4.1. Lemmas on symport like transporters

We define for $V \in \mathbb{R}^+$, $k \in \mathbb{R}^+$, $i \in \{1, 2\} \forall x \in \mathbb{R}$,

$$T_i(x) = V \frac{x^i}{x^i + k}.$$

**Lemma 1.** For the symport transporter of degree 1 (id est $i = 1$) and for $x \in \mathbb{R}^+$, then $T_i$ is $\frac{V}{k}$-Lipschitz continuous, nonnegative and bounded by $V$.

*Proof.* We set, $\forall x \in \mathbb{R}^+$,

$$T_1(x) = V \frac{x}{x + k}.$$

We therefore have $\forall x \in \mathbb{R}^+$,

$$0 \le T_1(x) \le V.$$

Let $x$ and $y$ be in $\mathbb{R}^+$. Then,

$$|T_1(x) - T_1(y)| = V \frac{k|x - y|}{(k + x)(k + y)} \le \frac{V}{k}|x - y|,$$

so that $T_1$ is $\frac{V}{k}$-Lipschitz continuous. $\square$

<figure><img src="image_4.png" /></figure>

**Lemma 2.** For the symport transporter of degree 2 (id est $i = 2$) and for $x \in [0, A]$, $A \in \mathbb{R}^+$ then $T_i$ is $\frac{2VA}{k}$-Lipschitz continuous, nonnegative and bounded by $V$.

*Proof.* We set, $\forall x \in \mathbb{R}^+$,

$$T_2(x) = V \frac{x^2}{x^2 + k}.$$
We therefore have $\forall x \in \mathbb{R}^+$,

$$0 \le T_2(x) \le V.$$

Let $x$ and $y$ be in $[0, A]$, $A \in \mathbb{R}^+$. Then,

$$ \begin{aligned} |T_2(x) - T_2(y)| &= |V \frac{k(x^2 - y^2)}{(x^2 + k)(y^2 + k)}|, \\ &\le \frac{V}{k} |x + y| |x - y|, \\ &\le \frac{2VA}{k} |x - y|, \end{aligned} $$

so that $T_2$ is $2VAk$-Lipschitz continuous. $\square$

In Figure 3 we give examples of symport like functions.

<figure><img src="image_5.png" /><figcaption>Figure 3: Examples of symport like functions $T_i$ with $V = 1$ and $k = 3$.</figcaption></figure>
### 4.2. Study of the bell-shaped transporter

We define for $\mu_p \in \mathbb{R}^+$, $\sigma_p \in \mathbb{R}^+ \forall x \in \mathbb{R}$,

$$H_p(x) = \exp\left(-\frac{(x - \mu_p)^2}{2\sigma_p^2}\right).$$

**Lemma 3.** For $x \in \mathbb{R}$ then $H_p$ is Lipschitz continuous, nonnegative and bounded by 1.

*Proof.* We set $\forall x \in \mathbb{R}$,

$$H_p(x) = \exp\left(-\frac{(x - \mu_p)^2}{2\sigma_p^2}\right).$$

Then its derivative is $\forall x \in \mathbb{R}$,

$$H'_p(x) = \frac{-1}{\sigma_p^2}(x - \mu_p) \exp\left(-\frac{(x - \mu_p)^2}{2\sigma_p^2}\right),$$

and the second derivative is $\forall x \in \mathbb{R}$,

$$H'_p(x) = \frac{1}{\sigma_p^2}(\frac{1}{\sigma_p^2}(x - \mu_p)^2 - 1) \exp\left(-\frac{(x - \mu_p)^2}{2\sigma_p^2}\right).$$

We therefore have the following variation table.
<table><thead><tr><th>x</th><th>-∞</th><th>μ<sub>p</sub> - σ<sub>p</sub></th><th>μ<sub>p</sub></th><th>μ<sub>p</sub> + σ<sub>p</sub></th><th>+∞</th></tr></thead><tbody><tr><th>H''<sub>p</sub>(x)</th><td>+</td><td>-</td><td>-</td><td>+</td><td></td></tr><tr><th>H'<sub>p</sub>(x)</th><td>0<sup>+</sup></td><td>1/σ<sub>p</sub> exp(-1/2)</td><td>0</td><td>-1/σ<sub>p</sub> exp(-1/2)</td><td>0<sup>-</sup></td></tr><tr><th>H'<sub>p</sub>(x)</th><td>+</td><td>+</td><td>-</td><td>-</td><td></td></tr><tr><th>H<sub>p</sub></th><td>0<sup>+</sup></td><td>1</td><td>0<sup>-</sup></td><td>0<sup>-</sup></td><td></td></tr></tbody></table>

The derivative $H'_p$ is bounded by $\frac{1}{\sigma_p} \exp(-\frac{1}{2})$ and therefore $H_p$ is Lipschitz continuous. □

In Figure 4 we give an example of a bell-shaped $H_p$ function.

### 4.3. Study of the switch transporter

We define for $V_0 \in \mathbb{R}^+$, $R_o \mathbb{R}^+ \forall x \in \mathbb{R}$,

$$G_o(x) = \frac{V_o}{\pi} \arctan(R_o - x) + \frac{V_o}{2}.$$

**Lemma 4.** For $x \in \mathbb{R}$ $G_o$ is $\frac{V_o}{\pi}$ Lipschitz continuous, nonnegative and bounded by $V_o$.

*Proof.* We set $\forall x \in \mathbb{R}$,

$$G_o(x) = \frac{V_o}{\pi} \arctan(R_o - x) + \frac{V_o}{2}.$$

Then its derivative is $\forall x \in \mathbb{R}$,

$$G'_o(x) = \frac{-V_o}{\pi} \frac{1}{1 + (R_o - x)^2}$$
<figure><img src="image_6.png" /><figcaption>Figure 4: Example of a bell-shaped function $H_p$ with $\mu_p = 5$ and $\mu_p = 0.3$.</figcaption></figure>

and we have $\forall x \in \mathbb{R}$,

$$0 \le G'(x) \le \frac{V_o}{\pi}$$

Therefore $G_0$ is $\frac{V_o}{\pi}$-Lipschitz continuous and we have the following variation table.

<table><tr><td>x</td><td>-∞</td><td>+∞</td></tr><tr><td>G'<sub>o</sub>(x)</td><td colspan="2">-</td></tr><tr><td>G<sub>o</sub>(x)</td><td>V<sub>o</sub></td><td>0+</td></tr></table>

This function is therefore nonnegative with an upper bound $V_0$.

<figure><img src="image_7.png" /></figure>

In Figure 5 we give an example of a bell-shaped $G_o$ function.
<figure><img src="image_8.png" /><figcaption>Figure 5: Example of a switch function $G_0$ with $V_o = 0.3$ and $R_o = 300$.</figcaption></figure>

### 4.4. Cross-results

**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 $\exists F \in \mathbb{R}^+/|f| \le F$. We also assume that the function $g$ $k_g$-Lipschitz such that $\exists G \in \mathbb{R}^+/|g| \le G$.

Let $x, y \in \mathbb{R}$. We have

$$|(fg)(x) - (fg)(y)| \le |f(x)||g(x) - g(y)| + |g(y)||f(x) - f(y)|, \\ \le Fk_g|x - y| + Gk_f|x - y|, \\ \le \max(Fk_g, Gk_f)|x - y|,$$

and therefore the function $fg$ is Lipschitz continuous. □

Results in this preliminary work have two aims : giving a first idea of each transporter dynamics and get mathematical results which are useful in view of the mathematical analysis.
# 5. Mathematical analysis

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

## 5.1. Nonnegativity, existence and bounds

### 5.1.1. Nonnegativity

Recall that an ODE system $x'(t) = f(x(t))$ on $\mathbb{R}^n$, $x = [x_1, \dots, x_n]$, $f = [f_1, \dots, f_n]$ is called quasipositive if the condition :

$$x \ge 0, x_k = 0 \Rightarrow f_k(x) \ge 0$$

is verified for all $k = 1, \dots, n$. Using Section 4, we know that system (1)-(4) is quasipositive. Hence solutions with nonnegative initial data $(\bar{C}; \bar{R}, \bar{I}; \bar{D})$ remain in $(\mathbb{R}^+)^4$ for all positive times.

### 5.1.2. Upper bounds

From 1-4, we have the following inequalities :

$$R'(t) \le V_s + V_f(C(t) - R(t)), \qquad (5)$$

$$C'(t) \le V_f(R(t) - C(t)) + V_p + V_o + FV_T, \qquad (6)$$

$$I'(t) \le FV_L - \alpha_i I(t), \qquad (7)$$

$$D'(t) \le FV_L - \alpha_d D(t). \qquad (8)$$

Inequalities 7 and 8 imply, using Gronwall's lemma,

$$I(t) \le I(0) + \frac{FV_L}{\alpha_i} := B_I,$$

and

$$D(t) \le D(0) + \frac{FV_L}{\alpha_d} := B_D,$$

**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 :

$$V_n + V_m \ge FV_T E_T + V_o$$

In that case, $\exists B_C, B_R \in \mathbb{R}^{+*}$

$$C \le B_C, \qquad (9)$$

$$R \le B_R. \qquad (10)$$
Proof. Setting $\forall x \in \mathbb{R}^+$

$$u(x) = V_o - V_n \frac{x}{x + k_n} - V_m \frac{x^2}{x^2 + k_m} + FV_T E_T$$

and using 1-4 we have

$$R'(t) + C'(t) = \frac{V_o}{\pi} \arctan(R_o - R(t)) + \frac{V_o}{2} - V_n \frac{C(t)}{C(t) + k_n} - V_m \frac{C(t)}{C(t) + k_m} \\ + f(t)V_T(E_T - \frac{C(t)}{C(t) + k_{c,T}}) \qquad (11)$$

$$\le u(C(t)). \qquad (12)$$

We want to study $u$ variations for $x \in \mathbb{R}^+$. Its derivative reads

$$u'(x) = -V_n \frac{k_n}{(x + k_n)^2} - V_m \frac{2xk_m}{(x^2 + k_m)^2} < 0,$$

so that $u$ is a strictly decreasing function and

<table><tr><td>x</td><td>0</td><td>+&infin;</td></tr><tr><td>u'(x)</td><td colspan="2">&minus;</td></tr><tr><td>u(x)</td><td colspan="2">V<sub>o</sub> + FV<sub>T</sub>E<sub>T</sub><br/>Vo + FV<sub>T</sub>E<sub>T</sub> &minus; (V<sub>n</sub> + V<sub>m</sub>).</td></tr></table>

Therefore, assuming that $V_n + V_m \ge FV_T E_T + V_o$, we have, using the intermediate value theorem, $\exists! y > 0$ such that $u(y) = 0$. Moreover for all $x \in \mathbb{R}^+$, $x > y$, we have $u(x) < u(y) = 0$. If this condition is not satisfied, $u$ is nonnegative.

Therefore using Equation 12 and assuming $V_n + V_m \ge 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 $\forall t \in \mathbb{R}^+$

$$V_s \frac{C(t)}{C(t) + k_s} + V_f(C(t) - R(t)) < V_p \frac{R(t)}{R(t) + k_{r,p}} \frac{I(t)}{I(t) + k_{i,p}} \exp(-\frac{(C - \mu_p)^2}{2\sigma_p^2}),$$

so that, using 1

$$C'(t) < \frac{V_o}{\pi} \arctan(R_o - R(t)) + \frac{V_o}{2} - V_n \frac{C(t)}{C(t) + k_n} - V_m \frac{C(t)}{C(t) + k_m} \\ + f(t) V_T \frac{D(t)}{D(t) + k_{D,T}} (E_T - \frac{C(t)}{C(t) + k_{c,T}}), \quad (13) \\ < V_o - V_n - V_m - FV_T E_T.$$

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

$$C(t) < \max(N_c, C(0)) := B_C.$$

Using 5, this implies that

$$R'(t) < V_s + V_f(B_c - R(t)),$$

so that, using Gronwall lemma,

$$R(t) \le R(0) + \frac{V_s}{V_f} + B_c := B_R,$$

and therefore $R$ and $C$ are bounded. $\square$

**Remark 1.** Lemma 6 ensures bounds on $C$ and $R$ under the condition

$$V_n + V_m \ge FV_T E_T + V_o.$$

This result has to be compared with Sobol indices results in Section 3. It is assumed that $V_n, V_m, V_o$ and $V_T$ have a strong impact on the dynamics in both studies.
### 5.1.3. Existence of the solution

Since we have nonnegativity and bounds of C, R, D and I then, based on Section 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

$$X(t) := (C(t); R(t); I(t); D(t)),$$

to have $\forall t \in \mathbb{R}^+$:

$$X'(t) = H(t, X(t)), \quad X(0) = X_0,$$

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 $\forall t \in \mathbb{R}^+$.

## 5.2. Steady-state study

Assuming $f = 0$, the steady-state of system 1-4 is given by

$$\begin{aligned} \tilde{I} &= 0, \\ \tilde{D} &= 0, \\ \tilde{R} &= \frac{V_s}{V_f} \frac{\tilde{C}^2}{\tilde{C}^2 + k_s} + \tilde{C}, \\ \frac{V_o}{\pi} \arctan(R_o - \tilde{R}) + \frac{V_o}{2} &= V_n \frac{\tilde{C}}{\tilde{C} + k_n} + V_m \frac{\tilde{C}^2}{\tilde{C}^2 + k_m}. \end{aligned}$$

We define the function S such that, $\forall x \in \mathbb{R}^2$,

$$S(x) = \frac{V_o}{\pi} \arctan\left(R_o - \frac{V_s}{V_f} \frac{x^2}{x^2 + k_s} - x\right) + \frac{V_o}{2} - V_n \frac{x}{x + k_n} - V_m \frac{x^2}{x^2 + k_m}.$$

We are looking for a solution at $E(\tilde{C}) = 0$. The related derivative reads

$$S'(x) = \frac{V_o}{\pi} \frac{-2k_s x \frac{V_s}{V_f} - (x^2 + k_s)^2}{(x^2 + k_s)^2 \left(1 + \left(R_o - \frac{V_s}{V_f} \frac{x^2}{x^2 + k_s} - x\right)^2\right)} - V_n \frac{k_n}{(x + k_n)^2} - V_m \frac{2k_m x}{(x^2 + k_m)^2} < 0.$$

We therefore have
<table><tr><td>x</td><td>0</td><td>+∞</td></tr><tr><td>S'(x)</td><td colspan="2">-</td></tr><tr><td>S(x)</td><td>V<sub>0</sub></td><td>0</td></tr></table>

Using the intermediate value theorem, ∃! $x_c > 0$ such that $S(x_c) = 0$. The related steady-state is therefore given by

$$ 
\begin{aligned}
\tilde{I} &= 0, \\
\tilde{D} &= 0, \\
\tilde{C} &= x_c, \\
\tilde{R} &= \frac{V_s}{V_f} \frac{x_c^2}{2 + k_s} + x_c.
\end{aligned}
 $$

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 steady-state is always included on the viability domain defined by the bounds. We know that $\bar{C}$ and $\bar{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 \le B_c = \max(C(0), N_c)$. To do so it is sufficient to prove that $x_c \le N_c$. We also want to prove that $\tilde{R} < B_r$.

**Remark 2.** This verification is only possible when upper bounds are exhibited, that is to say when $V_n + V_m \ge 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.
We know that, $\forall x \in \mathbb{R}$,

$$u(x) - S(x) = FV_T E_T + V_o \left(1 - \frac{1}{\pi} \arctan\left(R_o - \frac{V_s}{V_f} \frac{x^2}{x^2 + k_s} - \frac{1}{2}\right)\right),$$
$$u(x) > S(x).$$

Hence, recalling that $S(x_c) = 0$,

$$u(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

$$x_c < N_c \le B_c.$$

Hence $\bar{C}$ is always on the viability domain of $C$.

We will now prove that $\tilde{R} < B_r$. We know that

$$B_r = R(0) + \frac{V_s}{V_f} + B_c,$$

so that we have

$$\begin{aligned} \tilde{R} &= \frac{V_s}{V_f} \frac{x_c^2}{x_c^2 + k_s} + x_c, \\ &\le \frac{V_s}{V_f} + x_c, \\ &\le \frac{V_s}{V_f} + B_c, \\ &\le B_r, \end{aligned}$$

and $\bar{R}$ is always on the viability domain of $R$.

## 6. 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.

## 7. 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.

<figure><img src="image_9.png" /><figcaption>Figure 6: Representative cytosolic calcium oscillations. Experiments were done in human glioma stem cells with the calcium sensitive probe Fura2.</figcaption></figure>
The Figure shows the oscillation recorded during 60 to 90 seconds. 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.* [6] and Keener and Sneyd [18]. 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 [22]. 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.

* 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 *et al.* [19] adding the results of

* Politi *et al.* [25] on the PMCA,

* Behringer *et al.* [1] on the IP3R,

* Haeri *et al.* [13] on the cytosolic calcium impact on the IP3R transporter,

* Hofer *et al.* [15] on the IP3 dynamics.

These values help to find accurate model dynamics.
7.1. *Mechanical periodicity*

To test our model, we provide simulations setting *f* = *F*. In Figure 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

<figure><img src="image_10.png" /><figcaption>Figure 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.</figcaption></figure>

All the concentrations remain nonnegative through time. On the one
hand and since they have the same production rate and comparable degra-
dation 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.
This model supports the idea that spontaneous Ca²⁺ oscillations can be
generated without the aid of external stimulation [24]. However related os-
cillations occur in long time and look different that recordings obtained in
living cells, reported by Figure 6.
<figure><img src="image_11.png" /><figcaption>Figure 8: Portrait phases of calcium concentrations on the left and DAG/IP3 concentrations on the right. While DAG and IP3 increase (or decrease) at the same time, calciums concentration vary roughly in turn.</figcaption></figure>

To figure out how these mechanical periods may vary, we test several values of some maximal fluxes : $V_n$ for the NCX (Figure 9), $V_m$ for the PMCA (Figure 10), $V_T$ for the ROCC (Figure 11) and $V_f$ for the leak flux (Figure 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 Figure 14 and 15.
<figure><img src="image_12.png" /><figcaption>Figure 9: Trajectories of cytosolic calcium, ER calcium, IP3 and DAG with no external stimulation (F constant) and different values of Vn. Three kinds of dynamics set appart : the one for Vn < 0.25 having a "high" stable state, the one for 0.25 < Vn < 0.35 with a periodic behavior and the one for Vn > 0.35 having a "low" stable state. Figure fig:mecal gives reference parameters</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 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:mecal gives reference parameters</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 11: Trajectories of cytosolic calcium, ER calcium, IP3 and DAG with no external stimulation (F constant) and different values of V<sub>o</sub>. Again three kinds of dynamics set appart : the one for V<sub>o</sub> < 0.03 having a "low" stable state, the one for 0.03 < V<sub>o</sub> < 0.07 with a periodic behavior and the one for V<sub>o</sub> > 0.07 having a "high" stable state. Figure fig:mecal gives reference parameters</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 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:mecal gives reference parameters</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 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:mecal gives reference parameters</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 14: Trajectories of cytosolic 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:mecal gives reference parameters</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 15: Trajectories of ER calcium with no external stimulation (F constant) and different values of V<sub>s</sub> and k<sub>s</sub>. Not only is the value of V<sub>s</sub> by itself important but also the union of all given transport parameters (V<sub>s</sub> but also k<sub>s</sub> for the SERCA). Figure fig:mecal gives reference parameters</figcaption></figure>
Results from Figure 9 to 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 Section 3.

## 7.2. 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 $\forall t \in \mathbb{R}^+$,

$$f(t) = \frac{F}{2}(\cos(freq * t) + 1).$$

This function $f$ is periodic, nonnegative and bounded by $F$. It describes external proteic variations. It is displayed in Figure 16. In Figure 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.
<figure><img src="image_19.png" /><figcaption>Figure 16: External stimulation impacting ROCC and PLC dynamics given by the function $f$ for $F = 1.5$ and $freq = 0.42$.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 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 $\sim 60s = 1m$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 18: Portrait phases of calcium concentrations on the left and DAG/IP3 concentrations on the right. While DAG and IP3 increase (or decrease) at the same time, calciums concentration have a complex variation with two periods.</figcaption></figure>
Cytosolic calcium displays two kinds of periods. These results are in ac-
cordance with literature data (see Section 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 dif-
ferent values of *freq*.

<figure><img src="image_22.png" /><figcaption>Figure 19: Trajectories of ER calcium with external stimulation (*f* periodic) and different frequencies *freq*. Three kinds of dynamics set appart : the one for *freq* < 0.035 with smooth oscillations, the one with *freq* ∈ [0.35, 0.45] with two periods and the ones for *freq* > 0.45 with sharped oscillations.</figcaption></figure>
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 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*.

<figure><img src="image_23.png" /><figcaption>Figure 20: Example of cytosolic calcium trajectory with different values of *freq*. Figure fig:osc2 gives reference parameters.</figcaption></figure>
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 Figure 21 displays what happens when $f = 0$. In that case one can observe a steady-state as predicted by the mathematical analysis.

<figure><img src="image_24.png" /><figcaption>Figure 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 500s) while calcium concentrations reach more slowly a non-zero steady-state (time scale 100 000s ~ 28h). Figure 17 gives reference parameters</figcaption></figure>

# 8. 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 [18], 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 sepa-
rately. We then study their combined effect on calcium concentrations. We
obtain existence, uniqueness and bounds on the solutions for the related sys-
tem. 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.

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 [6, 18].

3. The periodic behavior is only true for a narrow set of parameters in-
cluding $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 pa-
rameter value outside this window of values makes the system reach a
steady-state.

4. Spontaneous periodic variations of calcium can be obtained in the cyto-
plasm 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 peri-
odic 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 be-
tween 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 steady-state. 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 tranport rate $R_o$, the maximal calicum, IP3 and DAG concentrations are not the same and possibly limit oscillations (Figure 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.

## References

[1] Behringer E. J., Scallan J. P., Jafarnejad M., Castorena-Gonzalez J. A., Zawieja S. D., Moore J. E., Michael J. Davis & Segal S. S. *Calcium and electrical dynamics in lymphatic endothelium*. The Journal of physiology, (2017), 595(24), p. 7347-7368.

[2] Berridge M. J. *Calcium oscillations*. *Journal of Biological Chemistry*, (1990), 265(17), p. 9583-9586
[3] Berridge M. J. *The inositol trisphosphate/calcium signaling pathway in health and disease*. Physiological reviews, (2016), 96(4), p. 1261-1296.

[4] Berridge M. J., Bootman M. D. & Roderick H. L. *Calcium signalling: dynamics, homeostasis and remodelling*. Nature reviews Molecular cell biology, (2003), 4(7), p. 517-529.

[5] Berridge M. J., Lipp P. & Bootman M. D. *The versatility and universality of calcium signalling*. Nature reviews Molecular cell biology, (2000), 1(1), p. 11-21.

[6] Borghans J. M., Dupont G. & Goldbeter A. *Complex intracellular calcium oscillations: A theoretical exploration of possible mechanisms*. Bio-physical chemistry, (1997),66(1), p. 25-41.

[7] Croisier H., Tan X., Perez-Zoghbi J. F., Sanderson M. J., Sneyd J. & Brook B. S. *Activation of store-operated calcium entry in airway smooth muscle cells: insight from a mathematical model*. PloS one, (2013), 8(7).

[8] Déliot N. & Constantin, B. *Plasma membrane calcium channels in cancer: Alterations and consequences for cell proliferation and migration*. Biochimica et Biophysica Acta (BBA)-Biomembranes, (2015), 1848(10), p. 2512-2522.

[9] Dietrich A., Kalwa H, Rost BR & Gudermann T. *The diacylglycerol-sensitive TRPC3/6/7 subfamily of cation channels: functional characterization and physiological relevance*. Pflügers Arch, (2005), 451, p. 72-80.

[10] Dupont G. *Modeling the intracellular organization of calcium signaling*. Wiley Interdisciplinary Reviews: Systems Biology and Medicine, (2014), 6(3), p. 227-237.

[11] Dupont G., Combettes L., Bird G. S., & Putney J. W. *Calcium oscillations*. Cold Spring Harbor perspectives in biology, (2011), 3(3).

[12] Fewtrell C. $Ca^{2+}$ oscillations in non-excitable cells. Annual Review of Physiology, (1993), 55(1), p. 427-454.

[13] Haeri H. H., Hashemianzadeh, S. M. & Monajjemi M. *A kinetic Monte Carlo simulation study of inositol 1, 4, 5-trisphosphate receptor (IP3R)*
*calcium release channel*. Computational biology and chemistry, (2007), 31(2), 99-109.

[14] Hirose K., Kadowaki S., Tanabe M., Takeshima H., & Iino, M. *Spatiotemporal dynamics of inositol 1, 4, 5-trisphosphate that underlies complex Ca2+ mobilization patterns*. Science, (1999), 284(5419), p. 1527-1530.

[15] Höfer T., Venance L., & Giaume C. *Control and plasticity of intercellular calcium waves in astrocytes: a modeling approach*. Journal of Neuroscience, (2002), 22(12), p. 4850-4859.

[16] Houart G., Dupont G. & Goldbeter A. *Bursting, chaos and birhythmicity originating from self-modulation of the inositol 1, 4, 5-trisphosphate signal in a model for intracellular Ca2+ oscillations*. Bulletin of mathematical biology, (1999), 61(3), p. 507-530.

[17] Iooss B, Lemaître P. *A review on global sensitivity analysis methods*. Uncertainty Management in Simulation-Optimization of Complex Systems. Springer (2015). p. 101-122.

[18] Keener J. P. & Sneyd J. *Mathematical physiology 1: Cellular physiology*. (2009).

[19] Lavrentovich M. & Hemkin S. *A mathematical model of spontaneous calcium (II) oscillations in astrocytes*. Journal of Theoretical Biology, (2008), 251(4), p. 553-560.

[20] Liu W., Low N. W., Feng B., Wang G. & Diniz da Costa J. C. *Calcium precursors for the production of CaO sorbents for multicycle CO2 capture*. Environmental science & technology, (2010), 44(2), p. 841-847.

[21] Mikoshiba, K. *Role of IP3 receptor signaling in cell functions and diseases*. Advances in biological regulation, (2015), 57, p. 217-227.

[22] Parekh A. B. *Decoding cytosolic Ca2+ oscillations*. Trends in biochemical sciences, (2011), 36(2), p. 78-87.

[23] Parekh A. B. & Putney Jr, J. W. *Store-operated calcium channels*. Physiological reviews, (2005), 85(2), p. 757-810.
[24] Parri H. R. & Crunelli V. *The role of Ca²⁺ in the generation of spontaneous astrocytic Ca²⁺ oscillations*. Neuroscience, (2003), 120(4), p. 979-992.

[25] Politi A., Gaspers L. D., Thomas A. P. & Höfer, T *Models of IP3 and Ca2+ oscillations: frequency encoding and identification of underlying feedbacks*. Biophysical journal, (2016),90(9), p. 3120-3133.

[26] Prevarskaya N., Skryma R. & Shuba Y. (2018). *Ion channels in cancer: are cancer hallmarks oncochannelopathies?* (2018), Physiological reviews, 98(2), p. 559-621.

[27] Terrié E., Coronas V., & Constantin B. *Role of the calcium toolkit in cancer stem cells*. Cell calcium, (2019), 80, p. 141-151.

[28] Thore S., Dyachok O., & Tengholm A. *Oscillations of phospholipase C activity triggered by depolarization and Ca2+ influx in insulin-secreting cells*. Journal of Biological Chemistry, (2004), 279(19), p. 19396-19400.

[29] Vazquez G., Lievremont J. P., Bird G. S. J. & Putney J. W. *Human Trp3 forms both inositol trisphosphate receptor-dependent and receptor-independent store-operated cation channels in DT40 avian B lymphocytes*. Proceedings of the National Academy of Sciences, (2001), 98(20), p. 11777-11782.

[30] Young R., Wu F., Van de Water N., Ames R., Gamble G., & Reid I. R. *Calcium sensing receptor gene A986S polymorphism and responsiveness to calcium supplementation in postmenopausal women*. The Journal of Clinical Endocrinology & Metabolism, (2003), 88(2), p. 697-700.