# Multiscale Modeling Indicates That Temperature Dependent [Ca2+]i Spiking in Astrocytes Is Quantitatively Consistent with Modulated SERCA Activity

Authors: Niko Komin; Mahsa Moein; Mark H. Ellisman; Alexander Skupin
Journal: Neural Plasticity
 (2015)

---

## Abstract

Changes in the cytosolic Ca2+ concentration ([Ca2+]i) are the most predominant active signaling mechanism in astrocytes that can modulate neuronal activity and is assumed to influence neuronal plasticity. Although Ca2+ signaling in astrocytes has been intensively studied in the past, our understanding of the signaling mechanism and its impact on tissue level is still incomplete. Here we revisit our previously published data on the strong temperature dependence of Ca2+ signals in both cultured primary astrocytes and astrocytes in acute brain slices of mice. We apply multiscale modeling to test the hypothesis that the temperature dependent [Ca2+]i spiking is mainly caused by the increased activity of the sarcoendoplasmic reticulum ATPases (SERCAs) that remove Ca2+ from the cytosol into the endoplasmic reticulum. Quantitative comparison of experimental data with multiscale simulations supports the SERCA activity hypothesis. Further analysis of multiscale modeling and traditional rate equations indicates that the experimental observations are a spatial phenomenon where increasing pump strength leads to a decoupling of Ca2+ release sites and subsequently to vanishing [Ca2+]i spikes.


---

## Body

## 1. Introduction

The human brain is by far the most complex organ; it hosts our consciences, is the source of our inspiration, and controls many of our vegetative functions. The extraordinary complexity of cellular interactions in terms of structural as well as dynamic properties enables the huge flexibility of our behavior and the amazing plasticity required for learning [1]. While traditional neuroscience has focused mainly on neurons for decades, the importance of the more abundant glia cells has only become evident more recently [2, 3].Glia cells are generally known for their supportive function for neurons. The minority of glia cells are oligodendrocytes that myelinate axons [4]. More recently, oligodendrocytes were also found to essentially support axon integrity by metabolic coupling through monocarboxylate transporters (MCTs) [5, 6]. Microglia represent the macrophages of the brain and play an important role in immune responses as well as in modifying the neuronal network structure by synapse pruning [7]. The majority of glia cells are astrocytes that are commonly considered as metabolic supporters of neurons. On the one hand astrocytes regulate the blood flow in response to neuronal activity [8]. Astrocytes take up approximately 80% of the metabolized glucose from the blood and perform glycolysis [9]. According to the astrocyte-neuron lactate shuttle (ANLS) hypothesis [10], astrocytes generate lactate from pyruvate, the end product of glycolysis, by the lactate dehydrogenase (LDH) and export it into the extracellular space by MCTs. From the extracellular space, neurons take up the lactate and generate pyruvate by the reversible LDH reaction. The neuronal pyruvate is subsequently used as energy substrate for the TCA cycle in mitochondria to generate adenosine triphosphate (ATP) that is extensively needed to run the ion pumps in neurons to keep the large electrochemical potential across the plasma membrane. On the other hand, astrocytes play also an important function for neuronal proteostasis as they remove degraded material from the tissue either directly by intercellular transport towards the blood stream or by controlling the cerebral fluid flux [11].Besides their supportive role in metabolism, astrocytes are now also known to modulate neuronal activity through active signaling mechanisms. Astrocytes are sealing up the dyadic cleft facilitating neurotransmitter diffusion and are essential for glutamate uptake and recycling [12]. Since astrocytes express also glutamate receptors, they can be stimulated by neuronal activity. Upon stimulation, astrocytes can secrete glutamate as well as ATP into the extracellular space including the dyadic cleft and thereby induce a local amplification mechanism. This neuron-astrocyte crosstalk has led to the picture of the tripartite synapse [13, 14].The most predominant and studied signaling mechanism in astrocytes is inositol 1,4,5-trisphosphate (IP3) mediated Ca2+ signaling. In general, Ca2+ signaling is a versatile and universal second messenger that translates extracellular signals into intracellular responses like molecular motor activation or modulation of gene expression [15–17]. Within the IP3 pathway, extracellular agonists are detected by plasma membrane receptors. Upon binding, these G-protein coupled receptors activate the phospholipase C (PLC) that produces IP3 at the plasma membrane from which it diffuses into the cytosol. There IP3 can bind to IP3 receptors (IP3Rs) localized in the membrane of the endoplasmic reticulum (ER). When IP3 and Ca2+ are bound to the receptor, it can be open and Ca2+ is released from the ER into the cytosol due to the large concentration differences between these two compartments of typically several orders of magnitude. From the cytosol the released Ca2+ is moved back into the extracellular space by plasma membrane calcium ATPase (PMCA) and into the ER by sacroendoplasmic reticulum calcium ATPases (SERCAs).The collection of the Ca2+ tool kit [18] can generate a wide spectrum of cellular responses ranging from long lasting high cytosolic Ca2+ levels ([Ca2+]i) to single [Ca2+]i transients. The most common Ca2+ signals in astrocytes are repetitive cytosolic [Ca2+]i transients often referred to as [Ca2+]i oscillations [18] that are typically elicited by neurotransmitters and hormones but astrocytes do also exhibit spontaneous [Ca2+]i spiking.This oscillation-like behavior is rather common for many cell and stimulation types and has led to the general assumption that Ca2+ signaling is frequency coded [17]. This perspective is challenged by more recent studies including ours that have demonstrated that [Ca2+]i spiking is indeed stochastic and therefore not truly oscillatory [19–21]. By systematic analysis of single cell [Ca2+]i time courses for hundreds of cells of different cell types and stimulation antagonists, we found a preserved linear relation between the average period (Tav) and the standard deviation (σ) of the inter-[Ca2+]i-spike intervals. Applying methods from the theory of stochastic processes, we have shown that this pathway specific linear σ-Tav relation characterizes the internal feedback strength [22] and represents an estimator for the information content [20, 23].To understand the random signaling mechanism and underlying principles, we developed a computational multiscale model that takes into account the compartmentalization of the cell into cytosol and ER, the stochastic binding of signaling molecules to individual IP3Rs, their spatial organization into channel clusters, and the spatial diffusion of Ca2+ and Ca2+ binding buffers [27]. The enabled mechanistic simulations have shown that the experimental observations are in accordance with the implemented hierarchical stochastic system where single channel fluctuations are carried onto the level of the cell by diffusion mediated Ca2+ induced Ca2+ release (CICR). The modeling based insights on the functional robustness have guided further experiments that have revealed a fold change behavior in Ca2+ signaling of stimulated HEK cells and hepatocytes and were used to develop and validate a universal encoding relation [28].The establishment of the signaling mechanism was driven to a large extent by experiments with astrocytes and microglia [21]. In particular the spontaneous activity of cultured astrocytes was exploited to obtain a baseline activity of the Ca2+ signaling machinery without potential antagonist induced feedback mechanisms. The astrocytic base level activity is characterized by a slope of 1 in the linear σ-Tav relation and corresponds to the theoretical minimal information gain. This is in accordance with the general assumption that stimulation should carry information that would correspond to a decreased slope of the σ-Tav relation and is observed in stimulated cells [20].The spontaneous Ca2+ transients in astrocytes occur independently of neuronal activity but are modulated by neuronal function like induced epileptic activity [29]. The astrocytic Ca2+ increase can in turn be sensed by neurons because it has been correlated with release of neuroactive substances [30, 31]. Thereby astrocytes can release glutamate [32, 33], D-Serine [34, 35], ATP [36–39], chemokines [40], and nitric oxide (NO) [41]. The Ca2+ signal forms are rather heterogeneous and include random profiles and rhythmic oscillations as well as bursting activity.Spontaneous[Ca2+]i activity has been reported for astrocytes in culture [42–44] and acutely isolated slices of several brain regions [31, 45–47]. Spontaneous astrocytic Ca2+ signaling is modulated by pathological events such as those in reactive astrocytes surrounding stab wounds in the neocortex that lack spontaneous Ca2+ activity or astrocytes in the neocortex under epileptic conditions that show more spontaneously active astrocytes [45].Despite all these insights, the common mechanism of this spontaneous activity is still not fully understood. In particular, the higher abundance of spontaneous[Ca2+]i activity in cultured astrocytes compared to astrocytes within acute brain slices is only barely investigated [24]. We have shown in a comparison study that these differences can be explained to a large extent by the different temperatures at which the astrocytic [Ca2+]i activity is typically measured. Here we revisit the experimental data and integrate it with our mechanistic multiscale modeling framework to quantitatively test our hypothesis on the molecular foundation of the temperature dependency.


## 2. Results

Motivated by the typically observed and reported significant differences in Ca2+ singling between cultured astrocytes and astrocytes in acute brain slices, we reconsidered the experimental workflow [24]. The most dramatic difference between the two experimental protocols is the different temperatures these two cell systems are typically measured at [21, 48]. In thermodynamics the general impact of temperature T on an individual reaction rate k is well approximated by the Arrhenius law by k=Ae-Eact/RT, where A denotes the approximate reaction rate, Eact describes the activation energy of the reaction, and R is the universal gas constant [49]. While the Arrhenius law describes the temperature dependency of an isolated reaction in a well-stirred reactor, the impact of temperature in more complex reaction systems may differ from the simple exponential relation. Alternatively it might be possible that the thermal energy in the measured temperature range is of the order of the activation energy, where the Arrhenius law shows a linear dependency. To explore if temperature could explain the Ca2+ signaling phenomenon, we summarize our careful analysis of temperature induced [Ca2+]i dynamics changes and use here computational modeling to investigate potential reasons for this phenomenon.


### 2.1. Experiments Show Strong Temperature Dependency in Cortical and Cultured Astrocytes

In the experimental study [24] we used cortical slices obtained from P 10–P 12 mouse brains and cultured astrocytes. The temperature was either changed abruptly or smoothly either by changing the heating temperature of the water bath for the perfusion buffer or by switching from one perfusion buffer to another, respectively. Thereby the temperature was measured within the recording chamber close to the sample cover slip.


#### 2.1.1. Higher Temperatures Reversibly Suppress [Ca2+]i Spiking Instantaneously

Starting from the standard protocols, we first measured Ca2+ activity at room temperature and separately at more physiological temperatures with Fluo4-AM loaded to the cells before (more detailed experimental description can be found in [24]). In cortical slices obtained from P 10–P 12 mouse brains measured at room temperature, astrocytes showed robust spontaneous signaling activity defined by long durations and large amplitudes of [Ca2+]i spikes. In contrast, when slices were measured at physiological temperatures, only very few cells showed spontaneous activity.To compare this behavior with primary astrocytes cultured from mouse neocortex, we measured Ca2+ activity at room temperature (20°C) and then abruptly changed the perfusion buffer to one of 28°C or of 30°C in another setup. For these cultured astrocytes we observed a similar significant decrease of cells exhibiting Ca2+ signals already within the first minute after the temperature shift (Figure 1). Out of the active cells under 20°C (72 ± 24% of all cells in 4 independent experiments) more than a quarter of cells have not shown any [Ca2+]i spikes under 28°C and under 30°C only 20% exhibited at least a single [Ca2+]i spike. Remarkably, this dramatic decrease was fully reversible when perfusing the cells again with buffer at room temperature. Analogous results about the reversibility were also found in experiments with acute brain slices and the reversibility was independent of the experimental order of temperature changes.Figure 1
Summary of experimental observations on temperature dependent[Ca2+]i activity [24]. (a) Fluorescence recordings from three individual astrocytes in three different experiments. The temperature was changed between 20°C and 30°C by using heated and nonheated perfusion buffer (indicated by the bar). Fluorescent amplitude ΔF is defined as the fluorescent temporal signal normalised by the initial signal strength. (b) [Ca2+]i traces from three different cultured astrocytes at 28°C and 20°C. Note that the transient fluorescence changes have a smaller spike width and occur less frequently at higher temperature. The recordings were separated by 5 min. (c) Average number of cells which exhibited a [Ca2+]i spike (average from four independent experiments). The number of responsive cells at 20°C was set to 100% and compared to the number of responsive cells at different temperatures or to the recovery at 20°C (bars indicate s.e.; P∗<0.05). ((d), (e)). Average spike width of the fluorescent signals of all experiments obtained from cultured astrocytes (d) and in slices (e) determined at different temperatures as indicated (error bars indicate s.e.). (f) From the data the duration of the spike width is plotted as a function of temperature for astrocytes from culture (black dots) and brain slices (circles). The data points were fitted by f(T)=620s·exp⁡(-0.109T/∘C) (P∗∗<0.005).
(a)
(b)
(c)
(d)
(e)
(f)The observed reduction in activity, starting at 28°C, was significantly smaller compared to recordings with the higher temperature of 30°C (Figure1), indicating that between 28°C and 30°C is a threshold-like temperature with regard to [Ca2+]i spike activity. The rather instantaneous response to temperature changes and the reversibility of the change in the activity pattern suggests that temperature is directly modifying the characteristic property of some molecular properties of the Ca2+ tool kit and not acting primarily on transcriptional regulation.


