Nonlinear Gap Junctions Enable Long-Distance Propagation of Pulsating Calcium Waves in Astrocyte Networks

A new paradigm has recently emerged in brain science whereby communications between glial cells and neuron-glia interactions should be considered together with neurons and their networks to understand higher brain functions. In particular, astrocytes, the main type of glial cells in the cortex, have been shown to communicate with neurons and with each other. They are thought to form a gap-junction-coupled syncytium supporting cell-cell communication via propagating Ca2+ waves. An identified mode of propagation is based on cytoplasm-to-cytoplasm transport of inositol trisphosphate (IP3) through gap junctions that locally trigger Ca2+ pulses via IP3-dependent Ca2+-induced Ca2+ release. It is, however, currently unknown whether this intracellular route is able to support the propagation of long-distance regenerative Ca2+ waves or is restricted to short-distance signaling. Furthermore, the influence of the intracellular signaling dynamics on intercellular propagation remains to be understood. In this work, we propose a model of the gap-junctional route for intercellular Ca2+ wave propagation in astrocytes. Our model yields two major predictions. First, we show that long-distance regenerative signaling requires nonlinear coupling in the gap junctions. Second, we show that even with nonlinear gap junctions, long-distance regenerative signaling is favored when the internal Ca2+ dynamics implements frequency modulation-encoding oscillations with pulsating dynamics, while amplitude modulation-encoding dynamics tends to restrict the propagation range. As a result, spatially heterogeneous molecular properties and/or weak couplings are shown to give rise to rich spatiotemporal dynamics that support complex propagation behaviors. These results shed new light on the mechanisms implicated in the propagation of Ca2+ waves across astrocytes and the precise conditions under which glial cells may participate in information processing in the brain.