#### 2.1.2. Temperature Changes Modulate Signal Form of [Ca2+]i Spike

Careful inspection of the time courses reveals also a strong temperature dependent on the duration and signal form of the Ca2+ transient in both cultured astrocytes and astrocytes in acute brain slices as shown in Figure 1. To quantify the effect on single [Ca2+]i-transients, we measured the spike width (SW) of individual Ca2+ transients by the full width half maximum (FWHM) method. Despite the cell specific heterogeneity, spike durations of the Ca2+ transients of astrocytes in acute brain slices decreased significantly when the temperature was raised from 24°C to 37°C. The decrease from 45 ± 17 s to 13 ± 4 s was fully reversible, since lowering the temperature back to 24°C increased the spike duration again to 45 ± 18 s (10 slices, 42 cells at T24°C  70 spikes; T37°C 53 spikes; T24°C 71 spikes). We additionally determined the spike duration at 31°C which averaged at 22 ± 5 s (9 slices, 29 cells at T31°C 43 spikes) (Figure 1).Within the analogous analysis of cultured astrocytes, we found that the spike duration changed from an average of 63 ± 18 s at 20°C to 21 ± 5 s at 30°C. This effect of temperature on spike duration was again reversible as spike duration increased again to 72 ± 20 s at 20°C after lowering the temperature from 30°C (4 independent experiments; atT30°C 15 cells, 26 spikes; at T20∘C 30 cells, 56 spikes) (Figure 1). Within the other experimental setup, the average spike duration was 27 ± 6 s for 28°C as compared to 75 ± 16 s at 20°C (4 experiments; 30 cells at T28°C 103 spikes; T20°C 105 spikes).The similarity of the signal form and particularly of the quantified spike widths of cells obtained from cultures and brain slices around 30°C points to temperature as the main reason for the typically reported differences in the signals. To analyze this behavior further we plotted the temperature dependence of the spike duration by pooling data from the two different experimental systems. The resulting data points can be fitted nicely by the exponential functionf(T)=ae-bT with a=620 s and b=0.109/°C indicating a functional relation between spike width and temperature.Within the experimental comparison study we also investigated the role of extracellular Ca2+ influx and intracellular Ca2+ release from intracellular stores. While Ca2+ free buffer has reduced the total number of responsive cells only slightly to 78 ± 34% (n=13 slices, averaged cell number/slice = 31 ± 15), blocking SERCA pumps by loading 1 μM thapsigargin for 15 min has abolished Ca2+ signals in more than 90% of cells at 31°C. Overall, these findings indicate a rapid change of the intracellular signaling mechanism. In particular, the astonishing reversibility of the Ca2+ activity and its signal forms when resetting the temperature to its previous value after tens of minutes suggests a molecular activity dependency on temperature.


### 2.2. Modeling the [Ca2+]i Temperature Dependence by Modulated SERCA Activity

As previously speculated, the experimentally reported temperature dependent activity of the SERCA pumps [26, 50] may explain the observed changes in [Ca2+]i signal forms. These studies report an increase of the pump strength with increasing temperature from 5°C to 38°C. Thus, the Ca2+ uptake at 5°C amounts only 5% of the uptake at 38°C. This can have a large effect on CICR that represents the spike generating mechanism because the IP3Rs are not homogeneously distributed within the cell but form channel clusters that are typically separated by 1 μm or more. This separation may lead together with Ca2+ binding buffers and increasing SERCA activity to steeper gradients close to open channel clusters and subsequently to a weaker coupling between release sites and hence to a reduced CICR.The data shown in Figure1 support this perspective. The individual spike form is influenced by the stronger pumping as demonstrated experimentally by the decreased spike width for larger temperatures. Moreover, the spike widths of the two different cell types could be fitted nicely by one exponential function of the temperature. To test this SERCA based hypothesis we apply here our developed physiological multiscale model that is able to generate the whole spectrum of known [Ca2+]i signals and was used to investigate the experimentally quantified stochastic nature of [Ca2+]i spiking including the functional robustness of the σ-Tav relation [27].


#### 2.2.1. Multiscale Modeling Approach

Based on the observations that Ca2+ spiking is stochastic and obeys a nontrivial signaling signature, our multiscale model considers the molecular dynamics of individual IP3Rs and Ca2+ diffusion as well as buffer dynamics [27]. Our modeling strategy for this coupled system is based on separation of the two involved length scales. On the microscopic scale we use a detailed model for the channels [25]. Those will be open and close in a stochastic manner due to stochastic binding of the signaling molecules IP3 and Ca2+. The stochastic channel dynamics acts as a noisy source term that is driving the three-dimensional reaction diffusion system (RDS) representing the cytosol. Thereby the model also takes reactions of Ca2+ with buffer molecules like BAPTA or calmodulin and the SERCA pumps into account. From this perspective the cellular dynamics is given by an in general nonlinear reaction diffusion system that considers the cytosolic and lumenal dynamics by (1a)∂[Ca2+]∂t=Dc∇2[Ca2+]+σ~l+Pc(r,t)·([E]-[Ca2+])-P~pCa2+2kd2+Ca2+2+∑j‍kj-CaBj-kj+BjCa2+,(1b)∂[E]∂t=DE∇2[E]-γCa2+2kd2+Ca2+2σ~l+Pcr,tE-Ca2+-P~pCa2+2kd2+Ca2+2+∑k‍kk-EGk-kk+GkE,(1c)∂[Bj]∂t=DBj∇2Bj-kj+BjCa2++kj-CaBj,(1d)∂[Gk]∂t=DGk∇2[Gk]-kk+[Gk][E]+kk-[EGk],where the first two equations describe the free Ca2+ concentration in the cytosol and the ER, respectively, and the other two equations correspond to a variety of cytosolic and lumenal buffers [Bj] and [Gk]. The first terms on the right hand side characterize diffusion with the individual diffusion coefficients Di. The following terms in (1a) and (1b) describe the leak flux (with permeability σl) and release through channels (Pc). The space dependence of Pc(rc,t) reflects the location of channel clusters and the time dependent stochastic opening and closing of the channels. This time dependent source term is governed by stochastic simulation of the channel states. γ is the volume ratio Vcyt/VER. The third term results from Ca2+ uptake into the ER from the cytosol by SERCA pumps with a flux strength Pp and the SERCA dissociation constant kd. The remaining terms in (1a)–(1d) describe the reactions of Ca2+ with cytosolic and lumenal buffers [Bk] and [Gj], respectively, where k+ denotes the capture rate and k- is the dissociation rate.The cellular Ca2+ dynamics based on the CICR mechanism results from the local properties of the IP3R channels and their spatial coupling by Ca2+ diffusion that is influenced by buffer reactions and SERCA pump activity. Thus, the cellular Ca2+ signaling is determined by the channel states and corresponding release fluxes that shape the signals in dependence on the properties of the cellular diffusion. In turn, the resulting Ca2+ concentration profiles influence the channel opening behavior due to its dependence on the local Ca2+ concentration.Standard simulation techniques for reaction diffusion systems like finite difference or finite element methods have to calculate the concentration profiles for the complete numerical mesh for convergence criteria. The large concentration gradients typically observed close to open channels require a fine spatial discretization that leads in particular in three spatial dimensions to large computational costs [51]. To circumvent this restriction and enable simulations of biological relevant time courses, we developed an analytic solution of the reaction diffusion equations in terms of coupled Green’s function that allows for calculating [Ca2+]i at specific spatial positions only such as channel cluster locations. These local concentrations subsequently determine potential changes in the stochastic channel dynamics leading to modified source terms in the reaction diffusion system.In order to solve the reaction diffusion system analytically we restrict the system to one mobile and one immobile buffer, introduce dimensionless variables, and linearize the corresponding equations [27]. To obtain an analytic solution of the resulting reaction diffusion system we have to determine the geometry and boundary conditions. As a natural choice we use a spherical geometry and, motivated by the observation that astrocytic Ca2+ signals do not change significantly in Ca2+ free medium [24], no-flux boundary condition that reflects the properties of the plasma membrane for our in silico cell. With these specifications we could develop the analytic solution by coupled Green’s functions [52, 53] describing Ca2+ and buffer dynamics where open channels correspond to source terms. For mathematical details of the implemented Green’s cell algorithm (GCA) see [27] where we have shown how microscopic channel fluctuations are transmitted onto the level of the cell by diffusion thereby matching the experimental findings on cell variability and buffer influence.For the IP3R description within the GCA, we use the DeYoung-Keizer model (DKM) that assumes three binding sites for each of the four subunits per channel [54]. Each of the three binding sites can be free or occupied leading to 23 different states Xijk and 12 possible transitions in dependence on the Ca2+ concentration C and IP3 concentration I that can be visualized on a cube as shown in Figure 2(a). A subunit is active in the state X110 only where the first index equals one if IP3 is bound, the second index corresponds to Ca2+ binding to an activating site, and the last index describes the occupation of the dominant inhibiting Ca2+ binding site. The activating binding sites for Ca2+ have a higher affinity compared to the dominant inhibiting site. This is a minimal choice to generate the nonlinearity in Ca2+, which is the source of the bell shaped open probability (Figure 2(b)) and the basis of CICR [25].Figure 2
CICR is based on the nonlinear open probability of IP3R. (a) Scheme of the DeYoung-Keizer model for a single subunit. A subunit is active, if IP3 (I) is bound and Ca2+ (C) is only bound to the activating site, that is, in state X110. A channel opens if at least 3 of its 4 subunits are active. Open IP3Rs release Ca2+ that can diffuse to other IP3R clusters and activate the channels there as well, eventually leading to cell wide Ca2+ release and generating a [Ca2+]i spike. From the cytosol Ca2+ is pumped back into the ER by SERCA that decreases diffusion. (b) From the Markov chain model of a single subunit the non-Markovian stationary open probability Po of a channel can be calculated in dependence on Ca2+ and IP3 concentration [25].
(a)
(b)


#### 2.2.2. Modeling the Impact of SERCA Pump Strength on Ca2+ Dynamics

Here we use our validated GCA implementation to test the hypothesis on the temperature dependent SERCA activity as a main reason for the observed differences in[Ca2+]i spiking. For the simulations, we use physiological cellular diffusion parameters listed in Table 1. Our in silico cell has 31 randomly distributed IP3R clusters each of them having a random number of channels between 2 and 16, yielding 258 channels in total. The clusters are separated by at least 1.5 μm and only located within a sphere with 80% of the cell radius R to guarantee convergence of Green’s function representation. To test the SERCA hypothesis, this prototype of a cell was simulated with different pump strengths Pp.Table 1
Physiologic parameters of the in silico cell used in simulations. Note that the GCA uses a linearized reaction diffusion system leading to rates forσl and Pp.
Physiological cell parameter values for GCA simulations
R
10μm
Cell radius
γ
5
Volume ratio between cytosol and ER
D
Ca
220μm2/s
Diffusion coefficient of cytosolicCa2+
D
E
70μm2/s
Diffusion coefficient of lumenalCa2+
D
B
95μm2/s
Diffusion coefficient of mobile buffer
[
C
a
2
+
]
0
55 nM
CytosolicCa2+ base level
[
E
]
0
800μM
Endoplasmic reticulumCa2+ base level
[
I
P
3
]
80μM
IP
3 concentration
[
B
]
T
52μM
Total mobile buffer concentration (EGTA)
k
B
+
0.3 (μMs)−1
On rate of the mobile buffer
k
B
-
1.5 s−1
Dissociation rate of the mobile buffer
[
B
i
]
T
30μM
Total immobile buffer concentration
k
B
i
+
600 (μMs)−1
On rate of the immobile buffer
k
B
i
-
100 s−1
Dissociation rate of the immobile buffer
σ
l
≈0.01 s−1
Leak flux constant given by variedPp, [Ca2+]0 and [E]0Figure3 shows the cell wide dynamics, where panel (a) exhibits the number of open channels and panel (b) shows the resulting cytosolic Ca2+ concentration for pump strengths Pp varying between 22 s−1 and 200 s−1 what is around the estimated physiological reference value of 86 s−1 [55]. A first conspicuous property shown in panel (a) is that the amplitudes of open channels during a spike are rather similar and vary between 25 and 30, rather independently of the pump strengths Pp. Nevertheless, Pp has a huge influence on the time course.Figure 3
Dependence of Ca2+ signals on the pump strength. The temperature dependent activity of SERCA leads to higher pump rates Pp for higher temperatures. This influences the Ca2+ signals as shown by the number of open channels (a) and the corresponding Ca2+ concentrations (b). For the two smallest strengths the simulations ended before the shown real time of 800 s caused by run time restrictions on the compute cluster.For low pump strengths, Ca2+ is removed too slowly from the cytosol for a coordinated signal. During the initial activation from base level the amplitude of open channels becomes maximal and is limited by inhibition according to the channel properties explained in Figure 2. When a channel recovers from inhibition due to the stochastic unbinding of Ca2+ from a dominantly inhibiting site of a subunit, the cytosolic Ca2+ concentration is still in the Ca2+ concentration range where the higher affinity for the Ca2+ activating binding site leads to a large open probability shown in Figure 2(b). Subsequently, the channel will open again, increase the local Ca2+ concentration, and reinduce channel inhibition. Repeating this scenario, the channels will mainly switch between activation and inhibition back and forth where a rather constant proportion of up to 50% of subunits are inhibited until randomly no reactivation occurs and the cell is reaching its resting state again. This permanent switching between activation and inhibition leads to the observed plateaus for the two lowest pump strengths.This behavior is also reflected by the resulting cytosolic Ca2+ concentrations. For Pp=22 s−1, the cell exhibits a plateau response with some superimposed [Ca2+]i spiking, which are often observed in experiments with astrocytes from acute brain slices [24, 56–58]. For Pp=32 s−1, the shorter plateaus of open channels are translated into a more oscillation-like behavior of the cytosolic Ca2+ concentration exhibiting rather diverse spike amplitudes. The concentration is in a similar range as the cell with the lower pump strength, indicating a regime where local channel inhibition and spatial coupling determined by pumps have an equal influence on the cellular dynamics.For increasingPp, the cells exhibit an even more spiking like behavior since long lasting channel activity is absent, and therefore [Ca2+]i shows more pronounced peaks. This can be again explained by the interplay of inhibition and Ca2+ removal. After a spike has occurred and most channels are inhibited, the SERCAs decrease the cytosolic Ca2+ concentration sufficiently fast leading to base level [Ca2+]i when channels recover from inhibition. Subsequently, the cell relaxes to its resting level, from which it can be activated again by CICR, that is, by triggering opening of a channel cluster. For very large pump strengths, the spikes of the number of open channels have still a similar height to those for smaller Pp, but their widths decrease drastically. This induces low and slim peaks of the cytosolic Ca2+ concentration.So far, the simulations seem to support the hypothesis of the decreased SERCA activity leading to an increased[Ca2+]i activity at lower temperatures. In terms of the experimental findings, a temperature increase in experiments and observed changes in the Ca2+ signals correspond to an increase of the pump rate from the region from 42 s−1 to 62 s−1 to 152  s−1 or higher.For a further detailed analysis, we determine the spike widths of the simulated cytosolic Ca2+ peaks as done for the experimental data in Figure 1. Figure 4(a) exhibits the dependence of the spike width on the pump strength Pp. The analysis reveals an exponential dependence on the pump rate Pp as shown by the line defined by SW=61s·e-0.012s·Pp. At second glance we observe that the spike width becomes most regular for Pp=62 s−1 as depicted by the small error bar that exhibits the smallest relative variation even compared to the smaller error bars for large Pp. This might indicate an optimal dynamical regime.Figure 4
Spike width and amplitude in dependence on the pump strength for the Ca2+ signal shown in Figure 3. (a) The spike width exhibits an exponential dependence on Pp. This is in good agreement with the experimental findings in Figure 1. (b) Also the amplitude dependence on Pp obeys an exponential relation (error bars denote s.e.).
(a)
(b)An analogue behavior can be found for the average amplitudes in Figure4(b). The dependence on Pp can be fitted by the exponential relation 2.9·e-0.011s·Pp, and again the mean amplitude for Pp=62 s−1 exhibits the smallest relative variations. The amplitude behavior obtained from simulations has a similar trend like those seen in the experiments but contradicts former theoretical results where the amplitude stayed constant [59].To analyze the potential optimal regime, we determine the oscillation characteristicsTav and σ. Indeed, Tav exhibits a minimum of 42 s for a pump rate of 62 s−1 as shown in Figure 5(a). Panel (b) displays the dependence of σ on Tav showing a slope close to one as in earlier experiments and simulations and indicates a deterministic time of 20 s. The spread of the data points is exclusively caused by the different pump strengths. Similar simulations with smaller buffer concentrations have shown an analogue behavior with a shift of the minimum of Tav to larger Pp. This is in line with our consideration, since a lower buffer concentration leads to an increased spatial coupling that can be compensated by higher pump rates.Figure 5
Dependence of the mean periodTav on the pump strength Pp. (a) Tav exhibits a nonlinear dependence on the pump strength with a minimal mean period of 42 s for Pp=62 s−1. This might illuminate a possible control mechanism cells can use to tune their oscillations. (b) The σ-Tav relation exhibits a population slope close to one and a deterministic time of approximately 20 s. (c) The height of the local Ca2+ concentration Ca2+ close to a cluster is hardly affected by Pp as demonstrated by the representative example of the same cluster in different simulations using different pump rates. Hence, the differences in the cytosolic Ca2+ oscillations in Figure 3(b) are spatially induced.
(a)
(b)
(c)For an additional test of our conception based on the spatial coupling of release sites, we analyze the local Ca2+ concentrations at channel clusters in Figure 5(c). For all pump strengths the concentration peaks are in the range of 100 μM independently of Pp. Hence, the different cell wide cytosolic Ca2+ signals shown in Figure 3(b) are not caused by different local properties but are caused by the spatial coupling of release sites. Thus, the results of the temperature experiments are consistent with the hypothesized spatial control mechanism.


#### 2.2.3. Rate Equation Based Modeling

To further investigate the role of the spatial aspect in Ca2+ signaling, we model the effect of SERCA activity on the Ca2+ profile by a traditional rate equation model implemented by ordinary differential equations (ODEs). These types of models do not consider any spatial aspects but treat the cell as a well-stirred homogeneous entity and were extensively used in the establishment of the Ca2+ signaling mechanism [60, 61]. A variety of models have been proposed and applied to diverse aspects of Ca2+ signaling where the common mechanism is Ca2+ cycling between the ER and the cytosol. Here we use a model based on previous studies which reflects the main Ca2+ fluxes according to Figure 2 as well as Ca2+ buffering by mitochondria. The governing rate equations for the cytosolic and ER Ca2+ concentration read(2a)dCacdt=fcytJext+ΔJmito+Jrelease+Jleak-Jserca-Jpmca,(2b)dCaerdt=ferγ-1Jserca-Jrelease-Jleak,where fcyt and fer are the fractions of free Ca2+ within the cytosol and the ER [62], γ is again the volume ratio between cytosol and ER, Δ denotes the mitochondria to cytosol volume ratio, and the corresponding fluxes J reflect typical physiological characteristics. The 2 most important fluxes are the Ca2+ release from the ER by IP3R (Jrelease) and the SERCA mediated Ca2+ removal from the cytosol into the ER (Jserca) given by (3a)Jrelease=PrIP32p12+IP32Cac2p22+Cac2·p34p44+Cac4Caer-Cac,(3b)Jserca=PpCac2p52+Cac2,respectively. The dominant flux for the [Ca2+]i spiking is the nonlinear release flux inducing the CICR by the bell shaped dependence on the Ca2+ concentration [63] here modeled in analogy to previous approaches [64]. The form of the complementary SERCA mediated Ca2+ flux into the ER (Jserca) reflects the fact that two Ca2+ ions are pumped into the ER by the use of one ATP molecule where the ATP concentration is assumed to be constant. The remaining fluxes are given in Table 2 where Jext and Jmito describe Ca2+ flux from extracellular space and mitochondria, respectively. Jleak corresponds to the Ca2+ leakage from ER and Jpmca is the flux into the extracellular space through plasma membrane Ca2+ ATPases. The corresponding parameters pi are also listed in Table 2.Table 2
Physiologic parameters of the in silico cell used for the rate equation simulations.
Additional flux definitions of the full rate equation model
Sodium-Ca2+ exchanger mediated flux
J
NaCa
=
p
5
C
a
m
C
a
m
+
p
6
p
7
2
p
8
2
+
C
a
c
2
Mitochondrial unitransporter mediated flux
J
uni
=
p
9
C
a
c
2
p
10
2
+
C
a
c
2
Effective mitochondrial flux
J
mito
=
J
NaCa
-
J
uni
Leak flux through ER membrane
J
leak
=
p
11
C
a
er
-
C
a
c
Influx from external space
J
ext
=
J
ext,bas
-
(
p
12
Ca
er
+
p
13
Ca
c
)
PMCA mediated flux
J
pmca
=
p
14
C
a
c
2
p
15
2
+
C
a
c
2
Parameter values for the rate equation calcium model
P
r = 0.03 ms−1
p
1 = 0.5 μM
p
2 = 0.9 μM
p
3 = 1.3 μM
p
4 = 0.2 μM
P
p = 0.06 μM ms−1
p
5 = 10 μM ms−1
p
6 = 0.6 μM
p
7= 1 μM
p
8 = 10 μM
p
9 = 0.2 μM ms−1
p
10 = 1 μM
p
11 = 10-3 ms−1
p
12 = 3·10-6 ms−1
p
13 = 15·10-6 ms−1
p
14 = 0.01 μM ms−1
p
15 = 0.7 μM
IP
3 = 0.4 μM
f
cyt
=
0.015
f
er = 0.5
J
ext,bas = 3.5·10-3 μM ms−1
Δ
=
0.07
γ
-
1
=
30To investigate the influence of SERCA pump strength on Ca2+ dynamics, we simulate the ODEs for a wide range of pump strength Pp and three different values of the Ca2+ release flux (Pr) from IP3R channels and analyze the spike characteristics (Figure 6). This parameter scan shows that the Ca2+ period within the rate equation model exhibits a similar behavior as the spatial GCA simulations in terms of a minimal period in dependence on the pump strength (Figure 6(a)). While this finding is not in clear favor for the perspective of spatial coupling of release sites as the [Ca2+]i spike generating mechanism, the results obtained with the rate equation model given by (2a) and (2b) about the spike width and amplitude exhibit a different behavior than the experimental data.Figure 6
Cell characteristics in the rate equation model. Dependence of the period (a), amplitude (b), and spike width (c) on the SERCA pump strengthPp for three different Ca2+ release rates Pr from IP3 channels.
(a)
(b)
(c)In the ODE model a higher SERCA pump rate causes more Ca2+ within the ER and also a higher Ca2+ release into the cytosol because IP3R flux is controlled to a large proportion by the concentration difference and thus by the ER concentration. A high SERCA pump rate leads to a larger spike width as observed in Figure 6(b). These higher fluxes in and out of the ER also yield higher spike amplitudes (Figure 6(c)). Both observations are in contrast to the experimental observations and the GCA simulations, potentially indicating the fact that a spatial distribution of channels is a necessity for the experimentally observed SERCA pump effects.