Introduction
The 20 th century witnessed crystallization of the neuronal doctrine, viewing neuron as the fundamental building block responsible for higher brain functions. Yet, neurons are not the only cells in the brain. In fact, almost 50% of the cells in the human brain are glial cells [1,2]. Due to their apparent lack of fast electrical excitability, the potential importance of glial cells in neural computation was downgraded in favor of the critical role played by these cells in neural metabolism. Recent experimental evidence however suggests that glial cells provide a role much more than support, including control of synapse function and formation, adult neurogenesis and regulation of cerebral blood flow (see e.g. [3] for a review). As a consequence, a new paradigm is emerging in brain science, according to which glial cells should be considered on a par with neurons.
In particular, astrocytes, the main type of glial cells in the cortex, have attracted much attention because they have been shown to communicate with neurons and with each other. Indeed, astrocytes can integrate neuronal inputs and modulate the synaptic activity between two neurons [4]. Neurotransmitters released from pre-synaptic neurons can bind to specific receptors on the astrocyte membrane and evoke Ca 2+ elevations in the astrocyte cytoplasm [5]. In turn, these activated astrocytes may release gliotransmitters, including glutamate and ATP, which feed back onto the synaptic terminals and modulate neuron responses [6].
Two main types of neuronal activity-dependent Ca 2+ responses are observed in astrocytes [7,8]: (1) transient Ca 2+ increases that are restricted to the very extremity of their distal processes [9,10] and (2) Ca 2+ elevations propagating along these processes as regenerative Ca 2+ waves, eventually reaching the cell soma. The latter kind of event can even propagate to neighboring astrocytes, thus forming intercellular Ca 2+ waves [11,12]. Although intercellular Ca 2+ waves have been extensively observed in astrocyte cultures [13,14], recent experimental evidence supports the possibility that they could also occur under physiological conditions [6], with propagation distances ranging from four [15] to up to 30 astrocytes [16]. These results therefore indicate that waves in astrocytes may represent an effective form of intercellular signaling in the central nervous system [17,18]. But further, they almost irresistibly bring about the hypothesis that this persistent astrocyte wave-based signaling could extend the repertoire of neural network communications, adding non-local interactions, both in space and in time [3].
In order to assess this hypothesis though, several aspects of Ca 2+ signaling in astrocytes remain to be elucidated. Experimental data suggest that a stimulus impinging on an astrocyte is preferentially encoded in the modulation of the frequency (FM) of astrocytic Ca 2+ oscillations [10]. This type of oscillations is often characterized by pulsating waves, i.e. the propagation of peak waveforms, with width smaller than period. However, the possibility of amplitude modulation (AM) or even coexisting AM and FM (AFM) encoding have also been inferred [19,20]. Actually, the frequency and amplitude of astrocytic Ca 2+ oscillations can be highly variable, depending on cell-specific properties such as Ca 2+ content of the intracellular stores, or the spatial distribution, density and activity of (sarco-)endoplasmic reticulum Ca 2+ -ATPase (SERCA) pumps [21,22]. Yet, the propagation of wave-like signalling in the context of such great variability is yet not fully understood [14].
Much effort has also been devoted to understand the mechanisms responsible for initiation and propagation of intercellular Ca 2+ waves. From a single-cell point of view, intracellular Ca 2+ dynamics in astrocytes is mainly due to Ca 2+ -induced Ca 2+ release (CICR) from the endoplasmic reticulum (ER) stores and its regulation by inositol trisphosphate (IP 3 ) [6]. But for the transmission of these internal signals from one astrocyte to the other, two possible routes have been uncovered. The first one involves the transfer of IP 3 molecules directly from the cytosol of an astrocyte to that of an adjacent one through gap junction intercellular hemichannels [23]. In the second route instead, propagation is mediated by extracellular diffusion of ATP which binds to plasma membrane receptors on neighboring astrocytes and regulates IP 3 levels therein [18,24]. Although these two routes need not be mutually exclusive, experiments indicated that intracellular propagation through gap junctions is likely the predominant signaling route in many astrocyte types [19,[25][26][27].
Albeit experimental protocols monitor wave propagation as variations of intracellular Ca 2+ , the molecule that is transmitted through gap junctions to neighboring astrocytes is not Ca 2+ , but IP 3 [27]. Indeed, when the IP 3 in a given cell increases, some of it can be transported through a gap junction to a neighbor astrocyte. This IP 3 surge in the neighbor cell can in turn trigger CICR, thus regenerating the original Ca 2+ signal. Yet, the transported IP 3 is required to reach a minimal threshold concentration to trigger CICR in the neighboring cell. If this threshold is not reached, propagation ceases [28]. In this regard, previous theoretical studies stressed the importance of a mechanism for at least partial regeneration of IP 3 levels [29,30]. Such a mechanism, coupled with IP 3 transport, could induce local IP 3 concentrations large enough to trigger CICR [30], thus enabling Ca 2+ wave propagation. Production of IP 3 by Ca 2+ -dependent PLCd has been suggested as a plausible candidate regeneration mechanism [29,31,32]. However, the intercellular latencies of the Ca 2+ waves simulated with this mechanism are hardly reconcilable with experimental observations, hinting a critical role for gap junction IP 3 permeability [29,30].
In the present study, we investigated the intercellular propagation of Ca 2+ waves through the gap-junctional route by a computer model of one-dimensional astrocyte network. To account for intracellular Ca 2+ dynamics, we adopted the concise realistic description of IP 3 -coupled Ca 2+ dynamics in astrocytes previously introduced in Ref. [33]. We specifically focused on the influence of gap junction linearity and internal Ca 2+ dynamics on the wave propagation distance. By means of bifurcation analysis and numerical solutions, we show that nonlinear coupling between astrocytes can indeed favor IP 3 partial regeneration thus promoting large-distance intercellular Ca 2+ wave propagation. Our study also shows that long-distance wave propagation critically depends on the nature of intracellular Ca 2+ encoding (i.e. whether Ca 2+ signals are FM or AM) and the spatial arrangement of the cells. Furthermore, our results suggest that, in the presence of weak coupling, nonlinear gap junctions could also explain the complex intracellular oscillation dynamics observed during intercellular Ca 2+ wave propagation in astrocyte networks [12].

The ChI model of intracellular Ca 2+ dynamics
We describe calcium dynamics in astrocytes by an extended version of the Li-Rinzel model [34], called the ChI model that we developed and studied in [33]. A detailed presentation of this model is also given in the Supplementary Information. Briefly, the ChI model accounts for the complex signaling pathway illustrated in Figure 1 that includes Ca 2+ regulation by IP 3 -dependent CICR as well as IP 3 dynamics resulting from PLCd-mediated synthesis and degradation by IP 3 3-kinase (3K) and inositol polyphosphate (IP) 5-phosphatase (5P). The temporal evolution of astrocytic intracellular calcium in our model is described by three coupled nonlinear equations: Author Summary In recent years, the focus of Cellular Neuroscience has progressively stopped only being on neurons but started to include glial cells as well. Indeed, astrocytes, the main type of glial cells in the cortex, dynamically modulate neuron excitability and control the flow of information across synapses. Moreover, astrocytes have been shown to communicate with each other over long distances using calcium waves. These waves spread from cell to cell via molecular gates called gap junctions, which connect neighboring astrocytes. In this work, we used a computer model to question what biophysical mechanisms could support long-distance propagation of Ca 2+ wave signaling. The model shows that the coupling function of the gap junction must be non-linear and include a threshold. This prediction is largely unexpected, as gap junctions are classically considered to implement linear functions. Recent experimental observations, however, suggest their operation could actually be more complex, in agreement with our prediction. The model also shows that the distance traveled by waves depends on characteristics of the internal astrocyte dynamics. In particular, long-distance propagation is facilitated when internal calcium oscillations are in their frequency-modulation encoding mode and are pulsating. Hence, this work provides testable experimental predictions to decipher long-distance communication between astrocytes.
in which the variables C, h, IP 3 represent the cell-averaged calcium concentration, the fraction of open IP 3 R channels on the ER membrane, and the cell-averaged concentration of IP 3 second messenger, respectively. Each one of these variables is coupled to others via the set of equations that describe contributions of different biochemical pathways, as described in details in Supplementary Information (equations S1-S4) alongside the complete mathematical analysis of the model features.
In a single-cell context, this model reproduces most of the available experimental data related to calcium oscillations in astrocytes. In particular, it faithfully reproduces the experimentally reported changes of oscillation frequency and wave shape caused by SERCA pump activity modulations [22].

Astrocyte coupling
Experimental evidence shows that chemical signaling between astrocytes usually takes the form of propagating Ca 2+ pulses that are elicited following the gap-junctional transfer of IP 3 second messenger molecules [13]. Intracellular IP 3 activates the CICR pathway, giving rise to the observed rapid transient elevations in cytosolic free calcium. We considered three scenarios to describe the exchange of IP 3 between a pair of adjacent astrocytes: (1) linear, (2) threshold-linear (composed of a linear term operating after a threshold) and (3) non-linear (here described as sigmoid) coupling (see Figure 2). The linear model is a simple diffusive coupling; however, threshold-linear and non-linear models both transfer IP 3 only when the IP 3 gradient between the two adjacent cells overcomes a threshold value.
Our investigation of nonlinear coupling case was motivated by the experimental observations suggesting that gap junction permeability in itself can be actively modulated by various factors, among them different second messengers. Indeed, there is growing evidence that gap junctions may have greater selectivity and more active gating properties than previously recognized [35]. Several signaling pathways are able to modulate junctional permeability. In particular, the conductance state of Cx43, the main type of connexin in astrocyte gap junctions [36] is regulated by phosphorylation by PKC, which is also involved in IP 3 degradation [37,38], as well as by intracellular Ca 2+ [39]. These data suggest that astrocyte gap junction gating could be coupled to intra-and inter-cellular IP 3 and Ca 2+ dynamics [40,41] in a nontrivial fashion. Accordingly, several previous simulation studies have explored the influence of complex (e.g. regulated by second messengers) gap junctions [42][43][44]. We explore here their effects on intracellular Ca 2+ wave propagation in astrocytes.
Linear coupling function. The linear model simply results from Fick's law of diffusion. The flux J i?j of IP 3 molecules (where i, j are indices of adjacent model astrocytes) is proportional to their concentration gradient: where D ij IP 3~I P 3 i {IP 3 j . Such coupling function is the standard model for a gap junction acting as a passive channel [30]. The coupling strength (or permeability) F depends on the number of gap junction channels and their unitary permeability, and in what follows, it will be considered as a parameter.
Non-linear coupling functions: sigmoid and thresholdlinear. Threshold-linear coupling only partially keeps the linear characteristics of the ''classical'' gap junction adding a threshold on IP 3 gradient below which the flux J i?j is zero. On the other hand, sigmoid coupling adds a further saturating threshold on the IP 3 gradient value, above which the IP 3 flux is constant.
Sigmoid coupling is defined as: where F is the coupling factor, IP 3 thr is the predetermined threshold value and IP 3 scale is the width of the transition zone in the sigmoid function (see Figure 2). In order to allow comparison of the effects of the threshold-linear coupling function with that of the sigmoid one, the slope of the threshold-linear function was chosen to be coincident with that of the sigmoid coupling function ( Figure 2). Threshold linear coupling is thus defined as: with the threshold function H x ð Þ~x if xw0, and 0 otherwise: The network model We consider chains of N astrocytes where each astrocyte is coupled to its two nearest neighbors via gap junctions. Each i-th astrocyte (i = 1,…,N) is associated with three variables C i , h i and IP 3 i , that are respectively the cytosolic Ca 2+ concentration, the ratio of open IP 3 Rs and the intracellular IP 3 concentration in this astrocyte. The dynamics of these internal variables is given by the ChI model (equations 1-3 and Supplementary Information for a detailed explanation): For all cells that are not at the boundaries of the astrocyte chain (i.e. Vi~2, . . . ,N{1): where the internal reaction term for IP 3 , i.e. Reflective (zero-flux) boundaries assume that IP 3 exiting cell 1 or N can only flow to cell 2 or N21, respectively. They are given by: We also considered the case where cells 1 and N are absorbing, namely they entrap incoming IP 3 fluxes. This is the case of absorbing boundary conditions in which IP 3 can flow from cell 2 to cell 1, but the reverse flux (from cell 1 to 2) is always null (and similarly for cells N and N21). Accordingly, equations read: Finally, with periodic boundary conditions, the 1D astrocyte chain actually takes the shape of a ring and the equations read:

Stimulation
To induce wave propagation in the astrocyte chain, one cell (referred to as the ''driving'' cell) is stimulated by a supplementary exogenous IP 3 input. This external stimulus is supplied through a (virtual) ''dummy'' cell, coupled to the driving cell by one of the coupling functions described above. In this sense, the dummy cell acts as an IP 3 reservoir in which the level of IP 3 is kept fixed to a constant value IP 3 bias . Let k be the coordinate of the stimulated cell (driving cell) within the 1D chain. In this study, we usually stimulate the first cell or the central one, that is k = 1 or k = (N+1)/2 (N odd). Hence, IP 3 dynamics in the k-th cell is given by Most simulations done in this work were driven by a constant value of IP 3 bias . In the last section though, a square positive wave stimulus was applied to the model.
Initial conditions for all cells were set in agreement with experimental values reported in astrocytes for Ca 2+ and IP 3 at basal conditions [45].

Numerics
The chain model consists of 3N non-linear ODEs, where the number of astrocytes in the chain, N, ranged from 1 to 100. Time solutions were obtained via numerical integration by a standard 4 th -order Runge-Kutta scheme with a time step of 10 ms as this value showed to be the best compromise between integration time and robustness of the results. The computational model was implemented in Matlab (2009a, The MathWorks, Natick, MA) and C. Bifurcation analysis was done using XPPAUT (http:// www.math.pitt.edu/,bard/xpp/xpp.html). Nonlinear time series analysis was performed using the TISEAN software package [46]. Table S1 in the Supplementary Information lists the values of the parameters used in the model.

Results
Before proceeding to study the propagation of calcium waves in spatially extended networks of astrocytes, it was necessary to understand the dynamical response of a single model cell in response to IP 3 stimulation. To this end, we performed a detailed bifurcation analysis of our model astrocytes. A wealth of dynamical regimes was discovered, allowing model astrocytes to encode information about IP 3 stimulus in amplitude-modulated (AM), frequency-modulated (FM) or mixed (AFM) modes, depending on parameter values (see Text S1 in the Supplementary Information and refs. [47,48]). We then proceeded to study the bifurcation diagrams for systems of coupled model astrocytes (utilizing different types of coupling as detailed in Methods). Briefly, the bifurcation analysis showed the existence of FM pulse-like oscillatory regimes at low IP 3 bias values, which can turn into complex oscillations for larger IP 3 bias values. Because stable oscillation regimes could coexist in the bifurcation diagrams with stable fixed points, it could not be predicted from these diagrams whether an IP 3 input to the cell would trigger pulse-like oscillations or not, i.e. whether it would switch the system from the fixed point to the oscillatory regime. Thus, we resorted to extensive numerical simulations to investigate under what conditions one could observe propagation of Ca 2+ waves along the astrocyte chain.
In agreement with previous studies (see [33] for a review), IP 3triggered CICR indeed allows intercellular Ca 2+ wave propagation in our modeling framework, as shown in Supplementary Information Figure S3. However, the range of wave propagation was usually restricted to and depended on the biophysical parameters that determine the profile of intracellular IP 3 dynamics [32]. In what follows, we delineate the role that these parameters play in Ca 2+ wave propagation, under different coupling modes (linear vs. nonlinear) and different encoding regimes (FM vs. AFM).

Calcium wave propagation along astrocyte chains
Linear vs. non-linear gap junctions. Propagation patterns both for the linear and nonlinear cases are presented in Figure 3 for the case of FM-encoding cells (N = 12). Here, a constant stimulus (IP 3 bias = 1.0 mM) was used and always applied to the first cell of the chain. Model analysis (see Methods) predicted that this level of IP 3 would trigger periodic Ca 2+ pulses at least in the stimulated cell and possibly in the other ones as indeed confirmed by simulations (see Astrocyte 1) both in the case of linear and nonlinear coupling. For the linear coupling case, we observed propagation failure at 6 th -7 th cell from the driving one (Figures 3a,b). By contrast, in the nonlinear coupling scenario, Ca 2+ pulses can propagate for the whole length of the chain (Figures 3c,d).
Analysis of the IP 3 pattern for the nonlinear coupling function ( Figure 3d) evidences a strong correlation between IP 3 and Ca 2+ pulses. The IP 3 pulses are followed in time by the Ca 2+ ones, suggesting that pulsed Ca 2+ propagation is mediated by the propagation of IP 3 across the cells. By contrast, in the case of linear coupling, the correlation between the propagating Ca 2+ pulses and the intracellular IP 3 signals is not so apparent (Figures 3a,b) as IP 3 seems to diffuse smoothly from the stimulated cell without any effective propagation pattern.
The observed difference in the propagation distance between linear and nonlinear gap-junction couplings can be understood from this analysis. Indeed, in the case of linear gap-junction coupled cells, the IP 3 arriving in cell i from cell i21 is transferred forward to cell i+1 before it can significantly accumulate in cell i. As a result, the IP 3 displays the almost diffusive pattern of Figure 3b, with a fast decay as the distance from the stimulated cell increases, and not real travelling wave structure in space. Hence even with large values of the coupling strength or stimulus intensity, beyond a limited number of cells away from the stimulated one, the IP 3 concentration becomes too small to trigger CICR. This stops Ca 2+ wave propagation. Conversely, with nonlinear gap junctions, IP 3 can accumulate in cell i (and trigger CICR) before it reaches the gap junction threshold and gets transferred to cell i+1. As a result, the IP 3 concentration evolves to the locally regenerative spatiotemporal pattern illustrated in Figure 3d) that allows Ca 2+ wave propagation over the whole network.
The distances (measured in units of number of cells) travelled by the propagating Ca 2+ waves as a function of the stimulation amplitude for an astrocyte chain of N = 25 FM-encoding cells are reported in Figure 4a. For linear gap junctions, the propagation distance increases with IP 3 bias , but never exceeds one third of the chain length. On the contrary, with nonlinear sigmoid coupling, Ca 2+ oscillations propagate along the whole chain as soon as the oscillatory regime is engaged (that is for IP 3 bias .0.72 mM, see Figure S2). Figure 4a further shows that threshold-linear gap junctions exhibit almost the same response as sigmoid ones but with different effective threshold IP 3 concentrations. Hence, these results indicate that the significant parameter for long-distance wave propagation through nonlinear gap junctions is the presence of an IP 3 concentration threshold below which the junction is closed (this property is shared by the two nonlinear models), rather than the saturation of the transport at high IP 3 concentrations (found only in the sigmoid model).
Because the effects due to the different shapes of coupling curves could be conflated in the above observations, we computed the dependence of wave propagation range on the maximal strength of coupling, for linear vs. nonlinear coupling cases (Figure 4c). For the most part of the range of examined coupling strengths, Ca 2+ wave propagation distance was significantly larger for the nonlinear case as compared to the linear case, ruling out the possibility that our findings are just a trivial confound. The only exceptions to this claim were noted for low F values (these dynamics at low coupling are studied thereafter) and for very small regions (around F = 1.5 and F = 2.5), were the propagation distances for both couplings were comparable. More importantly, Figure 4c evidences that linear gap junction fails to propagate long-distance Ca 2+ waves (except in the ''chaotic'' low F domain). Figure 4d illustrates propagation ranges in for FM cells when n d (max. rate of IP 3 production by PLCd) and r 5P (max. rate of IP 3 degradation by IP-5P) vary. We locate with black dots the (n d ,r 5P ) pairs for which Ca 2+ waves propagate across the whole cell chain. Clearly, long-range propagation is found for a wide region of this parameter space. As expected, larger IP 3 synthesis rates must be balanced by larger IP 3 degradation rate to allow long range propagation, hence the diagonal-like aspect of the black region in the panel.
These first results thus indicate that the propagation distance of Ca 2+ waves in our model is much smaller with linear gap junctions than with nonlinear ones. This observation remains valid when the number of cells in the chain is much larger (we have simulated up to 120 cells in the chains) or/and when up to the 20 first cells in the chain receive the stimulation simultaneously (not shown). The above results are also robust with respect to the changes in boundary conditions (see Methods). For instance, Figure S4 illustrates longdistance Ca 2+ wave propagation for a chain of N = 12 FM-encoding astrocytes with periodic or absorbing boundary conditions and gap junctions endowed with sigmoid-like coupling.
This confirms that the difference of propagation distance between linear and nonlinear gap junction-coupling is a robust and fundamental property of our model. Hence, the existence of a threshold concentration for cell-to-cell IP 3 diffusion, similar to the one displayed by nonlinear gap junctions may be a critical factor for long-distance propagation of Ca 2+ waves across astrocytes. In what follows, we examine the influence of a second physiological characteristic of Ca 2+ signaling in astrocytes, namely their stimulus encoding mode (FM-encoding or AFM-encoding chains).
AFM vs. FM cells. As illustrated in Figure 5, Ca 2+ waves do not propagate in our model of AFM-type astrocyte chains. In Figure 5, large Ca 2+ variations are observed only in the driving cell (cell 1 in Figure 5), whereas the other astrocytes exhibit subthreshold Ca 2+ changes or no Ca 2+ change at all. Importantly, Figure 4b shows that this observation is not restricted to the parameters of This intrinsic difference in the propagation properties brought about by AFM or FM modes can be explained on the basis of the single-cell bifurcation diagrams (Figures S1a, c). Indeed, in the AFM-encoding mode, the peak concentration of the IP3 and Ca 2+ oscillations decreases with decreasing IP 3 stimulations. Hence the IP 3 generation in AFM is such that a local depression of transmitted IP 3 will be accentuated in the next cell. Any decline of IP 3 production in a given cell will thus be transmitted outward and amplified along the chain, until the signal eventually fails. This phenomenon is not observed with FM cells because, by definition of the FM mode, the peak amplitude of the IP 3 oscillations in cell i is hardly dependent of the strength of the IP 3 stimulus coming from cell i21 (at least in the limit where the incoming IP 3 stimulus falls within the oscillatory range of cell i). This simple mechanism guarantees that peak IP 3 values in cell i will remain high even though the incoming stimulation is lower. In other words, the FM mode guarantees robust regeneration of the wave propagation.
Moreover, the range of IP 3 input that gives rise to oscillations in the AFM encoding regime is much narrower than in the FM case. Thus, a perturbation of the IP 3 stimulation from cell i to i+1 in the AFM mode is more likely to be enough to push cell i+1 outside of its oscillatory range, leading to termination of wave propagation in this cell. Importantly, our simulations with AFM-encoding cell chains reported propagation failures for all the F values ( Figure 4c) and (n d ,r 5P ) pairs that were tested (results not shown).
Therefore, these results suggest a neat functional difference between AFM and FM oscillations in astrocytes: while FM could support long distance propagation of pulse-like Ca 2+ waves, AFM is rather expected to give rise to localized Ca 2+ signalling with diffusion-like spatial patterns for IP 3 . Hence, any parameter relevant to intra-cellular Ca 2+ signaling and able to switch the cell between AFM and FM modes (e.g. the affinity or activity of the SERCA pumps) is predicted to play a key role in the inter-cellular propagation of Ca 2+ signals in astrocytes.
Propagation in composite chains. Because the astrocyte population within the brain is heterogeneous [49], the results reported above question the possibility of intercellular Ca 2+ wave propagation across astrocytes with different properties. Here we tackled this issue using composite astrocyte chains, namely chains constituted of both FM and AFM cells, and investigated under what conditions propagation is possible with nonlinear sigmoid gap junctions.
In Figure 6, we stimulated the first cell of the chain (cell 1) with a constant stimulus so as to initiate the Ca 2+ wave in this cell. The intensity of the stimulus was set close to the upper edge of the cell oscillatory range according to the bifurcation diagrams in Figures S1c, d so as to maximize the chance of wave propagation. Figure 6a illustrates the propagation of the Ca 2+ wave in a chain of alternating FM (black traces) and AFM (gray traces) cells and shows that propagation abruptly terminates at the second AFM cell in the chain (cell 4). Notably, closer inspection of IP 3 dynamics in the subsequent FM cell (i.e. cell 5 in Figure 6d) reveals that the IP 3 concentration intermittently passes across the predicted threshold for oscillations in this cell. However the time spent above the threshold is never large enough to trigger CICR, so that propagation halts. Extending propagation to further cells in the chain thus demands faster endogenous IP 3 production. This can for instance be obtained with larger values of the maximal rate of PLCd, n d , and/or smaller IP 3 degradation rates, n 3K and/or r 5P .
The former possibility is considered in Figures 6b,e. The simulation reported in these figures corresponds to the same conditions as in Figures 6a,d, except that n d is larger in AFM cells. Clearly, wave propagation now extends across the entire chain. Due to the increased rate of IP 3 production, all AFM cells in the chain in fact maintain intracellular IP 3 concentration either beyond or within the oscillatory range. Moreover, since this range essentially overlaps with the lower part of the oscillatory range for FM cells, the IP 3 transported from one AFM to the next FM cell in the chain can trigger CICR there, thus perpetuating propagation.
Another possible mechanism to facilitate wave propagation across AFM cells in heterogeneous conditions consists in increasing the frequency of the wave pulses in the FM cell preceding the AFM one. This effect is actually naturally obtained when several successive FM cells are placed between two AFM The presence of homogenous FM cell domains between AFM cells is therefore likely to enable long traveling distances for propagating Ca 2+ waves. One may even assume that if the number of successive FM cells in the FM domains is large enough, the Ca 2+ wave should propagate over the entire network, whatever its size. Although we did not further investigate this possibility, our simulations hint on the contrary that there likely exists an upper bound for the travelling distance because the frequency of propagating Ca 2+ waves in FM domains is not constant, but tend to decrease after each AFM cell, as can be seen by comparison of the pulse frequency in the Ca 2+ traces of FM cells in Figures 6a-c. Indeed such progressive decay of the pulse frequency along the chain eventually brings forth insufficient IP 3 diffusion through gap junctions, thus terminating propagation (results not shown).
We note that modifications of the IP 3 threshold for diffusion in the nonlinear gap junctions, IP 3 thr , should facilitate transmission of Ca 2+ pulses from cell to cell in the chain, thus increasing the frequency of propagation ( Figure S5), with the possibility of observing very different dynamics of propagation in chains of identical astrocytes. Such scenario supports the notion that although nonlinear gap junctions could explain long-distance propagation, their specific properties are expected to be critical factors for the dynamics of propagation. This aspect is further investigated in the next and last section of results of our study.

Propagation of complex waves
The Ca 2+ and IP 3 dynamics observed so far were all obtained using a rather high value of the coupling strength (F = 2.0 mM?s 21 ). In these conditions, the properties of the propagated waves are rather simple: a pulse-like (or not) wave front travels across astrocytes, with conserved shape and either stops after a few cells or invades the whole cell chain. However, our system is a spatially extended dynamical system with large numbers of degrees of freedom. Such systems (e.g. coupled map lattices) are known to manifest complex spatiotemporal behaviors when the coupling strength changes. To get an insight on the possible propagation behavior exhibited by our model with weaker coupling, we considered the dynamics with reduced levels of gap junction permeability (setting F = 0.23 mM?s 21 ). Figure 7 shows the Ca 2+ dynamics of 41 coupled FM cells and a square wave periodic stimulation applied to the central cell #21 (see figure caption for details). Visual inspection of the Ca 2+ traces in each cell (Figure 7a) indicates that such periodic (oscillatory) stimulation can trigger Ca 2+ waves that can propagate along the whole chain. Importantly, this figure also evidences the occurrence of occasional propagation failures that do not seem to result from a simple spatiotemporal pattern. Actually, observation of the temporal traces of each individual cell reveals the occurrence of pulse-like events showing up with no apparent regularity. Accordingly, the distribution of the time-intervals between two such pulses can be very broad for some cells, with large intervals often almost as probable as small ones (see Figure S6). Albeit consistently pulse-like, the shape of the propagated Ca 2+ waves is also quite variable. Closer inspection of the time series for the driving cell (i.e. cell 21) for instance shows that the generated Ca 2+ pulses vary from a single-peak waveform to multiple peaks per single pulse (Figure 7c). Furthermore, Figure 7d shows that the variability and complexity of the IP 3 signals is also very large. The lack of obvious regular behavior is particularly striking on movies showing the parallel temporal evolution of the Ca 2+ and IP 3 level in each cell, as in Video S1 in the Supplementary Information.
To further illustrate the complexity of the obtained dynamics, we plot in Figure 7b the trajectory of the system in the phase space of the driving cell. It is very tempting to compare the resulting trajectories to those observed with classical low-dimensional strange attractors. In this regard, preliminary analysis of the three time series of the driving cell using nonlinear time series analysis tools [46] suggested that the dynamics indeed corresponds to deterministic chaos, with sensitivity to initial conditions testified by a positive maximal Lyapunov exponent that we estimated between 0.020 and 0.050 s 21 (depending on the time series under consideration).
The apparent complexity of the dynamics is most likely due to some form of spatiotemporal chaos, the nature of which is beyond the scope of the current article and is left to future work. But whatever the response, these simulations evidence that complex Ca 2+ wave propagation patterns can manifest at low couplings, even with spatially homogeneous cell properties and in the absence of any stochasticity source.

Discussion
Calcium-mediated signalling is a predominant mode of communication between astrocytes [50]. Consequently, it is important to understand how different biophysical mechanisms determine the ability of these brain cells to communicate over long distances. Here, we used the computational modeling approach to study the properties of gap junction-mediated signaling in simple networks of realistically modeled astrocytes. Using numerical simulations and tools of bifurcation theory, we showed that longdistance regenerative Ca 2+ wave propagation is possible when the gap junctions are rendered by nonlinear permeability but only when most of the model astrocytes are tuned to encode the strength of incoming IP 3 signal into frequency modulated Ca 2+ oscillations.
There has been a long-standing debate over the nature and characteristics of intercellular Ca 2+ waves observed in astrocyte networks. The present article concerns about the purely intracellular route, which involves the transfer of IP 3 molecules directly from cytosol to cytosol through gap junctions [23]. In the extracellular route instead, propagation is mediated by extracellular diffusion of ATP and purinergic receptor activation [18,24]. Although these two routes need not be mutually exclusive, experiments indicated that their relative influences vary across the brain. Indeed, experimental evidence suggests that the purely intracellular route predominates in astrocytes of the neocortex [14,51] and the striatum [28] while the extracellular (purinergicdependent) route seems predominant in the CA1 hippocampus area as well as in the corpus callosum [51]. Hence, the results obtained in the present paper are expected to be relevant to the former structures. Their relevance to the case where the two routes coexist could however be tested by simple extensions of our current model, in the spirit of the recent Ref. [52].
A critical issue for the modeling studies of intercellular Ca 2+ waves is to explain the observed variability of Ca 2+ wave travelling distance [14]. Indeed, experimental measurements show travelling distances varying from 30 cells [16] (that is, often outside the imaging microscope field) down to 3-4 cells only [15]. Models featuring purely regenerative waves (e.g. traveling waves in the usual mathematical sense) easily account for long distance propagations but hardly account for the observed short ones. Conversely, nonregenerative models (e.g. purely diffusive ones) cannot explain long-range propagation. A possible solution was suggested by Höfer et al. [32]. In the model proposed by these investigators, long-range propagating Ca 2+ waves are obtained via IP 3 regeneration in each cell by Ca 2+ -activated PLCd. However, whenever PLCd maximal activity is lower, regeneration becomes partial and the Ca 2+ wave propagation distance decreases. Yet this model does not include Ca 2+ -dependent IP 3 degradation, which could be critical for the occurrence of IP 3 -mediated Ca 2+ oscillations [33]. This latter process in particular, can compete with PLCd-mediated IP 3 production, thus hindering IP 3 regeneration and Ca 2+ wave propagation. This calls for additional factors to be taken into account to explain intercellular Ca 2+ wave propagation.
A first prediction of our model is that, regenerative waves are possible in a network composed in its majority of astrocytes that encode information about incoming IP 3 signals in the frequency of their Ca 2+ oscillations (FM). Interestingly, the response of astrocytes in vivo to IP 3 stimulation is known to exhibit high variability, both in frequency and amplitude [53][54][55]. This variability could be due to cell-to-cell heterogeneity (extrinsic noise) in some of the CICR parameters. In particular, this could include variability of the expression of PLCd or of the affinity for Ca 2+ of the SERCA pumps. To our knowledge, the kinetic properties of SERCA2b have never been measured in astrocytes. However the hypothesis that SERCA2b affinity for Ca 2+ shows variability in vivo seems realistic, given the experimental literature. First, reports of experimental measurements of the SERCA2b affinity showed somewhat variable results, ranging from 170 [56] to 270 nM [57], albeit both studies used cDNA transfection in COS1 cells. Secondly, SERCA2b functionality can be directly modulated by quality-control chaperones of the ER, e.g. calreticulin and calnexin [58]. In particular, there exists strong indication that calreticulin may dynamically switch SERCA affinity for Ca 2+ from 170 to ,400 nM [59]. In this case, cell-to-cell variability in the concentration of calreticulin could result in the mixed AFM-FM cell networks studied here. Our observations then lead to the experimental prediction that such variability or heterogeneity of the astrocyte response would have a strong impact on the propagation of intercellular calcium waves between these cells. Notably, this scenario is also supported by several experimental studies [21,60]. In particular, calreticulin has been shown to regulate Ca 2+ wave propagation via direct interaction with SERCA2b thus modulation of Ca 2+ uptake by this pump [61].
In our model, the strength and the transfer properties of the gap junction coupling are critical permissive factors that allow longrange intercellular signaling between the astrocytes. In particular, nonlinear gap junctions were found to significantly enhance the range of Ca 2+ wave propagation (as opposed to the classic linear gap junctions that caused fast dissipation). Gap junctions with dynamic resistance are known to exist in cardiac networks [62,63] and in several other cells [35]. Yet there is currently no direct evidence for nonlinear transfer of second messenger molecules through gap junctions between astrocytes. Nonetheless, the activation of PKC, which is intimately related to IP 3 metabolism [33,64], is known to block astroglial gap junction communication and inhibit the spread of Ca 2+ waves therein [65]. Hence, in light of the existing knowledge regarding the control of gap junctional permeability by various signaling molecules [66], it is plausible to assume that some nonlinearity should exist in astrocytes too. The exact form of the nonlinearity of course will be dictated by the properties of the solute and the nature of its interaction with the membrane channels in the proximity of the gap junction complex. Meanwhile, the generic form of nonlinear coupling that we considered here allowed us to get a qualitative insight into the putative effect of nonlinear coupling on signal propagation in model astrocyte networks.
In the present study, we considered a simplified setup of 1D network implemented as a regular chain of coupled cells. Such 1D chains display attractive aspects. In particular, we could proceed to a numerical bifurcation study of these 1D coupled-cell systems (see Figure S2), which has proven invaluable for the interpretation of the simulation results. Such bifurcation analysis would hardly be possible in higher dimensions (e.g. 2D), because the number of cells one needs to account for in 2D is much larger than in 1D at constant propagation distance. Furthermore, a serious study of a 2D system must include the exploration of the influence of the coupling network topology [67], which adds further parameters to the study of the robustness of the model dynamical features. However, real astrocytes in tissues are believed to organize in quasi 2D networks with significantly more complex structure. Our model is thus a simplification of this quasi 2D reality. For instance, obstruction of wave propagation could dependent on the spatial dimension. Indeed in 2D or 3D reaction-diffusion systems or on random graphs, where the strength of the coupling or the local number of neighbors can vary across the network, the wave propagation distance can critically depend on the number of stimulated cells or the distribution of the number of coupled neighbors [68]. It is not yet clear whether our observation that linear gap junctions support only local wave propagation is restricted to regular 1D networks such as those used in the present work. Future works will be designed to tackle this issue. Nevertheless, in spite of its simplicity, this 1D model yields important predictions about the influence of the spatial arrangement of astrocytes. In particular, it shows that the distribution in space of heterogeneous gap junction permeabilities can result in rich dynamics [14,30]. Reducing the maximal strength of coupling between the model astrocytes imparted the individual cells with rich dynamics, possibly associated with spatiotemporal chaos. Keeping in mind that in reality the changes in gap junction permeability are mediated by the dynamic action of different effectors, we anticipate that a network of biological astrocytes could have the capacity to self-regulate the complexity of its dynamics. Whether or not this is the case, can be determined by experiments that selectively target the pathways of gap junction regulation.
Recent studies suggest that the astrocytes within the cortex form heterogeneous populations [69,70]. Therefore, we considered the case of intercellular Ca 2+ wave propagation in composite 1D networks, consisting of both FM-and AFM-encoding cells. Our simulations predict that the propagation dynamics and distance of intercellular Ca 2+ waves critically depends both on the encoding property of the cells and on their spatial arrangement. Interestingly, the cell bodies of neighboring astrocytes within the brain are believed to distribute in space in a nonrandom orderly fashion called ''contact spacing'' [71,72]. Our study thus suggests a possible link between contact spacing and intercellular Ca 2+ wave propagation in astrocyte networks. If, as suggested by our model, the spatial arrangement of the astrocytes, coupled to the heterogeneity of their response, conditions Ca 2+ wave propagation, then contact spacing may play a critical role in intercellular wave propagations in the brain and the related computational properties of astrocyte networks.
It is now widely accepted that astrocytes and neurons are interwoven into complex networks and are engaged in an intricate dialogue, exchanging information on molecular level [3]. By releasing different gliotransmitters (such as glutamate and ATP) astrocytes dynamically modulate the excitability of neurons and control the flow of information at synaptic terminals [4]. Diffusion of glutamate and/or ATP is limited due to the action of glutamate transporters and degradation of ATP, thus defining spatiotemporal range for the local effect of astrocyte on neurons and synapses [73]. On the other hand, long-range and temporally delayed regulation of neuronal and synaptic activity by astrocytes is believed to be mediated by intercellular Ca 2+ waves spreading through the astrocyte network. The connectivity of this astrocyte network is in turn defined by the patterns of electrical activity in neuronal network [74]. Thus, it appears that astrocytes and neurons are organized in networks that operate on distinct time scales and utilize the principles of feedback regulation to modulate the activities of each other. How such mutual regulation of neuronal and astrocytic networks affects the complexity of neuronal network dynamics in health and disease is a question that should be addressed by future combined experimental and modeling studies. Figure S1 Bifurcation analysis of an uncoupled (i.e. isolated) ChI astrocyte for AFM (a, c, e) and FM (b, d, f) encoding regimes. (a, b) 3D-rendering of bifurcation surfaces in the state space. AFM oscillation amplitude (c) and period (e) are controlled by a supercritical Hopf (H) bifurcation and a saddle-node limit cycle (SNC) bifurcation respectively. Conversely in FM-mode, the occurrence of a saddle-node on an invariant circle (SNIC) bifurcation accounts for the rise of arbitrarily-small frequency Ca 2+ oscillations (f) at almost constant amplitude (d). Legend: (a, b): black lines: stable fixed points; red dashed lines unstable fixed points; blue lines: bifurcating limit cycles; semi-transparent surfaces denote envelopes of stable (grey) and unstable (red) oscillations. (cf): green: IP3; orange: Ca 2+ (c) full lines: stable oscillations; dashed lines: unstable oscillations. Parameters as in Table 1. Although not apparent in the figure, for IP 3 bias values larger than <0.8 mM, the stable oscillations become far more complex than in the isolated case. This is due to a very rapid cascade of period-doubling bifurcations, which yields extremely complex limit cycles (with numerous folds) that could not be precisely rendered in the figure (see also Section III.1.b). Moreover, for IP 3 bias .1.1 mM, the amplitude of these limit cycles shrinks and numerical investigations evidenced the coexistence of multiple complex stable orbits. Legend: thin full lines locate unstable fixed points, and thick full lines stable one. Full (open) circles denote the envelopes of stable (unstable) limit cycles. Letters denote bifurcation type as in Figure SI1. Found at: doi:10.1371/journal.pcbi.1000909.s002 (1.66 MB TIF) Figure S3 IP 3 -triggered CICR-mediated propagation of a pulsed Ca 2+ wave within a chain of five FM ChI astrocytes (A1-A5). (a) An IP 3 stimulation of constant intensity (IP 3 bias = 0.8 mM) is applied to cell A1 from t = 10 s to t = 30 s. This increases IP 3 concentration, thus triggering CICR from the ER and the generation of a Ca 2+ pulse. (b) By means of communication through gap junctions, suprathreshold IP 3 from A1 can diffuse to A2, triggering CICR there. The process is essentially regenerative so that a Ca 2+ pulse almost identical to the original one can be observed in the arrival cells. (c) As soon as the IP 3 influx to one cell from its neighbors is not sufficient to trigger CICR, the propagation stops. This is indeed the case of cells A4 and A5. Cells were coupled by sigmoid gap junctions and experienced reflective boundary conditions. Video S1 Illustration of the chaotic-like behavior at low coupling. This movie shows the evolution with time of the calcium (upper panel, blue bars) and IP 3 (lower panel, green bars) concentration in each cell (x-axis) of a N = 24 cell chain. The total duration represents 540 seconds of real time. Parameter values are the same as in Figure 7 (in particular, F = 0.23 mM?s2 1 ). Stimulation triggered by IP 3 bias = 1.0 mM from t = 0 s to t = 540 s applied to the first cell.