## 2.1. Experiments Show Strong Temperature Dependency in Cortical and Cultured Astrocytes

In the experimental study [24] we used cortical slices obtained from P 10–P 12 mouse brains and cultured astrocytes. The temperature was either changed abruptly or smoothly either by changing the heating temperature of the water bath for the perfusion buffer or by switching from one perfusion buffer to another, respectively. Thereby the temperature was measured within the recording chamber close to the sample cover slip.


### 2.1.1. Higher Temperatures Reversibly Suppress [Ca2+]i Spiking Instantaneously

Starting from the standard protocols, we first measured Ca2+ activity at room temperature and separately at more physiological temperatures with Fluo4-AM loaded to the cells before (more detailed experimental description can be found in [24]). In cortical slices obtained from P 10–P 12 mouse brains measured at room temperature, astrocytes showed robust spontaneous signaling activity defined by long durations and large amplitudes of [Ca2+]i spikes. In contrast, when slices were measured at physiological temperatures, only very few cells showed spontaneous activity.To compare this behavior with primary astrocytes cultured from mouse neocortex, we measured Ca2+ activity at room temperature (20°C) and then abruptly changed the perfusion buffer to one of 28°C or of 30°C in another setup. For these cultured astrocytes we observed a similar significant decrease of cells exhibiting Ca2+ signals already within the first minute after the temperature shift (Figure 1). Out of the active cells under 20°C (72 ± 24% of all cells in 4 independent experiments) more than a quarter of cells have not shown any [Ca2+]i spikes under 28°C and under 30°C only 20% exhibited at least a single [Ca2+]i spike. Remarkably, this dramatic decrease was fully reversible when perfusing the cells again with buffer at room temperature. Analogous results about the reversibility were also found in experiments with acute brain slices and the reversibility was independent of the experimental order of temperature changes.Figure 1
Summary of experimental observations on temperature dependent[Ca2+]i activity [24]. (a) Fluorescence recordings from three individual astrocytes in three different experiments. The temperature was changed between 20°C and 30°C by using heated and nonheated perfusion buffer (indicated by the bar). Fluorescent amplitude ΔF is defined as the fluorescent temporal signal normalised by the initial signal strength. (b) [Ca2+]i traces from three different cultured astrocytes at 28°C and 20°C. Note that the transient fluorescence changes have a smaller spike width and occur less frequently at higher temperature. The recordings were separated by 5 min. (c) Average number of cells which exhibited a [Ca2+]i spike (average from four independent experiments). The number of responsive cells at 20°C was set to 100% and compared to the number of responsive cells at different temperatures or to the recovery at 20°C (bars indicate s.e.; P∗<0.05). ((d), (e)). Average spike width of the fluorescent signals of all experiments obtained from cultured astrocytes (d) and in slices (e) determined at different temperatures as indicated (error bars indicate s.e.). (f) From the data the duration of the spike width is plotted as a function of temperature for astrocytes from culture (black dots) and brain slices (circles). The data points were fitted by f(T)=620s·exp⁡(-0.109T/∘C) (P∗∗<0.005).
(a)
(b)
(c)
(d)
(e)
(f)The observed reduction in activity, starting at 28°C, was significantly smaller compared to recordings with the higher temperature of 30°C (Figure1), indicating that between 28°C and 30°C is a threshold-like temperature with regard to [Ca2+]i spike activity. The rather instantaneous response to temperature changes and the reversibility of the change in the activity pattern suggests that temperature is directly modifying the characteristic property of some molecular properties of the Ca2+ tool kit and not acting primarily on transcriptional regulation.


### 2.1.2. Temperature Changes Modulate Signal Form of [Ca2+]i Spike

Careful inspection of the time courses reveals also a strong temperature dependent on the duration and signal form of the Ca2+ transient in both cultured astrocytes and astrocytes in acute brain slices as shown in Figure 1. To quantify the effect on single [Ca2+]i-transients, we measured the spike width (SW) of individual Ca2+ transients by the full width half maximum (FWHM) method. Despite the cell specific heterogeneity, spike durations of the Ca2+ transients of astrocytes in acute brain slices decreased significantly when the temperature was raised from 24°C to 37°C. The decrease from 45 ± 17 s to 13 ± 4 s was fully reversible, since lowering the temperature back to 24°C increased the spike duration again to 45 ± 18 s (10 slices, 42 cells at T24°C  70 spikes; T37°C 53 spikes; T24°C 71 spikes). We additionally determined the spike duration at 31°C which averaged at 22 ± 5 s (9 slices, 29 cells at T31°C 43 spikes) (Figure 1).Within the analogous analysis of cultured astrocytes, we found that the spike duration changed from an average of 63 ± 18 s at 20°C to 21 ± 5 s at 30°C. This effect of temperature on spike duration was again reversible as spike duration increased again to 72 ± 20 s at 20°C after lowering the temperature from 30°C (4 independent experiments; atT30°C 15 cells, 26 spikes; at T20∘C 30 cells, 56 spikes) (Figure 1). Within the other experimental setup, the average spike duration was 27 ± 6 s for 28°C as compared to 75 ± 16 s at 20°C (4 experiments; 30 cells at T28°C 103 spikes; T20°C 105 spikes).The similarity of the signal form and particularly of the quantified spike widths of cells obtained from cultures and brain slices around 30°C points to temperature as the main reason for the typically reported differences in the signals. To analyze this behavior further we plotted the temperature dependence of the spike duration by pooling data from the two different experimental systems. The resulting data points can be fitted nicely by the exponential functionf(T)=ae-bT with a=620 s and b=0.109/°C indicating a functional relation between spike width and temperature.Within the experimental comparison study we also investigated the role of extracellular Ca2+ influx and intracellular Ca2+ release from intracellular stores. While Ca2+ free buffer has reduced the total number of responsive cells only slightly to 78 ± 34% (n=13 slices, averaged cell number/slice = 31 ± 15), blocking SERCA pumps by loading 1 μM thapsigargin for 15 min has abolished Ca2+ signals in more than 90% of cells at 31°C. Overall, these findings indicate a rapid change of the intracellular signaling mechanism. In particular, the astonishing reversibility of the Ca2+ activity and its signal forms when resetting the temperature to its previous value after tens of minutes suggests a molecular activity dependency on temperature.


## 2.1.1. Higher Temperatures Reversibly Suppress [Ca2+]i Spiking Instantaneously

Starting from the standard protocols, we first measured Ca2+ activity at room temperature and separately at more physiological temperatures with Fluo4-AM loaded to the cells before (more detailed experimental description can be found in [24]). In cortical slices obtained from P 10–P 12 mouse brains measured at room temperature, astrocytes showed robust spontaneous signaling activity defined by long durations and large amplitudes of [Ca2+]i spikes. In contrast, when slices were measured at physiological temperatures, only very few cells showed spontaneous activity.To compare this behavior with primary astrocytes cultured from mouse neocortex, we measured Ca2+ activity at room temperature (20°C) and then abruptly changed the perfusion buffer to one of 28°C or of 30°C in another setup. For these cultured astrocytes we observed a similar significant decrease of cells exhibiting Ca2+ signals already within the first minute after the temperature shift (Figure 1). Out of the active cells under 20°C (72 ± 24% of all cells in 4 independent experiments) more than a quarter of cells have not shown any [Ca2+]i spikes under 28°C and under 30°C only 20% exhibited at least a single [Ca2+]i spike. Remarkably, this dramatic decrease was fully reversible when perfusing the cells again with buffer at room temperature. Analogous results about the reversibility were also found in experiments with acute brain slices and the reversibility was independent of the experimental order of temperature changes.Figure 1
Summary of experimental observations on temperature dependent[Ca2+]i activity [24]. (a) Fluorescence recordings from three individual astrocytes in three different experiments. The temperature was changed between 20°C and 30°C by using heated and nonheated perfusion buffer (indicated by the bar). Fluorescent amplitude ΔF is defined as the fluorescent temporal signal normalised by the initial signal strength. (b) [Ca2+]i traces from three different cultured astrocytes at 28°C and 20°C. Note that the transient fluorescence changes have a smaller spike width and occur less frequently at higher temperature. The recordings were separated by 5 min. (c) Average number of cells which exhibited a [Ca2+]i spike (average from four independent experiments). The number of responsive cells at 20°C was set to 100% and compared to the number of responsive cells at different temperatures or to the recovery at 20°C (bars indicate s.e.; P∗<0.05). ((d), (e)). Average spike width of the fluorescent signals of all experiments obtained from cultured astrocytes (d) and in slices (e) determined at different temperatures as indicated (error bars indicate s.e.). (f) From the data the duration of the spike width is plotted as a function of temperature for astrocytes from culture (black dots) and brain slices (circles). The data points were fitted by f(T)=620s·exp⁡(-0.109T/∘C) (P∗∗<0.005).
(a)
(b)
(c)
(d)
(e)
(f)The observed reduction in activity, starting at 28°C, was significantly smaller compared to recordings with the higher temperature of 30°C (Figure1), indicating that between 28°C and 30°C is a threshold-like temperature with regard to [Ca2+]i spike activity. The rather instantaneous response to temperature changes and the reversibility of the change in the activity pattern suggests that temperature is directly modifying the characteristic property of some molecular properties of the Ca2+ tool kit and not acting primarily on transcriptional regulation.


## 2.1.2. Temperature Changes Modulate Signal Form of [Ca2+]i Spike

Careful inspection of the time courses reveals also a strong temperature dependent on the duration and signal form of the Ca2+ transient in both cultured astrocytes and astrocytes in acute brain slices as shown in Figure 1. To quantify the effect on single [Ca2+]i-transients, we measured the spike width (SW) of individual Ca2+ transients by the full width half maximum (FWHM) method. Despite the cell specific heterogeneity, spike durations of the Ca2+ transients of astrocytes in acute brain slices decreased significantly when the temperature was raised from 24°C to 37°C. The decrease from 45 ± 17 s to 13 ± 4 s was fully reversible, since lowering the temperature back to 24°C increased the spike duration again to 45 ± 18 s (10 slices, 42 cells at T24°C  70 spikes; T37°C 53 spikes; T24°C 71 spikes). We additionally determined the spike duration at 31°C which averaged at 22 ± 5 s (9 slices, 29 cells at T31°C 43 spikes) (Figure 1).Within the analogous analysis of cultured astrocytes, we found that the spike duration changed from an average of 63 ± 18 s at 20°C to 21 ± 5 s at 30°C. This effect of temperature on spike duration was again reversible as spike duration increased again to 72 ± 20 s at 20°C after lowering the temperature from 30°C (4 independent experiments; atT30°C 15 cells, 26 spikes; at T20∘C 30 cells, 56 spikes) (Figure 1). Within the other experimental setup, the average spike duration was 27 ± 6 s for 28°C as compared to 75 ± 16 s at 20°C (4 experiments; 30 cells at T28°C 103 spikes; T20°C 105 spikes).The similarity of the signal form and particularly of the quantified spike widths of cells obtained from cultures and brain slices around 30°C points to temperature as the main reason for the typically reported differences in the signals. To analyze this behavior further we plotted the temperature dependence of the spike duration by pooling data from the two different experimental systems. The resulting data points can be fitted nicely by the exponential functionf(T)=ae-bT with a=620 s and b=0.109/°C indicating a functional relation between spike width and temperature.Within the experimental comparison study we also investigated the role of extracellular Ca2+ influx and intracellular Ca2+ release from intracellular stores. While Ca2+ free buffer has reduced the total number of responsive cells only slightly to 78 ± 34% (n=13 slices, averaged cell number/slice = 31 ± 15), blocking SERCA pumps by loading 1 μM thapsigargin for 15 min has abolished Ca2+ signals in more than 90% of cells at 31°C. Overall, these findings indicate a rapid change of the intracellular signaling mechanism. In particular, the astonishing reversibility of the Ca2+ activity and its signal forms when resetting the temperature to its previous value after tens of minutes suggests a molecular activity dependency on temperature.


## 2.2. Modeling the [Ca2+]i Temperature Dependence by Modulated SERCA Activity

As previously speculated, the experimentally reported temperature dependent activity of the SERCA pumps [26, 50] may explain the observed changes in [Ca2+]i signal forms. These studies report an increase of the pump strength with increasing temperature from 5°C to 38°C. Thus, the Ca2+ uptake at 5°C amounts only 5% of the uptake at 38°C. This can have a large effect on CICR that represents the spike generating mechanism because the IP3Rs are not homogeneously distributed within the cell but form channel clusters that are typically separated by 1 μm or more. This separation may lead together with Ca2+ binding buffers and increasing SERCA activity to steeper gradients close to open channel clusters and subsequently to a weaker coupling between release sites and hence to a reduced CICR.The data shown in Figure1 support this perspective. The individual spike form is influenced by the stronger pumping as demonstrated experimentally by the decreased spike width for larger temperatures. Moreover, the spike widths of the two different cell types could be fitted nicely by one exponential function of the temperature. To test this SERCA based hypothesis we apply here our developed physiological multiscale model that is able to generate the whole spectrum of known [Ca2+]i signals and was used to investigate the experimentally quantified stochastic nature of [Ca2+]i spiking including the functional robustness of the σ-Tav relation [27].


### 2.2.1. Multiscale Modeling Approach

Based on the observations that Ca2+ spiking is stochastic and obeys a nontrivial signaling signature, our multiscale model considers the molecular dynamics of individual IP3Rs and Ca2+ diffusion as well as buffer dynamics [27]. Our modeling strategy for this coupled system is based on separation of the two involved length scales. On the microscopic scale we use a detailed model for the channels [25]. Those will be open and close in a stochastic manner due to stochastic binding of the signaling molecules IP3 and Ca2+. The stochastic channel dynamics acts as a noisy source term that is driving the three-dimensional reaction diffusion system (RDS) representing the cytosol. Thereby the model also takes reactions of Ca2+ with buffer molecules like BAPTA or calmodulin and the SERCA pumps into account. From this perspective the cellular dynamics is given by an in general nonlinear reaction diffusion system that considers the cytosolic and lumenal dynamics by (1a)∂[Ca2+]∂t=Dc∇2[Ca2+]+σ~l+Pc(r,t)·([E]-[Ca2+])-P~pCa2+2kd2+Ca2+2+∑j‍kj-CaBj-kj+BjCa2+,(1b)∂[E]∂t=DE∇2[E]-γCa2+2kd2+Ca2+2σ~l+Pcr,tE-Ca2+-P~pCa2+2kd2+Ca2+2+∑k‍kk-EGk-kk+GkE,(1c)∂[Bj]∂t=DBj∇2Bj-kj+BjCa2++kj-CaBj,(1d)∂[Gk]∂t=DGk∇2[Gk]-kk+[Gk][E]+kk-[EGk],where the first two equations describe the free Ca2+ concentration in the cytosol and the ER, respectively, and the other two equations correspond to a variety of cytosolic and lumenal buffers [Bj] and [Gk]. The first terms on the right hand side characterize diffusion with the individual diffusion coefficients Di. The following terms in (1a) and (1b) describe the leak flux (with permeability σl) and release through channels (Pc). The space dependence of Pc(rc,t) reflects the location of channel clusters and the time dependent stochastic opening and closing of the channels. This time dependent source term is governed by stochastic simulation of the channel states. γ is the volume ratio Vcyt/VER. The third term results from Ca2+ uptake into the ER from the cytosol by SERCA pumps with a flux strength Pp and the SERCA dissociation constant kd. The remaining terms in (1a)–(1d) describe the reactions of Ca2+ with cytosolic and lumenal buffers [Bk] and [Gj], respectively, where k+ denotes the capture rate and k- is the dissociation rate.The cellular Ca2+ dynamics based on the CICR mechanism results from the local properties of the IP3R channels and their spatial coupling by Ca2+ diffusion that is influenced by buffer reactions and SERCA pump activity. Thus, the cellular Ca2+ signaling is determined by the channel states and corresponding release fluxes that shape the signals in dependence on the properties of the cellular diffusion. In turn, the resulting Ca2+ concentration profiles influence the channel opening behavior due to its dependence on the local Ca2+ concentration.Standard simulation techniques for reaction diffusion systems like finite difference or finite element methods have to calculate the concentration profiles for the complete numerical mesh for convergence criteria. The large concentration gradients typically observed close to open channels require a fine spatial discretization that leads in particular in three spatial dimensions to large computational costs [51]. To circumvent this restriction and enable simulations of biological relevant time courses, we developed an analytic solution of the reaction diffusion equations in terms of coupled Green’s function that allows for calculating [Ca2+]i at specific spatial positions only such as channel cluster locations. These local concentrations subsequently determine potential changes in the stochastic channel dynamics leading to modified source terms in the reaction diffusion system.In order to solve the reaction diffusion system analytically we restrict the system to one mobile and one immobile buffer, introduce dimensionless variables, and linearize the corresponding equations [27]. To obtain an analytic solution of the resulting reaction diffusion system we have to determine the geometry and boundary conditions. As a natural choice we use a spherical geometry and, motivated by the observation that astrocytic Ca2+ signals do not change significantly in Ca2+ free medium [24], no-flux boundary condition that reflects the properties of the plasma membrane for our in silico cell. With these specifications we could develop the analytic solution by coupled Green’s functions [52, 53] describing Ca2+ and buffer dynamics where open channels correspond to source terms. For mathematical details of the implemented Green’s cell algorithm (GCA) see [27] where we have shown how microscopic channel fluctuations are transmitted onto the level of the cell by diffusion thereby matching the experimental findings on cell variability and buffer influence.For the IP3R description within the GCA, we use the DeYoung-Keizer model (DKM) that assumes three binding sites for each of the four subunits per channel [54]. Each of the three binding sites can be free or occupied leading to 23 different states Xijk and 12 possible transitions in dependence on the Ca2+ concentration C and IP3 concentration I that can be visualized on a cube as shown in Figure 2(a). A subunit is active in the state X110 only where the first index equals one if IP3 is bound, the second index corresponds to Ca2+ binding to an activating site, and the last index describes the occupation of the dominant inhibiting Ca2+ binding site. The activating binding sites for Ca2+ have a higher affinity compared to the dominant inhibiting site. This is a minimal choice to generate the nonlinearity in Ca2+, which is the source of the bell shaped open probability (Figure 2(b)) and the basis of CICR [25].Figure 2
CICR is based on the nonlinear open probability of IP3R. (a) Scheme of the DeYoung-Keizer model for a single subunit. A subunit is active, if IP3 (I) is bound and Ca2+ (C) is only bound to the activating site, that is, in state X110. A channel opens if at least 3 of its 4 subunits are active. Open IP3Rs release Ca2+ that can diffuse to other IP3R clusters and activate the channels there as well, eventually leading to cell wide Ca2+ release and generating a [Ca2+]i spike. From the cytosol Ca2+ is pumped back into the ER by SERCA that decreases diffusion. (b) From the Markov chain model of a single subunit the non-Markovian stationary open probability Po of a channel can be calculated in dependence on Ca2+ and IP3 concentration [25].
(a)
(b)


### 2.2.2. Modeling the Impact of SERCA Pump Strength on Ca2+ Dynamics

Here we use our validated GCA implementation to test the hypothesis on the temperature dependent SERCA activity as a main reason for the observed differences in[Ca2+]i spiking. For the simulations, we use physiological cellular diffusion parameters listed in Table 1. Our in silico cell has 31 randomly distributed IP3R clusters each of them having a random number of channels between 2 and 16, yielding 258 channels in total. The clusters are separated by at least 1.5 μm and only located within a sphere with 80% of the cell radius R to guarantee convergence of Green’s function representation. To test the SERCA hypothesis, this prototype of a cell was simulated with different pump strengths Pp.Table 1
Physiologic parameters of the in silico cell used in simulations. Note that the GCA uses a linearized reaction diffusion system leading to rates forσl and Pp.
Physiological cell parameter values for GCA simulations
R
10μm
Cell radius
γ
5
Volume ratio between cytosol and ER
D
Ca
220μm2/s
Diffusion coefficient of cytosolicCa2+
D
E
70μm2/s
Diffusion coefficient of lumenalCa2+
D
B
95μm2/s
Diffusion coefficient of mobile buffer
[
C
a
2
+
]
0
55 nM
CytosolicCa2+ base level
[
E
]
0
800μM
Endoplasmic reticulumCa2+ base level
[
I
P
3
]
80μM
IP
3 concentration
[
B
]
T
52μM
Total mobile buffer concentration (EGTA)
k
B
+
0.3 (μMs)−1
On rate of the mobile buffer
k
B
-
1.5 s−1
Dissociation rate of the mobile buffer
[
B
i
]
T
30μM
Total immobile buffer concentration
k
B
i
+
600 (μMs)−1
On rate of the immobile buffer
k
B
i
-
100 s−1
Dissociation rate of the immobile buffer
σ
l
≈0.01 s−1
Leak flux constant given by variedPp, [Ca2+]0 and [E]0Figure3 shows the cell wide dynamics, where panel (a) exhibits the number of open channels and panel (b) shows the resulting cytosolic Ca2+ concentration for pump strengths Pp varying between 22 s−1 and 200 s−1 what is around the estimated physiological reference value of 86 s−1 [55]. A first conspicuous property shown in panel (a) is that the amplitudes of open channels during a spike are rather similar and vary between 25 and 30, rather independently of the pump strengths Pp. Nevertheless, Pp has a huge influence on the time course.Figure 3
Dependence of Ca2+ signals on the pump strength. The temperature dependent activity of SERCA leads to higher pump rates Pp for higher temperatures. This influences the Ca2+ signals as shown by the number of open channels (a) and the corresponding Ca2+ concentrations (b). For the two smallest strengths the simulations ended before the shown real time of 800 s caused by run time restrictions on the compute cluster.For low pump strengths, Ca2+ is removed too slowly from the cytosol for a coordinated signal. During the initial activation from base level the amplitude of open channels becomes maximal and is limited by inhibition according to the channel properties explained in Figure 2. When a channel recovers from inhibition due to the stochastic unbinding of Ca2+ from a dominantly inhibiting site of a subunit, the cytosolic Ca2+ concentration is still in the Ca2+ concentration range where the higher affinity for the Ca2+ activating binding site leads to a large open probability shown in Figure 2(b). Subsequently, the channel will open again, increase the local Ca2+ concentration, and reinduce channel inhibition. Repeating this scenario, the channels will mainly switch between activation and inhibition back and forth where a rather constant proportion of up to 50% of subunits are inhibited until randomly no reactivation occurs and the cell is reaching its resting state again. This permanent switching between activation and inhibition leads to the observed plateaus for the two lowest pump strengths.This behavior is also reflected by the resulting cytosolic Ca2+ concentrations. For Pp=22 s−1, the cell exhibits a plateau response with some superimposed [Ca2+]i spiking, which are often observed in experiments with astrocytes from acute brain slices [24, 56–58]. For Pp=32 s−1, the shorter plateaus of open channels are translated into a more oscillation-like behavior of the cytosolic Ca2+ concentration exhibiting rather diverse spike amplitudes. The concentration is in a similar range as the cell with the lower pump strength, indicating a regime where local channel inhibition and spatial coupling determined by pumps have an equal influence on the cellular dynamics.For increasingPp, the cells exhibit an even more spiking like behavior since long lasting channel activity is absent, and therefore [Ca2+]i shows more pronounced peaks. This can be again explained by the interplay of inhibition and Ca2+ removal. After a spike has occurred and most channels are inhibited, the SERCAs decrease the cytosolic Ca2+ concentration sufficiently fast leading to base level [Ca2+]i when channels recover from inhibition. Subsequently, the cell relaxes to its resting level, from which it can be activated again by CICR, that is, by triggering opening of a channel cluster. For very large pump strengths, the spikes of the number of open channels have still a similar height to those for smaller Pp, but their widths decrease drastically. This induces low and slim peaks of the cytosolic Ca2+ concentration.So far, the simulations seem to support the hypothesis of the decreased SERCA activity leading to an increased[Ca2+]i activity at lower temperatures. In terms of the experimental findings, a temperature increase in experiments and observed changes in the Ca2+ signals correspond to an increase of the pump rate from the region from 42 s−1 to 62 s−1 to 152  s−1 or higher.For a further detailed analysis, we determine the spike widths of the simulated cytosolic Ca2+ peaks as done for the experimental data in Figure 1. Figure 4(a) exhibits the dependence of the spike width on the pump strength Pp. The analysis reveals an exponential dependence on the pump rate Pp as shown by the line defined by SW=61s·e-0.012s·Pp. At second glance we observe that the spike width becomes most regular for Pp=62 s−1 as depicted by the small error bar that exhibits the smallest relative variation even compared to the smaller error bars for large Pp. This might indicate an optimal dynamical regime.Figure 4
Spike width and amplitude in dependence on the pump strength for the Ca2+ signal shown in Figure 3. (a) The spike width exhibits an exponential dependence on Pp. This is in good agreement with the experimental findings in Figure 1. (b) Also the amplitude dependence on Pp obeys an exponential relation (error bars denote s.e.).
(a)
(b)An analogue behavior can be found for the average amplitudes in Figure4(b). The dependence on Pp can be fitted by the exponential relation 2.9·e-0.011s·Pp, and again the mean amplitude for Pp=62 s−1 exhibits the smallest relative variations. The amplitude behavior obtained from simulations has a similar trend like those seen in the experiments but contradicts former theoretical results where the amplitude stayed constant [59].To analyze the potential optimal regime, we determine the oscillation characteristicsTav and σ. Indeed, Tav exhibits a minimum of 42 s for a pump rate of 62 s−1 as shown in Figure 5(a). Panel (b) displays the dependence of σ on Tav showing a slope close to one as in earlier experiments and simulations and indicates a deterministic time of 20 s. The spread of the data points is exclusively caused by the different pump strengths. Similar simulations with smaller buffer concentrations have shown an analogue behavior with a shift of the minimum of Tav to larger Pp. This is in line with our consideration, since a lower buffer concentration leads to an increased spatial coupling that can be compensated by higher pump rates.Figure 5
Dependence of the mean periodTav on the pump strength Pp. (a) Tav exhibits a nonlinear dependence on the pump strength with a minimal mean period of 42 s for Pp=62 s−1. This might illuminate a possible control mechanism cells can use to tune their oscillations. (b) The σ-Tav relation exhibits a population slope close to one and a deterministic time of approximately 20 s. (c) The height of the local Ca2+ concentration Ca2+ close to a cluster is hardly affected by Pp as demonstrated by the representative example of the same cluster in different simulations using different pump rates. Hence, the differences in the cytosolic Ca2+ oscillations in Figure 3(b) are spatially induced.
(a)
(b)
(c)For an additional test of our conception based on the spatial coupling of release sites, we analyze the local Ca2+ concentrations at channel clusters in Figure 5(c). For all pump strengths the concentration peaks are in the range of 100 μM independently of Pp. Hence, the different cell wide cytosolic Ca2+ signals shown in Figure 3(b) are not caused by different local properties but are caused by the spatial coupling of release sites. Thus, the results of the temperature experiments are consistent with the hypothesized spatial control mechanism.


### 2.2.3. Rate Equation Based Modeling

To further investigate the role of the spatial aspect in Ca2+ signaling, we model the effect of SERCA activity on the Ca2+ profile by a traditional rate equation model implemented by ordinary differential equations (ODEs). These types of models do not consider any spatial aspects but treat the cell as a well-stirred homogeneous entity and were extensively used in the establishment of the Ca2+ signaling mechanism [60, 61]. A variety of models have been proposed and applied to diverse aspects of Ca2+ signaling where the common mechanism is Ca2+ cycling between the ER and the cytosol. Here we use a model based on previous studies which reflects the main Ca2+ fluxes according to Figure 2 as well as Ca2+ buffering by mitochondria. The governing rate equations for the cytosolic and ER Ca2+ concentration read(2a)dCacdt=fcytJext+ΔJmito+Jrelease+Jleak-Jserca-Jpmca,(2b)dCaerdt=ferγ-1Jserca-Jrelease-Jleak,where fcyt and fer are the fractions of free Ca2+ within the cytosol and the ER [62], γ is again the volume ratio between cytosol and ER, Δ denotes the mitochondria to cytosol volume ratio, and the corresponding fluxes J reflect typical physiological characteristics. The 2 most important fluxes are the Ca2+ release from the ER by IP3R (Jrelease) and the SERCA mediated Ca2+ removal from the cytosol into the ER (Jserca) given by (3a)Jrelease=PrIP32p12+IP32Cac2p22+Cac2·p34p44+Cac4Caer-Cac,(3b)Jserca=PpCac2p52+Cac2,respectively. The dominant flux for the [Ca2+]i spiking is the nonlinear release flux inducing the CICR by the bell shaped dependence on the Ca2+ concentration [63] here modeled in analogy to previous approaches [64]. The form of the complementary SERCA mediated Ca2+ flux into the ER (Jserca) reflects the fact that two Ca2+ ions are pumped into the ER by the use of one ATP molecule where the ATP concentration is assumed to be constant. The remaining fluxes are given in Table 2 where Jext and Jmito describe Ca2+ flux from extracellular space and mitochondria, respectively. Jleak corresponds to the Ca2+ leakage from ER and Jpmca is the flux into the extracellular space through plasma membrane Ca2+ ATPases. The corresponding parameters pi are also listed in Table 2.Table 2
Physiologic parameters of the in silico cell used for the rate equation simulations.
Additional flux definitions of the full rate equation model
Sodium-Ca2+ exchanger mediated flux
J
NaCa
=
p
5
C
a
m
C
a
m
+
p
6
p
7
2
p
8
2
+
C
a
c
2
Mitochondrial unitransporter mediated flux
J
uni
=
p
9
C
a
c
2
p
10
2
+
C
a
c
2
Effective mitochondrial flux
J
mito
=
J
NaCa
-
J
uni
Leak flux through ER membrane
J
leak
=
p
11
C
a
er
-
C
a
c
Influx from external space
J
ext
=
J
ext,bas
-
(
p
12
Ca
er
+
p
13
Ca
c
)
PMCA mediated flux
J
pmca
=
p
14
C
a
c
2
p
15
2
+
C
a
c
2
Parameter values for the rate equation calcium model
P
r = 0.03 ms−1
p
1 = 0.5 μM
p
2 = 0.9 μM
p
3 = 1.3 μM
p
4 = 0.2 μM
P
p = 0.06 μM ms−1
p
5 = 10 μM ms−1
p
6 = 0.6 μM
p
7= 1 μM
p
8 = 10 μM
p
9 = 0.2 μM ms−1
p
10 = 1 μM
p
11 = 10-3 ms−1
p
12 = 3·10-6 ms−1
p
13 = 15·10-6 ms−1
p
14 = 0.01 μM ms−1
p
15 = 0.7 μM
IP
3 = 0.4 μM
f
cyt
=
0.015
f
er = 0.5
J
ext,bas = 3.5·10-3 μM ms−1
Δ
=
0.07
γ
-
1
=
30To investigate the influence of SERCA pump strength on Ca2+ dynamics, we simulate the ODEs for a wide range of pump strength Pp and three different values of the Ca2+ release flux (Pr) from IP3R channels and analyze the spike characteristics (Figure 6). This parameter scan shows that the Ca2+ period within the rate equation model exhibits a similar behavior as the spatial GCA simulations in terms of a minimal period in dependence on the pump strength (Figure 6(a)). While this finding is not in clear favor for the perspective of spatial coupling of release sites as the [Ca2+]i spike generating mechanism, the results obtained with the rate equation model given by (2a) and (2b) about the spike width and amplitude exhibit a different behavior than the experimental data.Figure 6
Cell characteristics in the rate equation model. Dependence of the period (a), amplitude (b), and spike width (c) on the SERCA pump strengthPp for three different Ca2+ release rates Pr from IP3 channels.
(a)
(b)
(c)In the ODE model a higher SERCA pump rate causes more Ca2+ within the ER and also a higher Ca2+ release into the cytosol because IP3R flux is controlled to a large proportion by the concentration difference and thus by the ER concentration. A high SERCA pump rate leads to a larger spike width as observed in Figure 6(b). These higher fluxes in and out of the ER also yield higher spike amplitudes (Figure 6(c)). Both observations are in contrast to the experimental observations and the GCA simulations, potentially indicating the fact that a spatial distribution of channels is a necessity for the experimentally observed SERCA pump effects.


## 2.2.1. Multiscale Modeling Approach

Based on the observations that Ca2+ spiking is stochastic and obeys a nontrivial signaling signature, our multiscale model considers the molecular dynamics of individual IP3Rs and Ca2+ diffusion as well as buffer dynamics [27]. Our modeling strategy for this coupled system is based on separation of the two involved length scales. On the microscopic scale we use a detailed model for the channels [25]. Those will be open and close in a stochastic manner due to stochastic binding of the signaling molecules IP3 and Ca2+. The stochastic channel dynamics acts as a noisy source term that is driving the three-dimensional reaction diffusion system (RDS) representing the cytosol. Thereby the model also takes reactions of Ca2+ with buffer molecules like BAPTA or calmodulin and the SERCA pumps into account. From this perspective the cellular dynamics is given by an in general nonlinear reaction diffusion system that considers the cytosolic and lumenal dynamics by (1a)∂[Ca2+]∂t=Dc∇2[Ca2+]+σ~l+Pc(r,t)·([E]-[Ca2+])-P~pCa2+2kd2+Ca2+2+∑j‍kj-CaBj-kj+BjCa2+,(1b)∂[E]∂t=DE∇2[E]-γCa2+2kd2+Ca2+2σ~l+Pcr,tE-Ca2+-P~pCa2+2kd2+Ca2+2+∑k‍kk-EGk-kk+GkE,(1c)∂[Bj]∂t=DBj∇2Bj-kj+BjCa2++kj-CaBj,(1d)∂[Gk]∂t=DGk∇2[Gk]-kk+[Gk][E]+kk-[EGk],where the first two equations describe the free Ca2+ concentration in the cytosol and the ER, respectively, and the other two equations correspond to a variety of cytosolic and lumenal buffers [Bj] and [Gk]. The first terms on the right hand side characterize diffusion with the individual diffusion coefficients Di. The following terms in (1a) and (1b) describe the leak flux (with permeability σl) and release through channels (Pc). The space dependence of Pc(rc,t) reflects the location of channel clusters and the time dependent stochastic opening and closing of the channels. This time dependent source term is governed by stochastic simulation of the channel states. γ is the volume ratio Vcyt/VER. The third term results from Ca2+ uptake into the ER from the cytosol by SERCA pumps with a flux strength Pp and the SERCA dissociation constant kd. The remaining terms in (1a)–(1d) describe the reactions of Ca2+ with cytosolic and lumenal buffers [Bk] and [Gj], respectively, where k+ denotes the capture rate and k- is the dissociation rate.The cellular Ca2+ dynamics based on the CICR mechanism results from the local properties of the IP3R channels and their spatial coupling by Ca2+ diffusion that is influenced by buffer reactions and SERCA pump activity. Thus, the cellular Ca2+ signaling is determined by the channel states and corresponding release fluxes that shape the signals in dependence on the properties of the cellular diffusion. In turn, the resulting Ca2+ concentration profiles influence the channel opening behavior due to its dependence on the local Ca2+ concentration.Standard simulation techniques for reaction diffusion systems like finite difference or finite element methods have to calculate the concentration profiles for the complete numerical mesh for convergence criteria. The large concentration gradients typically observed close to open channels require a fine spatial discretization that leads in particular in three spatial dimensions to large computational costs [51]. To circumvent this restriction and enable simulations of biological relevant time courses, we developed an analytic solution of the reaction diffusion equations in terms of coupled Green’s function that allows for calculating [Ca2+]i at specific spatial positions only such as channel cluster locations. These local concentrations subsequently determine potential changes in the stochastic channel dynamics leading to modified source terms in the reaction diffusion system.In order to solve the reaction diffusion system analytically we restrict the system to one mobile and one immobile buffer, introduce dimensionless variables, and linearize the corresponding equations [27]. To obtain an analytic solution of the resulting reaction diffusion system we have to determine the geometry and boundary conditions. As a natural choice we use a spherical geometry and, motivated by the observation that astrocytic Ca2+ signals do not change significantly in Ca2+ free medium [24], no-flux boundary condition that reflects the properties of the plasma membrane for our in silico cell. With these specifications we could develop the analytic solution by coupled Green’s functions [52, 53] describing Ca2+ and buffer dynamics where open channels correspond to source terms. For mathematical details of the implemented Green’s cell algorithm (GCA) see [27] where we have shown how microscopic channel fluctuations are transmitted onto the level of the cell by diffusion thereby matching the experimental findings on cell variability and buffer influence.For the IP3R description within the GCA, we use the DeYoung-Keizer model (DKM) that assumes three binding sites for each of the four subunits per channel [54]. Each of the three binding sites can be free or occupied leading to 23 different states Xijk and 12 possible transitions in dependence on the Ca2+ concentration C and IP3 concentration I that can be visualized on a cube as shown in Figure 2(a). A subunit is active in the state X110 only where the first index equals one if IP3 is bound, the second index corresponds to Ca2+ binding to an activating site, and the last index describes the occupation of the dominant inhibiting Ca2+ binding site. The activating binding sites for Ca2+ have a higher affinity compared to the dominant inhibiting site. This is a minimal choice to generate the nonlinearity in Ca2+, which is the source of the bell shaped open probability (Figure 2(b)) and the basis of CICR [25].Figure 2
CICR is based on the nonlinear open probability of IP3R. (a) Scheme of the DeYoung-Keizer model for a single subunit. A subunit is active, if IP3 (I) is bound and Ca2+ (C) is only bound to the activating site, that is, in state X110. A channel opens if at least 3 of its 4 subunits are active. Open IP3Rs release Ca2+ that can diffuse to other IP3R clusters and activate the channels there as well, eventually leading to cell wide Ca2+ release and generating a [Ca2+]i spike. From the cytosol Ca2+ is pumped back into the ER by SERCA that decreases diffusion. (b) From the Markov chain model of a single subunit the non-Markovian stationary open probability Po of a channel can be calculated in dependence on Ca2+ and IP3 concentration [25].
(a)
(b)


## 2.2.2. Modeling the Impact of SERCA Pump Strength on Ca2+ Dynamics

Here we use our validated GCA implementation to test the hypothesis on the temperature dependent SERCA activity as a main reason for the observed differences in[Ca2+]i spiking. For the simulations, we use physiological cellular diffusion parameters listed in Table 1. Our in silico cell has 31 randomly distributed IP3R clusters each of them having a random number of channels between 2 and 16, yielding 258 channels in total. The clusters are separated by at least 1.5 μm and only located within a sphere with 80% of the cell radius R to guarantee convergence of Green’s function representation. To test the SERCA hypothesis, this prototype of a cell was simulated with different pump strengths Pp.Table 1
Physiologic parameters of the in silico cell used in simulations. Note that the GCA uses a linearized reaction diffusion system leading to rates forσl and Pp.
Physiological cell parameter values for GCA simulations
R
10μm
Cell radius
γ
5
Volume ratio between cytosol and ER
D
Ca
220μm2/s
Diffusion coefficient of cytosolicCa2+
D
E
70μm2/s
Diffusion coefficient of lumenalCa2+
D
B
95μm2/s
Diffusion coefficient of mobile buffer
[
C
a
2
+
]
0
55 nM
CytosolicCa2+ base level
[
E
]
0
800μM
Endoplasmic reticulumCa2+ base level
[
I
P
3
]
80μM
IP
3 concentration
[
B
]
T
52μM
Total mobile buffer concentration (EGTA)
k
B
+
0.3 (μMs)−1
On rate of the mobile buffer
k
B
-
1.5 s−1
Dissociation rate of the mobile buffer
[
B
i
]
T
30μM
Total immobile buffer concentration
k
B
i
+
600 (μMs)−1
On rate of the immobile buffer
k
B
i
-
100 s−1
Dissociation rate of the immobile buffer
σ
l
≈0.01 s−1
Leak flux constant given by variedPp, [Ca2+]0 and [E]0Figure3 shows the cell wide dynamics, where panel (a) exhibits the number of open channels and panel (b) shows the resulting cytosolic Ca2+ concentration for pump strengths Pp varying between 22 s−1 and 200 s−1 what is around the estimated physiological reference value of 86 s−1 [55]. A first conspicuous property shown in panel (a) is that the amplitudes of open channels during a spike are rather similar and vary between 25 and 30, rather independently of the pump strengths Pp. Nevertheless, Pp has a huge influence on the time course.Figure 3
Dependence of Ca2+ signals on the pump strength. The temperature dependent activity of SERCA leads to higher pump rates Pp for higher temperatures. This influences the Ca2+ signals as shown by the number of open channels (a) and the corresponding Ca2+ concentrations (b). For the two smallest strengths the simulations ended before the shown real time of 800 s caused by run time restrictions on the compute cluster.For low pump strengths, Ca2+ is removed too slowly from the cytosol for a coordinated signal. During the initial activation from base level the amplitude of open channels becomes maximal and is limited by inhibition according to the channel properties explained in Figure 2. When a channel recovers from inhibition due to the stochastic unbinding of Ca2+ from a dominantly inhibiting site of a subunit, the cytosolic Ca2+ concentration is still in the Ca2+ concentration range where the higher affinity for the Ca2+ activating binding site leads to a large open probability shown in Figure 2(b). Subsequently, the channel will open again, increase the local Ca2+ concentration, and reinduce channel inhibition. Repeating this scenario, the channels will mainly switch between activation and inhibition back and forth where a rather constant proportion of up to 50% of subunits are inhibited until randomly no reactivation occurs and the cell is reaching its resting state again. This permanent switching between activation and inhibition leads to the observed plateaus for the two lowest pump strengths.This behavior is also reflected by the resulting cytosolic Ca2+ concentrations. For Pp=22 s−1, the cell exhibits a plateau response with some superimposed [Ca2+]i spiking, which are often observed in experiments with astrocytes from acute brain slices [24, 56–58]. For Pp=32 s−1, the shorter plateaus of open channels are translated into a more oscillation-like behavior of the cytosolic Ca2+ concentration exhibiting rather diverse spike amplitudes. The concentration is in a similar range as the cell with the lower pump strength, indicating a regime where local channel inhibition and spatial coupling determined by pumps have an equal influence on the cellular dynamics.For increasingPp, the cells exhibit an even more spiking like behavior since long lasting channel activity is absent, and therefore [Ca2+]i shows more pronounced peaks. This can be again explained by the interplay of inhibition and Ca2+ removal. After a spike has occurred and most channels are inhibited, the SERCAs decrease the cytosolic Ca2+ concentration sufficiently fast leading to base level [Ca2+]i when channels recover from inhibition. Subsequently, the cell relaxes to its resting level, from which it can be activated again by CICR, that is, by triggering opening of a channel cluster. For very large pump strengths, the spikes of the number of open channels have still a similar height to those for smaller Pp, but their widths decrease drastically. This induces low and slim peaks of the cytosolic Ca2+ concentration.So far, the simulations seem to support the hypothesis of the decreased SERCA activity leading to an increased[Ca2+]i activity at lower temperatures. In terms of the experimental findings, a temperature increase in experiments and observed changes in the Ca2+ signals correspond to an increase of the pump rate from the region from 42 s−1 to 62 s−1 to 152  s−1 or higher.For a further detailed analysis, we determine the spike widths of the simulated cytosolic Ca2+ peaks as done for the experimental data in Figure 1. Figure 4(a) exhibits the dependence of the spike width on the pump strength Pp. The analysis reveals an exponential dependence on the pump rate Pp as shown by the line defined by SW=61s·e-0.012s·Pp. At second glance we observe that the spike width becomes most regular for Pp=62 s−1 as depicted by the small error bar that exhibits the smallest relative variation even compared to the smaller error bars for large Pp. This might indicate an optimal dynamical regime.Figure 4
Spike width and amplitude in dependence on the pump strength for the Ca2+ signal shown in Figure 3. (a) The spike width exhibits an exponential dependence on Pp. This is in good agreement with the experimental findings in Figure 1. (b) Also the amplitude dependence on Pp obeys an exponential relation (error bars denote s.e.).
(a)
(b)An analogue behavior can be found for the average amplitudes in Figure4(b). The dependence on Pp can be fitted by the exponential relation 2.9·e-0.011s·Pp, and again the mean amplitude for Pp=62 s−1 exhibits the smallest relative variations. The amplitude behavior obtained from simulations has a similar trend like those seen in the experiments but contradicts former theoretical results where the amplitude stayed constant [59].To analyze the potential optimal regime, we determine the oscillation characteristicsTav and σ. Indeed, Tav exhibits a minimum of 42 s for a pump rate of 62 s−1 as shown in Figure 5(a). Panel (b) displays the dependence of σ on Tav showing a slope close to one as in earlier experiments and simulations and indicates a deterministic time of 20 s. The spread of the data points is exclusively caused by the different pump strengths. Similar simulations with smaller buffer concentrations have shown an analogue behavior with a shift of the minimum of Tav to larger Pp. This is in line with our consideration, since a lower buffer concentration leads to an increased spatial coupling that can be compensated by higher pump rates.Figure 5
Dependence of the mean periodTav on the pump strength Pp. (a) Tav exhibits a nonlinear dependence on the pump strength with a minimal mean period of 42 s for Pp=62 s−1. This might illuminate a possible control mechanism cells can use to tune their oscillations. (b) The σ-Tav relation exhibits a population slope close to one and a deterministic time of approximately 20 s. (c) The height of the local Ca2+ concentration Ca2+ close to a cluster is hardly affected by Pp as demonstrated by the representative example of the same cluster in different simulations using different pump rates. Hence, the differences in the cytosolic Ca2+ oscillations in Figure 3(b) are spatially induced.
(a)
(b)
(c)For an additional test of our conception based on the spatial coupling of release sites, we analyze the local Ca2+ concentrations at channel clusters in Figure 5(c). For all pump strengths the concentration peaks are in the range of 100 μM independently of Pp. Hence, the different cell wide cytosolic Ca2+ signals shown in Figure 3(b) are not caused by different local properties but are caused by the spatial coupling of release sites. Thus, the results of the temperature experiments are consistent with the hypothesized spatial control mechanism.


## 2.2.3. Rate Equation Based Modeling

To further investigate the role of the spatial aspect in Ca2+ signaling, we model the effect of SERCA activity on the Ca2+ profile by a traditional rate equation model implemented by ordinary differential equations (ODEs). These types of models do not consider any spatial aspects but treat the cell as a well-stirred homogeneous entity and were extensively used in the establishment of the Ca2+ signaling mechanism [60, 61]. A variety of models have been proposed and applied to diverse aspects of Ca2+ signaling where the common mechanism is Ca2+ cycling between the ER and the cytosol. Here we use a model based on previous studies which reflects the main Ca2+ fluxes according to Figure 2 as well as Ca2+ buffering by mitochondria. The governing rate equations for the cytosolic and ER Ca2+ concentration read(2a)dCacdt=fcytJext+ΔJmito+Jrelease+Jleak-Jserca-Jpmca,(2b)dCaerdt=ferγ-1Jserca-Jrelease-Jleak,where fcyt and fer are the fractions of free Ca2+ within the cytosol and the ER [62], γ is again the volume ratio between cytosol and ER, Δ denotes the mitochondria to cytosol volume ratio, and the corresponding fluxes J reflect typical physiological characteristics. The 2 most important fluxes are the Ca2+ release from the ER by IP3R (Jrelease) and the SERCA mediated Ca2+ removal from the cytosol into the ER (Jserca) given by (3a)Jrelease=PrIP32p12+IP32Cac2p22+Cac2·p34p44+Cac4Caer-Cac,(3b)Jserca=PpCac2p52+Cac2,respectively. The dominant flux for the [Ca2+]i spiking is the nonlinear release flux inducing the CICR by the bell shaped dependence on the Ca2+ concentration [63] here modeled in analogy to previous approaches [64]. The form of the complementary SERCA mediated Ca2+ flux into the ER (Jserca) reflects the fact that two Ca2+ ions are pumped into the ER by the use of one ATP molecule where the ATP concentration is assumed to be constant. The remaining fluxes are given in Table 2 where Jext and Jmito describe Ca2+ flux from extracellular space and mitochondria, respectively. Jleak corresponds to the Ca2+ leakage from ER and Jpmca is the flux into the extracellular space through plasma membrane Ca2+ ATPases. The corresponding parameters pi are also listed in Table 2.Table 2
Physiologic parameters of the in silico cell used for the rate equation simulations.
Additional flux definitions of the full rate equation model
Sodium-Ca2+ exchanger mediated flux
J
NaCa
=
p
5
C
a
m
C
a
m
+
p
6
p
7
2
p
8
2
+
C
a
c
2
Mitochondrial unitransporter mediated flux
J
uni
=
p
9
C
a
c
2
p
10
2
+
C
a
c
2
Effective mitochondrial flux
J
mito
=
J
NaCa
-
J
uni
Leak flux through ER membrane
J
leak
=
p
11
C
a
er
-
C
a
c
Influx from external space
J
ext
=
J
ext,bas
-
(
p
12
Ca
er
+
p
13
Ca
c
)
PMCA mediated flux
J
pmca
=
p
14
C
a
c
2
p
15
2
+
C
a
c
2
Parameter values for the rate equation calcium model
P
r = 0.03 ms−1
p
1 = 0.5 μM
p
2 = 0.9 μM
p
3 = 1.3 μM
p
4 = 0.2 μM
P
p = 0.06 μM ms−1
p
5 = 10 μM ms−1
p
6 = 0.6 μM
p
7= 1 μM
p
8 = 10 μM
p
9 = 0.2 μM ms−1
p
10 = 1 μM
p
11 = 10-3 ms−1
p
12 = 3·10-6 ms−1
p
13 = 15·10-6 ms−1
p
14 = 0.01 μM ms−1
p
15 = 0.7 μM
IP
3 = 0.4 μM
f
cyt
=
0.015
f
er = 0.5
J
ext,bas = 3.5·10-3 μM ms−1
Δ
=
0.07
γ
-
1
=
30To investigate the influence of SERCA pump strength on Ca2+ dynamics, we simulate the ODEs for a wide range of pump strength Pp and three different values of the Ca2+ release flux (Pr) from IP3R channels and analyze the spike characteristics (Figure 6). This parameter scan shows that the Ca2+ period within the rate equation model exhibits a similar behavior as the spatial GCA simulations in terms of a minimal period in dependence on the pump strength (Figure 6(a)). While this finding is not in clear favor for the perspective of spatial coupling of release sites as the [Ca2+]i spike generating mechanism, the results obtained with the rate equation model given by (2a) and (2b) about the spike width and amplitude exhibit a different behavior than the experimental data.Figure 6
Cell characteristics in the rate equation model. Dependence of the period (a), amplitude (b), and spike width (c) on the SERCA pump strengthPp for three different Ca2+ release rates Pr from IP3 channels.
(a)
(b)
(c)In the ODE model a higher SERCA pump rate causes more Ca2+ within the ER and also a higher Ca2+ release into the cytosol because IP3R flux is controlled to a large proportion by the concentration difference and thus by the ER concentration. A high SERCA pump rate leads to a larger spike width as observed in Figure 6(b). These higher fluxes in and out of the ER also yield higher spike amplitudes (Figure 6(c)). Both observations are in contrast to the experimental observations and the GCA simulations, potentially indicating the fact that a spatial distribution of channels is a necessity for the experimentally observed SERCA pump effects.


## 3. Discussion

Previous experimental findings have shown a clear temperature dependency of Ca2+ signals in astrocytes and indicated the different temperatures at which cultured astrocytes and astrocytes in acute brain slices are typically measured as a potential explanation for the reported differences in Ca2+ signaling [24]. Based on these observation and experimental studies reporting a temperature dependent SERCA activity, we used here multiscale modeling to investigate if this activity change may explain the experimental finding.For this purpose we simulated in silico cells with our GCA implementation that allows for controlling all physiological parameters and obtaining mechanistic insights. From the activity patterns of identical in silico cells that only differ in the SERCA activity shown in Figure3, we can directly conclude that the interplay of the IP3R properties with respect to CICR and inhibition together with the SERCA pumps can generate a wide spectrum of cellular Ca2+ signals. In particular, [Ca2+]i spiking occurs for an intermediate pump rate where Ca2+ is removed fast enough from the cytosol to prevent too frequent repetitive stimulation but where the spatial coupling between the release sites is preserved to allow for cell wide coordinated signals.As a more quantitative analysis we analyzed the spike width in the simulation data in analogy to the experimental findings where we found that the experimental spike width followed an exponential dependency on temperature independently if the width was determined from cultured astrocytes or from astrocytes in acute brain slices (Figure1). Remarkably, we obtained again an exponential dependency of the spike width from the simulations as shown in Figure 4. The modeling based exponential relation does not depend explicitly on the temperature but depends on the pump strength Pp that was varied in simulations. From the obtained two exponential dependencies of the spike width on the temperature in experiments g(T)=aexp⁡-bT and on the SERCA activity in simulations f(Pp)=αexp⁡-βPp we can infer the relation between these two quantities by setting g(T)∝f(Pp) as(4)PpT=bβT+ln⁡α/aβ.This relation predicts a linear relation between the SERCA activityPp and the temperature T as shown in Figure 7 where the model inference is shown by the line. Note that the proportional assumption can lead to an additional scaling factor in the second term of (4) but not in the first one that describes the slope of the relation. A linear relation was also measured experimentally in tuna by others [26]. Within this study no absolute values of the pump strength were measured but relative changes for different temperatures compared to a reference at 30°C were determined as shown in Figure 7 by the data points. The qualitative agreement between the directly measured relative changes (data points in Figure 7) and the spike width based model inference (line in Figure 7) is another strong evidence that modified SERCA activity may induce the change in astrocytic Ca2+ signals at different temperatures.Figure 7
SERCA pump activity depends linearly on the temperature as shown independently by our model inference (line) and by experimental data from tuna species (data points from [26], Figure 2(b)). For the experimental data points the pump strength is set to 100% at 30°C and used for normalisation for the lower temperatures leading to relative activity values shown on the left vertical axis. As a comparison the inference from our simulation (4) is shown by the straight line with the absolute values given on the right vertical axis.It is worth emphasizing that we varied in our multiscale simulations the pump strength motivated by the experimentally reported temperature dependent SERCA activity but did not include any specific relation. From the unbiased parameter scan onPp we first observed the general tendency of less coordinated Ca2+ signals for higher temperature in accordance with the experimental data. More strikingly, we found an exponential dependency on Pp of the spike width in the simulation data in analogy to the exponential dependency of the spike width on the temperature in experiments. The exponential relation was not included explicitly within the model but followed naturally from the physiological simulations indicating the importance of the spatial signal organization. With the two obtained exponential relations for the spike width g(T) from experiments and f(Pp) from simulations, we then inferred the linear relation between SERCA activity Pp and temperature T shown by the line in Figure 7. This inferred linear relation is in great qualitative accordance with the previous directly measured temperature dependency of the SERCA activity shown by the data dots in Figure 7. While the previous study has motivated the modeling strategy of modulating the pump strength Pp, no specific relation has been assumed but is only caused by the spatial signal mechanism.Strikingly, the two independently obtained relations for the SERCA pump activity both exhibit a linear dependency on the temperature questioning the applicability of the Arrhenius law in this case.The finding either indicates an activation energyEact in the order of the thermal energy RT or a more complex reaction system than what the Arrhenius law is considering. While the activation energy of SERCA is larger than the thermal energy as indicated by the use of ATP for the pump reaction, the underlying effective reaction scheme is more complex. The pump reaction consists of several (potentially interdependent) binding steps of Ca2+ and ATP to the SERCA protein. Since each of the binding and unbinding processes will be affected by temperature changes, the resulting reaction may exhibit already different characteristics from the collision of free particles assumed in Arrhenius law. More importantly, we determine here the cell wide averaged SERCA activity. The modulated SERCA activity is changing the spatial coupling between the release sites which leads to a modified hierarchical signal organisation by changing diffusion properties within the diffusion mediated CICR mechanism. This renders the reaction scheme even more complex and may be the main reason for the deviation from Arrhenius law.When comparing the different scales for the pump strength-temperature dependence in Figure7, it becomes obvious that the used SERCA activity parameters Pp in the simulations were smaller than those we would expect from the inference. Furthermore, the slope of the relation from inference is twice as large as the linear fit to the experimental data. This quantitative discrepancy is potentially caused by slightly different spatial characteristic like Ca2+ buffer content and physiologic settings of the cells. Nevertheless, the integrative approach shows how molecular properties can be characterized from global cell behavior by mechanistic modeling.To further evaluate the spatial characteristic of the phenomenon, we were also analyzing local concentrations at channel clusters. Their amplitudes were rather unaffected by the pump strength but the collective IP3R activity was significant modulated. This is a further evidence that it is the cell level that shapes the [Ca2+]i signals and not local properties. In line with these findings, we found that the rate equation based model (2a) and (2b) is unable to recapitulate the experimental finding when modulating the pump activity. This model approach neglects the spatial scale and thus the spatial coupling between release sites. The failure to model the experimental observations can be also seen as an indirect support for the spatial characteristic of Ca2+ signals. In general, the observed experimental behaviors might be also obtained by specific realizations of rate equation models that might consider slightly different flux definitions or integrate further pathways like ATP consumption by SERCA pumps and facilitated ATP production by Ca2+ triggered mitochondrial activity. These models are typically rather focused on specific scenarios and often not able to generate the wide spectrum of observed Ca2+ signals. In contrast, we have shown that the multiscale approach of the GCA is able to generate all these experimental observations in dependence on reasonable physiological parameters [27].Interestingly, we found in the multiscale GCA simulations a regime of optimal signaling in terms of regularity of Ca2+ spikes. The nonlinear dependence of Tav on the pump strength might indicate a possible control mechanism. A recent study [65] has reported how frequency modulation may control gene regulation. If we therefore assume that the ability to spike and to use frequency coding is the purpose of the Ca2+ signaling pathway, cells can control this behavior by the expression level of SERCA. From our findings here, we would expect a negative feedback of SERCA expression on fast oscillations and a positive feedback on slow oscillations with high concentration peaks and plateaus.The impact of temperature on Ca2+ signaling is probably not an astrocyte-specific phenomenon. In cardiac muscle the frequency of Ca2+ spikes decreases dramatically at 37°C compared with 22°C. In other cell types like in rabbit renal tubules, hepatocytes, and parenchymal and endothelial cells acute hypothermia affects the intracellular Ca2+ homeostasis and goes along with a rise in cytosolic Ca2+ levels [66–68]. Due to the wide distribution of SERCA and the universality of Ca2+ signaling, this mechanism might be used in many temperature sensitive cellular processes, which could make SERCA a key temperature sensor.



---

*Source: 683490-2015-08-04.xml*