A phenomenological cluster-based model of Ca2+ waves and oscillations for Inositol 1,4,5-trisphosphate receptor (IP3R) channels
aa r X i v : . [ q - b i o . S C ] S e p A phenomenological cluster-based model of Ca waves and oscillations for Inositol1,4,5-trisphosphate receptor (IP R) channels
Svitlana Braichenko ∗ and Atul Bhaskar Faculty of Engineering and the Environment, University of Southampton, SO17 1BJ, Southampton, UK
Srinandan Dasmahapatra
School of Electronics and Computer Science, University of Southampton, SO17 1BJ, Southampton, UK (Dated: September 14, 2018)Clusters of IP receptor channels in the membranes of the endoplasmic reticulum (ER) of manynon-excitable cells release calcium ions in a cooperative manner giving rise to dynamical patternssuch as Ca puffs, waves, and oscillations that occur on multiple spatial and temporal scales.We introduce a minimal yet descriptive reaction-diffusion model of IP receptors for a saturatingconcentration of IP using a principled reduction of a detailed Markov chain description of individualchannels. A dynamical systems analysis reveals the possibility of excitable, bistable and oscillatorydynamics of this model that correspond to three types of observed patterns of calcium release – puffs,waves, and oscillations respectively. We explain the emergence of these patterns via a bifurcationanalysis of a coupled two-cluster model, compute the phase diagram and quantify the speed of thewaves and period of oscillations in terms of system parameters. We connect the termination oflarge-scale Ca release events to IP unbinding or stochasticity. PACS numbers: 82.20.Wt, 82.40.Bj, 82.40.Ck
I. INTRODUCTION
Ionic calcium (Ca ) is a second messenger involvedin many processes such as fertilization, cell prolifera-tion, differentiation, embryonic development, secretion,muscular contraction, immune response, brain functions,chemical sensing, light transduction, etc. [1, 2]. Whilemany types of Ca release channels are implicated inthese processes, in this paper we are concerned with thedynamics of Ca released into the cytosol from inositol1,4,5-triphosphate receptor (IP R) channels [2] that aretriggered by IP secondary messenger molecules. IP Rchannels are usually present in non-excitable cells andare mainly localized in the endoplasmic reticulum (ER)membrane where they are believed to form tight clus-ters to enable Ca signaling. Experimental observa-tions [3, 4] suggest that the distribution of clusters variesby cell type. For example, SH-SY5Y neuroblastoma cellscontain around 4 to 6 channels per cluster (a few con-tain up to 10 channels), HeLa cells clusters on averagecontain 2 to 3 channels and Xenopus oocytes clusterscontain around 20 channels each [5]. However, theoreti-cal and experimental investigations of clustering and themechanisms by which it regulates Ca signaling are in-complete and is the subject of ongoing research.Ca signaling events corresponding to release fromCa channels [1] constitute a multiscale hierarchy con-sisting of three distinct types of events taking place atdifferent spatial and temporal scales involving individualchannels, clusters of channels and groups of channel clus-ters. The events occurring at the smallest scale ( ∼
10 nm ∗ Corresponding author; [email protected] and ∼
10 ms) are called blips . They are the buildingblocks of signals on larger scales. Blips occur when Ca is released from a single channel into the cytosol. Co-ordinated Ca release from a cluster of channels – acollection of co-occurring blips – is called a puff and ap-pears at the second level of the hierarchy ( ∼
100 nm and ∼
100 ms). This cooperation between signaling eventsis regulated by Ca in a concentration-dependent feed-back mechanism. Low Ca concentrations (maximumactivity for 200 −
500 nM) diffuse around the cytosol andcause further release from neighboring channels, whereasCa in high cytosolic concentrations inhibits furtherCa release. This feedback mechanism is called Ca induced Ca release (CICR). The final level in the hier-archy ( ∼ µ m and ∼ spikes , waves , and oscillations . Spikes correspond to local orglobal transient releases from a group of clusters. Ca waves are formed as sequential releases travel from ex-cited clusters to neighboring ones. Ca oscillations arerepetitive Ca spikes or waves. All the events on thehighest hierarchical level that we consider (spikes, waves,oscillations) consist of the coordinated release of Ca puffs from a group of clusters.Ca oscillations are experimentally observed in var-ious types of cells [6]. Many models associate the oc-currence of Ca oscillations with the IP regulation [6–10]. Some studies connect the emergence of Ca os-cillations with the oscillating level of IP triggered bybinding to various enzymes, such as PLC (phospholipaseC) or IP kinase [6, 9]. Others report that Ca oscil-lations may occur even if the IP concentration is non-oscillatory [7, 8]. Various studies explain oscillatory Ca behavior as a sequence of stochastic spikes [11–13]. Apartfrom oscillations and spikes propagating wavefronts havealso been much studied. These are implicated, for in-stance, in oocyte fertilization, and are generated when asperm cell makes contact with an oocyte.There is an extensive literature on the deterministicmodeling of Ca front propagation. References [7, 8, 14–17] consider the averaged Ca dynamics from a well-mixed ER membrane without addressing the clusteringof channels. Channel clusters, including inhomogeneouscluster distributions [18], are considered in determinis-tic [19–21] and stochastic models [19, 22] of Ca wavepropagation. The current paper studies the determinis-tic properties of Ca waves and oscillations occurringin the ER membrane containing multiple clustered chan-nels. The models in [20, 21] inject a burst of Ca ionsat each cluster which diffuses between clusters. Whilethe burst size is represented by a parameter in theirmodels, we choose to represent the release mechanismby a reduced version of the DYK [7] model proposedby R¨udiger [23]. This gives our model access to sin-gle channel characteristics [7] and cluster properties [23]providing more mechanistic insight into the link betweensingle-channel properties and features of waves and os-cillations at higher levels in the spatial hierarchy. Weuse dynamical systems analysis to explore the differentdynamical regimes at the single-cluster level and relatethese properties to larger spatial scales. Thus, unlike thefire-diffuse-fire model [20, 21] we can study the charac-teristics of global Ca release events such as waves andoscillations, including their inhibition and termination,using channel and cluster parameters. We also build onthe results of a detailed hybrid stochastic-deterministicsingle-cluster study [24] and its extension on a grid ofclusters [25] that show different durations of Ca eventsdepending on the IP level.In detail, we introduce a spatial model of Ca wavesand oscillations building upon qualitative insights fromdynamical systems theory applied to reduced, few vari-able models [8, 23] of IP R channel subunits [7] inte-grated into channels and heterogeneous clusters of chan-nels [24]. We assume the IP concentration to be sus-tained at a high level so that it increases the probabil-ity of the occurrence of Ca waves and oscillations, ex-tending the observation of excitability of channel clustersin [23]. Further, we account for variability of [IP ] lev-els by including the IP unbinding. The extended modelalso gives insight on the numbers of activatable chan-nels (all four subunits has IP bound [26]) in agreementwith [24, 25]. The model we propose is of a deterministicreaction-diffusion type where diffusion of Ca betweenclusters of channels smooths the Ca releases from indi-vidual activated clusters modeled by reaction dynamics.Here we study the dependence of the velocity and periodof waves and oscillations on the diffusion coefficient, equi-libration rate and a number of channels in each cluster.We observe a curious spatially alternating wave patternin a regime of parameter space that is reminiscent of theoccurrence of intracellular spatially organized Ca al-ternans [27]. The termination of large-scale events suchas waves and oscillations has been a topic of considerable discussion [28–32]. We provide evidence for two differentmechanisms – a stochastic termination and a determin-istic IP unbinding – that can be accommodated in ourmodels in order to account for the decay of Ca excita-tions in simulation traces. II. METHODS
A kinetic model of channel activation and deactivationis the basis for models of Ca dynamics, and many mod-els [8, 23] have been proposed that reduce the detailedkinetic scheme laid out by De Young and Keizer – the so-called DYK model [7] as a starting point. Our approachto simplifying the complexity of the problem is based onseveral basic models briefly discussed in this section. Theassumption of a well-mixed membrane where spatial ef-fects are ignored is a shared limitation of these models[7, 8]. The spatial extension of the three-state model [23]is proposed in the last subsection. For further detailsplease refer to the Appendices and Supporting Material. A. DYK and derived models of Ca signaling The DYK model [7] is one of the first models whichsuccessfully accounts for key empirical findings of IP Rchannels. In the model each channel consists of four sub-units, at least three of which need to be activated forchannel opening, reflecting the three conductance lev-els experimentally observed in [33]. Each subunit has3 binary state variables i, j, k ∈ { , } , correspondingto 3 binding sites – for IP , activating Ca and in-hibiting Ca – with a state transition diagram shownin Fig. 1 (a) and described in its caption. The frac-tion of subunits in a given state ( ijk ) is represented bythe variable x ijk , such that they sum up to unity, i.e. P i =0 P j =0 P k =0 x i,j,k = 1. The presence of IP at eachchannel enhances Ca release, while calcium either am-plifies or inhibits further Ca release (called calcium-induced calcium release, CICR).The DYK model assumes that all 4-subunit channelsare “well-mixed” in the membrane and dispenses withspatial effects such as clustering. The complexity of thestate space of the model makes it computationally dif-ficult to obtain properties of the spatial distribution ofchannel clusters in the ER membrane. Hence, multipleattempts at reducing its complexity have been made inthe last few decades, the majority of which share thecommon feature of IP -dependent activation and CICR.In one such reduced model [8], the variables corre-sponding to the DYK cube [Fig. 1 (a)] are dispensed withby assuming that activation by IP and Ca are fast,compared with the slow inhibition by Ca . The stateindices ( i, j,
0) and ( i, j,
1) distinguish the two states ofthe site where the inhibiting Ca binds. The fractions x ij and x ij of states with variable IP and activatingCa binding sites on the left and the right faces of the (d) Clusterlevel(b) Subunitlevel IP Ca act Ca inh (c) k i- k a + c k a - k a + c k a - k i- az gh' k i + ~ k i + c h k a + c k a z k k i + c k a - k i + ~ (a) b a c active a p a p a p a p b b b
110 111 b b a c a c b a c b a c b a c b a c b a c b R I GH TLE F T U PP E R a c b a c b a c a c b b U PP E R FIG. 1: (Color online) (a) DYK cube.
Schematicdiagram of the DYK eight-state gating model of theactivation of a single subunit of an IP R channel; theeight corners of the schematic will be referred as theDYK cube [7]. Each vertex of the cube directlycorresponds to a specific state of a subunit labeled withnumbers ijk , where i corresponds to the IP activatingsite, j to the Ca activating site and k to the Ca inhibiting site. Each binding site can be either unbound( i, j, k = 0), or bound ( i, j, k = 1). The green arrow indicates binding (unbinding) of IP to the binding sitewhile red dashed and blue dot-dashed arrows correspondto Ca activating and inhibiting transitionsrespectively. The transitions between the states aregoverned by the second-order ( a l p or a l c , where p =[IP ], c = [Ca ]), and the first-order ( b l ) rateconstants, where l = 1 , (b) Upper DYK plane. When p is high, the most of the transitions take place inthe upper plane of the cube. (c) Four-state model. The scheme of the model [23] for high p at the level of acluster of channels. The fractions of the states in acluster are marked as a, g, z, h ′ instead of110 , , ,
101 for subunits. The transition ratesbetween channel states are k ± a,i and ˜ k + i = k + i c s toindicate the elimination of detailed balance in a clusterof channels [34]. (d) Three-state model. The finalmodel includes a compound state h formed from thestates g and h ′ with the effective rates k , .DYK cube [Fig. 1 (a)] are slow variables in the slow-fast reduction. The pairs of fractions ( x j , x j ) and( x j , x j ) are determined from the steady state of thefast IP binding and unbinding processes. The dynam-ics of the reduced model is governed by the compounddynamics of the left-right faces and Ca . B. Single-cluster models
Another reduction of the dynamics of a group of DYKsubunits is introduced in [23] where, in contrast to [8],dynamical variables capture the collective behavior of acluster of IP R channels. This model shows that for largevalues of IP the dynamics is excitable, and explains theoccurrence of Ca puffs. As this model is the start-ing point of our analysis we present the reduction of theIP dynamics explicitly in Appendix A for completeness.Here we discuss its fundamental assumptions. The re-duction scheme follows a similar slow-fast approach [8]but for the case of high p = [IP ]. This assumption leadsto fast transitions with rates a , p [Fig. 1 (a)] and allchannel subunits are forced to have the IP binding siteoccupied leaving only 1 jk states to be dynamical (theupper plane of the DYK cube in Fig. 1 (a)). The dynam-ical variables of the four subunits shown in Fig. 1 (b) aregathered together to introduce cluster variables a , g , h ′ ,and z . These represent the fractions of channels withina cluster in the corresponding states shown in Fig. 1 (c).The transition rates k − a , k ± i are derived from the param-eters in the DYK model [7]. The transition rate k + a isobtained from the requirement that at least three sub-units need to be activated to open a channel (A28). Thenumerical values of the coefficients [23] are provided inTable I. A further simplification – corresponding to thepassage from Fig. 1 (c) to (d) – is introduced [23] by re-placing the pair of states ( g, h ′ ) representing IP -bound,inactivating Ca -bound states with a compound state h whose kinetics has effective rates k and k [(A32),(A33)]. This reduces the cluster dynamics to a two-statemodel as the subunit fractions within a cluster sum tounity, z = 1 − a − g − h ′ = 1 − a − h . A significant mod-eling assumption is made whereby the rates e k + i and k + i c in Fig. 1 (c) are taken to be different [34], in contrast tothe equal rates a c on the upper plane of Fig. 1 (a) madein [7], a condition that implies detailed balance.The three-state model of a cluster of IP channelsthus obtained exhibits excitable and bistable dynamicalregimes [23]. In our study, we modify the model slightly,by setting the parameters for a single cluster of channelsto fit the open probability distribution of IP channels inthe presence of Ca in Xenopus oocytes [24]. This leadsto changes in the transition rates a and b as shown inTab. I. Consequently, an oscillatory regime emerges in ad-dition to the excitable and bistable regimes found in [23].Later in this paper, we spatially extend this model in acluster-based manner to study the interplay between thedifferent regimes of clusters in producing wave and oscil-lation dynamics in a membrane.The dynamics of a cluster is governed by the set ofdynamical variables x ( t ) = { a ( t ), h ( t ), c ( t ) } correspond-ing to the fraction in the cluster of opened channels a ( t ),inhibited channels h ( t ), and the cytosolic Ca concen-tration c ( t ). The dynamical variables in the phase space,as shown in Fig. 1 (d), are governed by a system of cou-pled ordinary differential equations [23] ˙ x = F ( x ( t )) withappropriate initial conditions:d a d t = k + a ( c ) c (1 − a − h ) − k − a a + k ( c ) h − ˜ k + i a = △ f ( a, h, c ) , (1)d h d t = k + i c (1 − a − h ) − k ( c ) h − k ( c ) h + ˜ k + i a = △ g ( a, h, c ) . (2)The rates k a ( c ), k ( c ) and k ( c ) are given in eqs. (A28),(A32) and (A33), respectively [23, 32]. It is assumedin [34] that after a blip Ca levels at a channel poreremain high while neighboring channels in a cluster relaxto levels around c d . According to [35] [Ca ] at a nan-odomain drops rapidly from several hundreds of µ M to afew µ M. Thus, at the timescales of the dynamics studiedin [23] and in our paper, the calcium-dependent inacti-vation rate is taken to be frozen at ˜ k + i = k + i c s as shownin the transition scheme Fig. 1 (c), (d). The dynamics ofthe c variable accounts for homeostatic mechanisms suchas pumps and buffers in the simplest way, whereby Ca levels relax back to a steady-state level c d ( a ) at rate λ :d c d t = − λ [ c − c d ( a )] , (3)where λ = 10 s − [23] sets the fast time-scale so thatTABLE I: Parameters used in the model in [23] and thesingle cluster model of this paper. Parameter R¨udiger [23] Present studyCa activation bindingrate a , µ M − × s −
100 100Ca activation unbindingrate ( b , k − a ) a , s −
20 20Ca inhibition bindingrate ( a , k + i ) a , µ M − × s − . . a Ca inhibition unbindingrate ( b , k − i ) a , s − . . b Local [Ca ] at openedchannel ˜ k + i = k + i c s , s −
30 6 b Rest level of [Ca ] c , µ M 0 .
025 0 . − . α, µ M 0 .
74 0 . N λ, s − Characteristic of the step ǫ . . a Here we show the transition rates for subunit and channel asshown in Figs. 1 (b) and (c) respectively. These rates areassumed to be the same in this study. b These values are determined from patch-clamp data on theprobability distribution for channel opening typical for
Xenopus oocytes at high [IP ] [24] c d ( a ) is slaved to the dynamics of a ( t ). The activationbarrier for a cluster of channels is modeled separatelyby setting a threshold below which the cytosolic level ofcalcium is c and above which the concentration of Ca is approximately linear in the number of opened channels N a : c d ( a ) = c + 0 . αN a (cid:8) N a − /ǫ ] (cid:9) , (4)where ǫ sets the scale for smoothing the transition fromzeroth to first order dependence on opening and α is aconstant defining the strength of the coupling betweenchannels in a cluster. Such a step is proposed in [23] toavoid Ca release from the inactive clusters where thenumber of opened channels N a <
1. This modeling stepis necessary because a and h denoting fractions of openedor inhibited channels are continuous variables.The assumption of fast equilibration to steady-statelevel c d in (3) flattens out c ( t ) from its blip amplitude c s of cytosolic Ca at a cluster. For λ = 10 s − [23],it takes ∼ −
10 ms for Ca levels to reach cytosolicconcentrations in the 10 − nM range. In what follows,we introduce a diffusion term to enable Ca to activatequiescent clusters by CICR thereby having blips combineto form puffs as shown in Fig. 2 (b). This is reminiscentof the fire-diffuse-fire model [21].For a single cluster in a small domain of the mem-brane a computational study [24] using a hybrid reaction-diffusion model – a Markov chain description as perDYK [7] and partial differential equations for diffusingCa ions – shows the dependence of the transition frompuffs to waves on the numbers of channels in a clusterwhich initiate Ca release events under different IP loads. The complexity of the analysis prompts us to lookfor a simplified model that could probe the temporal pat-terns observed and provide qualitative explanations fortheir origin via dynamical systems theory. We extendthis study to a larger domain containing multiple clusterseach described by (1)–(2) [23]. While retaining the struc-ture of the model given in [23] we chose the parametersof [24] in order to study the effect of different numbers ofchannels in each cluster and the interplay between them.This allows us not only to study the transition from puffsto waves but to obtain the characteristics of waves andoscillations that emerge.This is the justification of applying the principles ofslow-fast reduction to the single channel and single clus-ter models – we can elucidate how these reduced modelsmay be coupled to account for physiological behavior thatis manifest spatially. C. Cluster-based spatial model
In this section, we lay down our assumptions for build-ing a spatially dependent model for a phenomenologi-cal representation of a Ca wave as a set of sequentialCa release from clusters [Fig. 2 (c)] using features ofearlier models [7, 8, 15, 21, 23, 24]. An examination of t = 1 s E m e b a e cytosol (a) BLIP t = 10 ms m
FIG. 2: (Color online) Schematic representation of thespatial and temporal hierarchy of Ca release events. (a) A blip as the release from a single channel. The green plane represents a part of the ER membrane, the red arrow indicates Ca release from the openedchannel indicated by the dotted circle . (b) A puff isrepresented by multiple blips occurring in a singlecluster of coupled channels. The coupling betweenchannels (marked by black arrows ) is defined by theconstant α in the nonlinear dependence (4). (c) Thewave propagation as a sequence of puffs ( red lines )caused by raised levels of [Ca ]. In our modelingapproach, we consider clusters of channels to belocalized at single points in space. The activation ofneighboring clusters is caused by Ca diffusion fromthe initial Ca release localized in the center ( blackdotted lines ).the dynamical regimes of individual cluster-level reduc-tions that account for the blip to puff phenomena enablesus to propose a model of diffusively coupled clusters inthe ER membrane where Ca release occurs. We ana-lyze the clusters separately using dynamical systems the-ory. We confirm the existence of two qualitatively differ-ent Ca behaviors in clusters with different numbers ofchannels, viz. , excitable (puff regime) and bistable (waveregime) [23]. We shall probe whether introducing diffu-sive coupling and linear relaxation (3) (that can subsumeeffects including buffering [16] or other modes of reduc-ing spatial gradients) can facilitate the occurrence of os-cillations in Ca release in the ER membrane. This isincorporated into the model by modifying the equilibra-tion eq. (3) by introducing a one-dimensional diffusive transport term ∂c∂t = D ∂ c∂x − λ X i H ( r cl − | x − x i | ) (cid:0) c − c d ( a i ) (cid:1) , (5)where c d ( a i ) is given by eq. (4), D is an effective diffusionconstant, x is a spatial coordinate, H ( x ) is the Heavisidestep-function ensuring that calcium releases into the cy-toplasm occur only at the clusters located at x i , i = 1 , L ,and r cl is the cluster radius. H ( x ) = ( , if x ≥ , , otherwise . (6)For the convenience of notation, we shall assume that the i dependence on a i is implicit hereafter.The model proposed in the current paper is of the fire-diffuse-fire type studied before [21]. Unlike previous stud-ies, we propose a direct link between the Ca level c andthe proportion of opened channels a through nonlineareq. (4). This also helps us to associate cellular behaviorof Ca releases on a spatial scale larger than the clustersize r cl to microscopic channel characteristics [7].System (1)–(2) together with eq. (5) is a reaction-diffusion system. The first term in the eq. (5) correspondsto the smoothing of a Ca front and the spreading of[Ca ] from the active clusters to the neighboring onesas sketched in Fig. 2 (c). In the second, reaction term, c d (4) is a function of the number of opened channels andaccounts for the opening of IP R channels in the clus-ters. The equilibration of cytosolic Ca to c d at rate λ is a linear homeostatic reaction term that subsumesthe action of pumps, leaks, and buffers, which drives thecytosolic [Ca ] to the steady state value c d .We explore next how this model is able to generatewaves and oscillations and characterize these phenomenain terms of the model parameters. III. RESULTS
In subsection III A, we apply the three-state model toa case of high [IP ] and corroborate the existence of atransition from excitable to bistable behavior in a clus-ter of channels [23]. By continuously changing the restCa concentration parameter c in a cluster with noopened channels, we uncover the existence of two Hopfbifurcations. This illustrates the possibility that a sin-gle cluster could exhibit oscillations and channel activityand homeostatic mechanisms could regulate the state ofthe system to undergo dynamical regime shifts.Instead of letting c be a fixed external control param-eter that sets the level of Ca ions in the appropriaterange for oscillatory behavior in the bifurcation analysis,we allow Ca levels to be dynamically set by diffusion,with the effective diffusion constant D being the controlparameter that sets the Ca levels locally and affectingthe physiological output. Diffusion levels gradients; inorder to ascertain whether altering the levels between acluster primed for bistability and that which is excitablymonostable, we set up in subsection III B a two-clustermodel to study its phase diagram by bifurcation analysis.Once again, we find the existence of two Hopf bifurcationpoints as D is altered and obtain the corresponding limitcycle trajectories. We associate those trajectories withthe emergence of Ca oscillations.In the subsection III C, we apply our cluster-basedmodel to study Ca release from a domain in the ERmembrane with clusters containing different numbers ofchannels. We demonstrate how the interplay betweenexcitable and bistable clusters is a mechanism for theemergence of Ca oscillations within the ER membrane.We also show that the cluster-based model is capable ofexhibiting Ca waves that propagate throughout themembrane.In the subsection III D, we show the role of the IP unbinding in the termination of the Ca release from abistable cluster. We account for the transitions betweenthe upper and lower planes of the DYK cube for varying[IP ] levels. A. Ca dynamics at the scale of a single cluster We study the three-state model in eqs. (1), (2) with pa-rameters as shown in the last column of Table I. We shalltreat as control parameters for our bifurcation analysis N , the number of channels in a cluster and c , the restconcentration of Ca ions when none of the channelsare opened, with (0 . ≤ c ≤ . µ M. The choice oftransition rates k + i = 0 .
02 ( µ M · s) − and k − i = 1 .
56 s − were determined by patch-clamp data on the probabil-ity distribution for channel opening typical for Xenopus oocytes at high [IP ] [24]. This has implications on thenumber of channels to be opened within each cluster inshaping the model building and results below.We assume that cytosolic [Ca ] equilibrates to thesteady-state value c d very fast ( λ = 10 s − ) in eq. (3).Therefore, we are able to represent the model of a sin-gle cluster dynamics by eqs. (1), (2) whose solutions arenumerically obtained by standard Runge-Kutta methodsand shown in the ( a , h ) phase plane in Figs. 3 (a), (b),and (c). In the next subsections, this value will be re-duced and the corresponding fast-slow decoupling will nolonger be valid.Figs. 3 (a) and (b) show system behavior as determinedby the orientation and position of the ˙ a = f ( a, h ) = 0( green dashed line ) and ˙ h = g ( a, h ) = 0 ( red solid line )nullclines and we mark by ‘X’ the locations of the stablefixed points where the nullclines intersect. The brokenarrow represents a perturbation of the system away fromthe stable fixed point closest to the origin. After exci-tation, the system evolves following the trajectory whichis guided by the nullclines ( green dashed and red solidlines ). If the initial condition for the dynamical vari-ables in Fig. 3 (a) lies to the right of the middle green (c) (a) N = 5, c = 0.025 (cid:1) M a h trajectorya-nullclineh-nullcline (b) N = 6, c = 0.025 (cid:0) M a (cid:2) ((cid:3)(cid:4) N = 5, c = 0.16 (cid:5) M h (cid:6) a a ah[Ca ], (cid:7) M N = 5, c = 0.16 (cid:8) M t(cid:9)(cid:10) FIG. 3: (Color online) The different regimes of Ca release in the two-dimensional three-state model (1)–(2)assuming c = c d given by (4) [23]. (a) The phasediagram of the excitable cluster ( N = 5), the red solidline is h -nullcline for inhibition, the green dashed line is a -nullcline for opened channel with at least three activeCa bound, the black solid line is a trajectory startedfrom a = 0 . h = 0 . black arrows ) and finished at thestable fixed point ( black cross ). (b) The phase diagramof the bistable cluster ( N = 6), the diagram depicts twostable fixed points. The trajectory starts at a = 0 . h = 0 . (c) Oscillations in the three-state model observed for N = 5and c = 0 . µ M. The only unstable fixed point ismarked as the open circle in the zoomed figure. (d)
The temporal behavior of oscillations for a fraction ofchannels in open state denoted a (red solid) , a fractionof closed channels denoted h (blue dot-dashed) , andcytosolic [Ca ] denoted c (green dashed) . dashed line , the state variables return to the (only) fixedpoint following the extended trajectory in the phase planeshown in black . In Fig. 3 (a) the trajectory comes back tothe only fixed point and demonstrates that the monos-table state is an excitable one, while in Fig. 3 (b) thesystem settles into the fixed point corresponding to thehigher values of ( a, h ), a bistable state. The excitabledynamics creates a short-lasting puff when the Ca con-centration is pulse-like. In contrast, the bistable clusterdynamics ensures that the Ca concentration does notreturn to a base level but stays elevated for a longer pe-riod of time, the behavior associated with a long-lastingpuff. The termination of such a long-lasting puff is at-tributed to IP unbinding which is discussed in subsec-tion III D. Also, we shall show later the possible role ofthe intra-cluster, inter-channel coupling represented by α in eq. (4) in accounting for termination of puffs.We confirm the transition from monostable excitableand bistable behavior [23] in Figs. 3 (a) and (b). How-ever, our use of parameters adapted from [24] revealsthe bistable state occurs when the number of channels is N ≥ N ≥ N ∼ , thehybrid reaction-diffusion simulation [24]. Having con-firmed that changing the total number of channels thatmay be opened alters the dynamics, we then analyzethe stability of the lower fixed point in the phase planeupon changes of parameter c . As seen in Figs. 3 (a)and (b), the rest concentration of Ca ions is takenas c = 0 . µ M. We change c in a continuous man-ner (0 . ≤ c ≤ . µ M to probe qualitative shifts insystem dynamics. The bifurcation diagram is shown inFig. 4 (a) and the real and imaginary parts of the eigen-values ( λ , ) of the Jacobian of the system in Figs. 4 (c)and (d). We find two Hopf bifurcation points ( ℜ ( λ , λ ) = 0 and ℑ ( λ , λ ) = 0) in the range of c that deter-mine the onset or disappearance of limit cycles. c is an external control parameter in the bifurcationanalysis with values in a physiologically plausible range.It is the cytosolic [Ca ] at a cluster before a Ca release event commences when no channels are opened.The oscillatory behavior of the dynamical variables ofthe cluster with initial condition c = 0 . µ M for theconcentration of Ca is shown in Figs. 3 (c) and (d) intwo different representations.To study the stability of the limit cycles emerging fromthe Hopf bifurcation we perform a detailed dynamicalsystems analysis, presenting the results of the continua-tion and bifurcation analysis performed using the Mat-Cont [36] package for Matlab in Fig. 4. This analysiskeeps track of how a continuous change in parameterssuch as c alters the fixed point shown in Fig. 3 (a)and, thus, the eigenvalues of the Jacobian around thealtered fixed point change as well (see, e.g., Fig. 3 (c)).In particular, we calculate the first Lyapunov coefficient l (see [37]) derived from the normal form of the sys-tem of eqs. (1), (2), at each Hopf point (denoted H U and H L ), assuming c = c d as in eq. (4). If l is positive(negative) the Hopf bifurcation is subcritical (supercriti-cal) and the limit cycle is unstable (stable). In our case, l = − . × at H U , the upper point in Fig. 4 (b) (in-set) where a stable limit cycle emerges via a supercriticalHopf bifurcation and l = 5 . × at the lower pointH L . We start the continuation analysis from a maximumvalue of c = 0 . µ M where the fixed point is stable [top ofthe blue curve in Fig. 4 (a)] as shown by the negative realand zero imaginary parts of the eigenvalues of Jacobian[right hand sides of Figs. 4 (c) and (d)]. Upon decreas-ing c we reach the H U ( c = 0 . µ M) point where thestable limit cycle emerges; we observe limit cycles withlarger amplitudes with the further decrease in c ( bluecyclic trajectories in Fig 4 (b)). The cyclic trajectory (c) a H U H L LPC (a) H L LPC H U ℜℜ ℜ H L H U LPC (cid:11)(cid:12)(cid:13) ℐℑℒ ℓ ℐℑℒ ℓ ℐℑℒ ℓ s(cid:14)(cid:15)(cid:16)(cid:17)(cid:18) (cid:19)(cid:20)u(cid:21)(cid:22)(cid:23)(cid:24)(cid:25)(cid:26)(cid:27) (cid:28)(cid:29) H L LPC a h FIG. 4: (Color online) (a)
Bifurcation diagram of thetwo-dimensional system at N = 5 plotted versus controlparameter c , black crosses labeled by H L and H U correspond to Hopf bifurcation points ( ℜ λ = 0, ℑ λ = 0), LPC corresponds to the value of c where thefold bifurcation of two limit cycles occurs. (b) Thestable limit cycle emerges at a supercritical Hopf pointlabeled as H U . The continuation of the stable limitcycle shown by blue solid lines . With the decrease in c the unstable limit cycle occurs from subcritical H L point as indicated by red dashed lines of smallamplitude in zoomed view. Black dotted (continued) for the intermediate value of c = 0 . µ M is shown inFig. 3 (c). After reaching H L point at c = 0 . µ MFIG. 4: line corresponds to the fold (limit point)bifurcation of two limit cycles marked as LPC. (c)
Thereal parts of the first two eigenvalues of the system,which are zeros at Hopf points. (d)
The imaginaryparts of the eigenvalues of the system are non-zeros atHopf points. Two saddle-node bifurcation points arealso present in the diagram at c ≈ . µ M, however,we do not show them in the plot as these points are notrelevant to the analysis below and will no longer bediscussed.the unstable limit cycle occurs as marked by red dashedcycles of small amplitude in zoom of Fig. 4 (b). Thestable cycle emerging from H U collides with the unstableone emerging from H L close to c = 0 . µ M. This foldbifurcation of limit cycles is labeled by LPC (marked by black dotted trajectory ). The limit cycle solution ceasesto exist while decreasing c value lower then LPC point.Thus, we conclude that the stable limit cycle solutionsexist in a small range of [Ca ] concentration betweenH U and LPC points.The bifurcation diagrams obtained from DYK [7] andLi-Rinzel [8] models are different from the one presentedabove, owing to different assumptions about the speedof dynamical variables and different reduction schemesemployed. In [7] and [8] [Ca ] and [IP ] are dynamical,whereas we keep [IP ] at a high constant level. More-over, the earlier models assume a well-mixed membraneand do not account for the clustering of channels thatwas introduced in [23] setting the inactivation rate to befrozen at ˜ k + i = k + i c s . We study the different Ca dy-namics associated with IP R clusters and physiologicalparameters for which the transitions for different kindsof behavior may occur.The phenomena described above shows the Ca oscil-lations in a cytosolic [Ca ] at a single cluster. The am-plitudes of oscillations lie within a physiologically plausi-ble range and their periods are similar to the oscillationsfound in several cells types [38]. However, various stud-ies [12, 13] imply that Ca oscillations are collective innature, and emerge when multiple clusters are involved.In our model, we, therefore, cannot allow individual clus-ters to oscillate autonomously. We shall show how it isonly via the coupling between clusters that oscillationsemerge, and switching off the inter-cluster coupling leadsto the loss of rhythmic behavior. Thus, we next studythe behavior of Ca release from arrays of clusters. B. Cluster-based model of a discrete chain ofclusters
Here we study the dynamics of a chain of clusters cou-pled by diffusing Ca ions. Firstly, we start by analyz-ing a simple example to illustrate the effect of couplingbetween neighboring clusters. Consider L clusters cou- pled with the diffusion term in eq. (5). The discrete formof eq. (5) for the i th cluster is expressed in terms of thelocal Ca concentration c i whose rate of change is pro-portional to the spatial gradient of the diffusion current J i for i = 1 , . . . , L . The continuum diffusion equation isderived from ∂c i ∂t = 1 ℓ ( J i − J i − ) = 1 ℓ h D (cid:18) c i − − c i ℓ (cid:19) − (7) − D (cid:18) c i − c i +1 ℓ (cid:19) i . In our model with a discrete distribution of clusters, weshall take ℓ ≈ . µ m as the optimal distance betweenclusters adapted for the current model and D , the ef-fective diffusion constant for Ca ions is taken to be D = 30 µ m / s − [39]. Here we assume the clusters to bearranged in a ring topology, c i + L = c i . We numericallysolve the system (1), (2) and (5) assuming the diffusionterm as in eq. (8) and model parameters as in the pre-vious section. Only the equilibration rate is modifiedfor these simulations and set to λ = 230 s − , instead of λ = 10 s − that was used earlier.Here we consider the discrete lattice model for L = 2, atypical unit cell of a discrete chain of clusters. The modelparameters for both clusters are the same apart from theactivatable numbers of channels in clusters, which are N = 9 ( N = 5) in the first (second) cluster. In this case,our model is given byd a i d t = f ( a i , h i , c i ) , i = 1 , h i d t = g ( a i , h i , c i ) , i = 1 , c d t = − λ [ c − c d ( a )] − Dℓ ( c − c ) , (10)d c d t = − λ [ c − c d ( a )] − Dℓ ( c − c ) . (11)We present the results of this unit cell L = 2 model inFig. 5 (a), (b) and (c). Note, that the unit cell case is aparticular example of a chain of clusters of length L . Weuse the discrete chain approximation to illustrate the be-havior of the clusters in a “one-dimensional membrane”.Unlike the single-cluster approach, calcium concentra-tion c in the cytosol is explicitly added to the model, theequilibration rate λ is reduced and the diffusion term isintroduced in eqs. (10), (11). The assumption that c willbe “slaved” to c d is inapplicable; we need to account for aseparate dynamics of Ca cytosolic concentration awayfrom the steady-state value c d .The oscillatory trajectories are exhibited in clusters 1[Fig. 5 (a)] and 2 [Fig. 5 (b)] as the consequence of adiffusive interaction between the clusters. Notably, theconcentration of Ca at cluster 1 exhibits small ampli-tude oscillations [Fig. 5 (c) top ] around the fixed pointwith high c [here the fixed point is associated with thehigher fixed point in Fig. 3 (b) for the single cluster case H L H U H L H U H L H U (cid:30)(cid:31) !" (c) c , μ M t, s c , μ M )*+, - . t, s / 379: ;<= c>?@ABC D h a c , μ M N = 5 c , μ M EFGHIJK L h a M N
D = 30 μm /s, λ = 230 s -1 ℜℜℜ
ℑℑℑ
FIG. 5: (Color online) The initiation of Ca oscillations in the array of two clusters D = 30 µ m / s, λ = 230 s − . (a) The trajectory in cluster 1 with N = 9 channels. (b) The oscillatory trajectory in cluster 2 ( N = 5) triggered bythe diffusion from the first cluster. (c) Oscillating Ca profiles in both clusters. (d) Bifurcation diagram of thesystem of two coupled clusters plotted versus diffusion coefficient D . The points H L and H U represent Hopfbifurcation points ( ℜ λ = 0, ℑ λ = 0). (e) The continuation of a limit cycle (blue lines) performed by analogy toFig. 4 (b) varying D . Similarly, the black markers correspond to the fold bifurcation of limit cycles (LPC). (f ) Thereal ( top ) and imaginary ( bottom ) parts of the eigenvalues of the system.but now we solve the 6–variable system]. The influx fromthe cluster 1 causes the emergence of Ca oscillationsin cluster 2 [Fig. 5 (c) bottom ]. In order to study this ef-fect, we perform the continuation analysis of the systemof two coupled clusters given by eqs. (8)–(11) similarly toa single cluster approach described earlier. The dimen-sionality of the system used in the bifurcation analysisremains 6, but we represent the results as two (three)dimensional plots to depict the governing dynamics ateach cluster. We plot the bifurcation diagram shown inFig. 5 (d) which represents the dependence of the a valueof the fixed point versus D as the control parameter. Theeffect of Ca flux from cluster 1 to cluster 2 is pro-portional to D which controls the strength of coupling.We observe a similar diagram as shown in Fig. 4. Wealso observe H U and H L points where stable and unsta-ble limit cycles occur, respectively. From the eigenvaluesplotted in Fig. 5 (f), we confirm that H U and H L areHopf bifurcation points. The continuation of the sta-ble limit cycle occurring at H U [Fig. 5 (e)] shows thatthe periodic solutions in the cluster 2 exist in the range(29 ≤ D ≤ µ m / s. This brings the range of concen-trations of Ca to the one shown in Fig. 4 for the same λ . Thus, we conclude that the emergence of oscillationsin the second cluster depends on the concentration influxfrom cluster 1 to cluster 2 which raises the level of con-centration in cluster 2 to a range that drives oscillationsin the single cluster approach.Having accounted for the ability of diffusive couplingsto bring about oscillations in a coupled cluster systemwhen the individual clusters exhibit non-oscillatory dy-namics on their own, we extend the model to a continuumdescription. In the next subsection, we apply this mech-anism to a domain of the membrane in order to quali-tatively demonstrate the emergence of Ca waves andoscillations in the ER membrane. C. Calcium Waves and Oscillations in thecluster-based model
As shown above, the behavior of the three-state modelfor a single cluster depends on the number of channels N in a cluster that are IP bound and may be activated ordeactivated by Ca , namely, excitable for small N tran-sitioning for large N to bistable dynamics. While the0ER membrane would typically contain clusters of differ-ent numbers of channels, we study the simplest heteroge-neous scenario with one bistable cluster with N = 9 andthe multiple excitable clusters with N = 5. Here andfurthermore we refer to a N = 9 cluster as bistable anda N = 5 cluster as excitable as per their behavior in theisolated single cluster case, even though it is the entiresystem whose stability matters. The choice of this figuresis based on the study [24], where the average number ofactivatable channels in Xenopus oocytes ER membranestend to N = 5 with the increase in [IP ], the distribu-tion of numbers of channels appears to be very close toa Poisson distribution [3, 24]. We chose this numbers tobe consistent with [24] and propose the modeling of morecomplex heterogeneous systems for future study. In thechosen situation the excitable behavior may be associ-ated with short-lasting puffs while the bistable behaviormay be associated with waves or long-lasting puffs, wherethe exit from the long-lasting puffs that is thought to bedriven by dissociation of IP or other mechanisms whichare discussed in subsection III D. The exit from the longelevated high Ca concentrations might be performedthrough the other mechanisms such as change in couplingstrength α , maximal Ca elevation c s or transition rate k + i .In this section, we model a domain of an ER membranewhich is assumed to contain multiple excitable clusters( N = 5) and the only one bistable ( N = 9) cluster in or-der to study the interplay between excitable and bistableclusters. We demonstrate the emergence of oscillationsand quantify the characteristics of the phenomena suchas front velocity and period.We numerically solve the eqs. (1), (2) and (5) inthe chosen one-dimensional domain with open bound-aries and Neumann boundary conditions c ′ x ( x = 0 , L ) =0. Let us assume the part of ER membrane with theinitial distribution of cytosolic [Ca ] c ( x,
0) = c + A exp {− [ x / (2 σ )] } , where A = 2 µ M, σ = 0 . µ m, whichtravels from the left part of a membrane as shown inFig. 6 (a). Calcium release upon activation occurs onlyat channel clusters i in eq. (5) via c d ( a ) and cytosolicconcentration of Ca is raised above basal level c bydiffusion between clusters a constant distance ℓ = 3 . µ mapart. The results depicted in Fig. 6 are obtained forequilibration rate λ = 230 s − . The effect of changingthe values of D , λ , and ℓ upon qualitatively differingCa release response is discussed further in this section.In Figs. 6 (c) and (d) we observe the emergence ofCa oscillations at the bistable cluster ( black dot ) trig-gered by the initial front [Fig. 6 (a)]. The front occursfrom an initial Ca increase close to the point x = 0 µ mand diffuses throughout the membrane. While the ini-tiation event is set by choosing initial condition in ourdeterministic model it reflects the potential occurrenceof an initial peak by either an internal or external stimu-lus or a stochastic fluctuation. After passing the bistablecluster the initial front raises the concentration at thiscluster to the value of c which corresponds to the up- secondary front secondary frontinitial front OPQRST initial front[Ca ]N=5N=9 UVW secondary frontsecondary front
XYZ x, x, x, t = 3.4 st = 4.4 st = 6.4 st = 9.4 s
FIG. 6: (Color online) Ca front propagation fromthe spatial model for two types of clusters( D = 30 µ m / s, λ = 230 s − ). The orange dots represent excitable clusters with N = 5 channels as longas the large black dot corresponds to a bistable clusterwith N = 9 channels. (a) The propagation of a frontcaused by the raised initial [Ca ] in the left part of themembrane. (b) The Ca peak occurs at the bistablecluster. (c) The diffusion from the peak causesinitiation of a secondary wave. (d)
The frontpropagates after that subsequent front occurs. Thewaves occur periodically from the bistable cluster.per fixed point in Fig. 5 (a). This behavior induces along-lasting state which is terminated by IP unbinding,a feature which we discuss later. Therefore, in Fig. 6 (b)a residual peak occurs at cluster N = 9. Thereafter, twosecondary fronts emerge from the residual peak as shownin Fig. 6 (c). As we consider deterministic oscillationsthis process repeats periodically with a constant period[Fig. 6 (d)].This behavior can be observed because of two differentregimes underlying Ca release in clusters: excitable forshort-lasting puffs and bistable for long-lasting events.Our model shows that bistable clusters can produce long-lasting puffs that might cause the emergence of Ca oscillations in cellular membranes.The diffusion coefficient D and relaxation parameter1 λ incorporate spatial effects without explicitly modelingbuffers in the membrane. As these are likely to influencethe characteristics of the Ca waves and oscillations wevary D to study how these characteristics change. Fur-thermore, by changing equilibration rate λ we can modelthe inclusion of a linearized reaction Ca leak flux andpump action terms in the eq. (5). There are differenttimescales for the activation and inhibition of the chan-nels “hidden” in eqs. (1), (2). The action of pumps occurson much slower timescales compared to λ chosen.TABLE II: Parameters used in the cluster-basedreaction-diffusion model. Parameter Value Description c . µ M Rest level of [Ca ] N bound) λ
230 s − Equilibration rate D µ m / s Diffusion coefficient ℓ . µ m Inter-cluster distance r cl . µ m Cluster radius In Fig. 7 (a) we show the values of the average time cal-culated between two subsequent Ca release events ateach cluster. This value appears to be almost constant.We map the change in a period of Ca oscillations tophysiological ranges of D and λ [Fig. 7 (a)]. Dark regionscorrespond to a non-oscillatory response, while oscilla-tions occur within the purple triangular zone. A curiouseffect is observed at the lower edge of the triangle in therange (100 ≤ λ ≤ − . The light spots marked bywhite circles correspond to the region on the diagramwhere oscillations exhibit periods ≈
10 s and belong toCa alternans – a pulse travels alternately to the right,followed by another to the left and so on (Fig. S2). Un-der given initial conditions, this effect can be explainedby the fact that the concentration diffusing from the ir-regular cluster is sufficient to excite a neighboring clusteronly from one side at a time. Thus, we observe the oscil-lations occurring on alternate sides of the bistable clustersequentially. In Fig. 7 (b) we show the velocity of prop-agation of the initial front for a similar range of D and λ : the dark regions corresponding to abortive waves orpuffs and the light trapezium corresponds to propagatingCa waves.The distance between clusters in the various cell typeslies in range 1 − µ m [40]. Figs. 7 (c) and (d) show thedependence of the periods and the velocities on λ for dif-ferent values of the distance ℓ between the clusters. Inboth graphs, we consider the effective diffusion constantto be fixed at a typical value of D = 30 µ m / s. The pe-riods of the oscillations [Fig. 7 (c)] generally reach largervalues for long distances between clusters. This can beexplained by the fact that the front travels a longer dis-tance in order to reach neighboring clusters. Interest- [\] l^_ ‘ λ , s -1 ) b e l o c i t y , μ m / s (b)(c)(a) D, μm /s d , s - Velocity, μm/s D, μm /s
10 5020 30 4001038006004002000 45 036 ef gij kmn , s - D, μm /s o , s - pqrvwxyz{|}~(cid:127)(cid:128) (cid:129)(cid:130)(cid:131)(cid:132) (cid:133)(cid:134) (cid:135)(cid:136)(cid:137)(cid:138)(cid:139)(cid:140) (cid:141) (cid:142)(cid:143)(cid:144) (cid:145) (cid:146)(cid:147)(cid:148) (cid:149) (cid:150)(cid:151)(cid:152) (cid:153)(cid:154)(cid:155)(cid:156)(cid:157)(cid:158)(cid:159)(cid:160) ¡¢£⁄ ¥ƒ§¤'“«‹› log ( fi , s -1 ) fl(cid:176)–†‡·(cid:181)¶ • ‚ μ „” » μ …‰ (cid:190) μ ¿(cid:192) ` μ ´ˆ ˜ μ ¯ FIG. 7: (Color online) (a)
The phase diagram whichrepresents different kinds of Ca dynamics such aspuffs, waves, and oscillations. The domain with zeroperiods [ dark-blue zone in (a)] contains both puffs andwaves in a co-existing manner. The triangular light-bluedomain with non-zero periods contains oscillatorybehavior of the system. The white circles indicate theregions with Ca alternans. (b) The phase diagramrepresenting two regimes of Ca dynamics such aspuffs and waves in the case when all the clusters areexcitable. Zero velocity domain in dark-blue corresponds to puffs and abortive waves, while thenon-zero-velocity domain (shown in blue to red colors)represents waves. (c) The dependence of the period on λ for different inter-cluster distances, D = 30 µ m / s. (d) The dependence of velocity on λ , the distances, andthe diffusion coefficient are the same as in (c).ingly, the oscillatory behavior exists in wider ranges of λ for smaller distances between the clusters [Fig. 7 (c)]; os-cillations are more robust to the equilibration process ifthe distances between clusters are smaller. The same pro-cesses define the velocities [Fig. 7 (d)]; these are higherfor shorter distances because a front reaches a neighbor-ing cluster faster and is able to trigger a response there.2The data with comparatively large periods ≈ − ℓ = 5 µ m, ℓ = 4 µ m and ℓ = 3 µ m in Fig. 7 (c) exhibits thenon-regular alternans emerging similarly as in Fig. 7 (a). D. IP dependence The assumption of high levels of [IP ] leads to a sus-tained Ca release from the bistable cluster. In the cur-rent subsection, we aim to link the termination of thisevent to the change in the [IP ] levels [ red dashed lines in Figs. 8 (d), (e)] on timescales slower than the oscil-lations found in previous subsections [ blue solid lines inFigs. 8 (d), (e)].The extension of the model is sketched in the Fig. 8 (a).In the case of two coupled clusters the IP dependentmodel (8), (9) reads asd a i d t = f p ( a i , h i , y i , c i ) , i = 1 , h i d t = g p ( a i , h i , y i , c i ) , i = 1 , y i d t = C + ( c i )(1 − y i ) − C − ch ( c i , p ) y i , (14)where f p and g p are given in eqs. (A34) and (A35), C + and C − ch are given in eqs. (A23) and (A29), respectively.The IP changes as p ( t ) = p , if t ≤ t dur , p a σ σ +2 D p t e − x − y Dpt +2 σ , otherwise , (15)where the initial value of [IP ] at a cluster p = 3 . µ M,the duration of the IP stimulus is t dur = 10 s. We as-sume the shape of the initial [IP ] distribution to beGaussian of amplitude p a = 4 . µ M and with variance σ = 1 µ m (see Supporting Material Sec. S3). Afterstimulus terminates, [IP ] diffuses in infinite 2D domain( D p = 10 µ m / s [41]). We calculate time traces shownby red dashed lines in Fig. 8 (d) and (e) at clusters posi-tioned at x = ± ( ℓ + r cl ) / y = 0 (here the distancebetween clusters ℓ = 1 . µ m, cluster radius r cl = 0 . µ m).Ca dynamics remains unaltered as in eqs. (10), (11).The trajectories in ( a, h, y ) phase space for both clus-ters are exhibited in Figs. 8 (c), (d). The global oscil-lations occur in our model at [IP ] ≈ − µ M, which iscompatible with the levels considered in [42]. The termi-nation of the oscillations [solid blue lines in Fig. 8 (e), (d)]at both clusters is caused by the decrease in [IP ] undera level of several tens of nM, which is typically takenas the rest value of the [IP ] [24]. It might be usefulto reconsider the influence of the calcium pumps in thetermination process as we had argued that their effectsare subsumed in a linearized equilibration rate λ . Wehave checked that the effect of introducing an explicitpump term [7] on the wave termination behavior can beaccommodated by adjusting λ . IV. DISCUSSION AND CONCLUSIONS
In this paper, we presented a reaction-diffusion modelof a membrane containing clusters of IP receptor chan-nels of the fire-diffuse-fire [21] type but with a mech-anistic description of channel firing similar to the spa-tial Fitzhugh-Nagumo model [43, 44]. This mechanisticmodel was built on top of a dynamical systems analysisof a previously proposed simplification [23] of the DYKmodel [7]. We noticed that incorporating the transitionrates given in [24] into the model revealed the emergenceof Ca oscillations even in a single cluster which inspiredthe membrane level formulation of a model of Ca puffs,waves, and oscillations whose characteristics were quanti-fied under varying cluster configurations and physiolog-ical conditions. Our model allows us to observe Ca waves with the velocities in agreement with experimentsby Marchant et al. [42] for various distances between clus-ters and equilibration rates [Figs. 7 (b) and (d)].Here we propose the mechanism of initiation of oscil-latory behavior in the membrane that consists of the fol-lowing sequence. • Upon receiving an IP stimulus the membrane con-figures IP receptor channels to form clusters withdifferent numbers of channels [3] with bound IP that are primed for activation by Ca binding. • If the number of primed channels in a cluster ishigher than a threshold value, then the cluster isbistable. When an initial surge of Ca passes thebistable cluster [Fig. 6 (b)] it drives Ca concen-tration there into the higher fixed point of the dy-namics as shown in Fig. 5 (a). Sustained Ca el-evation at a bistable cluster corresponds to a long-lasting release event. • The residual [Ca ] at the bistable cluster diffusesto the neighboring clusters. • This drives [Ca ] at the neighboring clusters intoan oscillatory regime as in Fig. 5. These oscillationsspread to further clusters due to diffusive coupling.It is generally known that high Ca concentrationsmight lead to cell death and this renders the presence ofbistability and consequent sustained Ca release in ourmodel an unhappy feature. Due to fast luminal Ca re-filling [45] sustained Ca release, such as what followsfrom a transition to the upper fixed point of the bistableregime, cannot be terminated by ER depletion at thetimescales considered by current paper. The mechanisticorigins of our model enable us to propose possible phys-iological possibilities that might explain the terminationof Ca release from the bistable cluster. • IP unbinding. In the full DYK approach, R¨udigeret al. [32] attribute the termination of long-lastingCa release events to IP unbinding. A simplifiedmodel with inclusion of the lower plane of the DYK3 ˘˙¨ (cid:201) a ˚ c ¸ (cid:204) ˝˛ˇ c — a - (cid:209)(cid:210) (cid:211) ~ a hz C - c h ( c , p ) C + ( c ) (e) c l u s t er D = 30 μm /s, λ = 230 s -1 (d) c l u s t er t, s t, s y N=9 [Ca ], μ M[IP ], μ M N=5(b) h a cluster 1 N=9 h a y (c) cluster 2 N=5 FIG. 8: (Color online) The IP dependent model for two coupled clusters. (a) The scheme of an extension of thethree-state model introduced in Fig. 1 (c). The transitions between the upper and the lower plane of the DYK cubeare governed by the steady-state rates C + and C − ch in analogy to Li-Rinzel [8]. (b) , (c) The trajectories in the( a, h, y ) phase space, where y corresponds to the fraction of the channels in a cluster occupying the lower plane ofthe DYK cube, N is the number of activatable channels in a cluster. (d) , (e) The oscillatory regime correspondingto Fig. 5 ( blue solid lines ) for high [IP ] is terminated due to the decrease in the [IP ] ( red dashed lines ).cube shown in subsection III D has also displayedthe termination of a Ca wave in the bistable clus-ter in Fig. 8. Thus, we argue that neglecting IP unbinding completely might lead to results whichseem unrealistic. Therefore, the limits of the ap-plicability of our model are bound to the change ofthe IP levels in the cell. • The uncoupling of IP R channels.
Our bifurca-tion analysis shows that the system ceases to bebistable if the coupling between channels α de-creases. Hence, if α decreases on timescales slowerthan oscillation periods observed in the model, thelong-lasting event can terminate (see Fig. 9). Inthis paper, we have considered α to be constant,which might not always be the case in real clus-ters. It is known from experiments that the cou-pling (uncoupling) of IP R channels is caused byIP binding (unbinding) [3, 46, 47]; the mechanismswhich connect channels remain largely unknown,however. If an increase in p leads to a fast (hun-dreds of ms according to [47]) increase in couplingstrength α could acquire a value as taken in thispaper. Following the initial stimulus that triggersthis increase of IP and its consequent effect on α , agradual diminution of IP levels could still restrictthe dynamical states to the upper place of the DYK cube, while reducing the coupling between channelsin a cluster. This would lead to the termination ofpuffs on the timescales defined by the change in[IP ] levels. However, only decreasing α under 0 . concen-tration at a single cluster, leads to the exit fromthe bistable state. A mechanism which causes de-coupling of channels in large clusters has been dis-covered in the RyR clusters [48] although a similarproperty has not been reported for IP receptors. • ER depletion.
In this study, we assume the sameavailability of Ca in the luminal pool in (4).Even though the ER pool can be large and itsfast refilling can prevent local depletion [45], otherstudies suggest that the large pool can be com-pletely depleted by high-amplitude sustained Ca releases [49]. • Stochastic termination.
The bistability found inthe single-cluster model can also be terminatedstochastically. Fig. 10 shows the distributionof inter-puff intervals (IPIs) obtained from theLangevin model based on eqs. (1)–(2). We con-sider the additive Wiener noise as shown in (B1),(B2) and compare the results with experimentalpuff data [5]. The duration of the events lies in4 (a)(b)(c)
FIG. 9: (a)
The change of the coupling strength α ineq. (4) on a slow timescale caused by the uncoupling ofIP R channels due to IP unbinding. (b) Thetermination of the sustained Ca release at thebistable cluster shown in Fig. 6. (c) The termination ofthe oscillatory Ca release in one chosen excitablecluster shown in the Fig. 6.a range of hundreds of ms which is also consistentwith the experimental data available.There are other consequences for the fixed levels ofIP in the model. Due to the fixed high level of [IP ]assumed, our model exhibits periods of oscillations lowerthan the wide range found in typical experimental stud-ies. Although the high IP assumption narrows the rangeof periods observed, it helps us to study the main char-acteristics of waves and oscillations phenomenologically.The variability of the periods can be achieved by incor-porating more complex features of Ca signals such asstochasticity [12, 13] and array enhanced coherence res-onance [11] under dynamical [IP ] loads.There are other characteristics of the emergent Ca waves and oscillations in our model that matches thosereported in the extensive experimental literature [38, 42,50, 51]. In particular, Marchant and Parker [51] haveshown that an increase in [IP ] shortens intervals betweenCa release events. Our model agrees with these exper-iments by exhibiting oscillations with frequency increas-ing with increase in [IP ]. The periods of oscillationsare also consistent with R¨uckl. et al. [24]. The ampli-tudes of oscillations are in a range of observed puffs andwaves [52]. Also, the effect of wave collision found exper-imentally is present in our model.In this paper, a reduced description of clustered chan-nels with a parametric representation of the dynamics ofsubunits of channels is present alongside membrane-leveldiffusion-enabled interactions. While the simplifications simulat ions N=7experiment Dickinson et al. [5] N o . o f ob s e r v a ti on s IPI, s120
FIG. 10: (Color online) Inter-puff intervals of thestochastic single-cluster model (B1), (B2), (3)containing N = 7 channels with additive Wiener noise. Grey bars are fitted with the blue solid curve andrepresent the distribution of IPIs obtained from themodeling of 622 puffs.
White bars with black solid fitare reproduced from the experiments reported in [5].introduced in the reduction scheme can lead to model be-havior inconsistent with empirical observations, the dy-namical systems analysis presented here as a means oflinking phenomena at different scales also enables us topropose explanatory mechanisms that could be accom-modated in more elaborate models.See Supplemental Material for the extended bifurca-tion analysis of the single-cluster model (Subsec. III A),the figure illustrating alternans (Subsec. III C), thechange of IP (Subsec. III D) derived from diffusion equa-tion. ACKNOWLEDGMENTS
S. B. would like to acknowledge Alex Vasilev for inspi-ration for this research. Also, Andrii Iakovliev for usefulideas and proofreading the manuscript. Current researchhas been funded by Erasmus Mundus ACTIVE grant andFaculty of Engineering and the Environment at the Uni-versity of Southampton.
Appendix A: DYK reduction to the three-statemodel
The kinetics governing the upper plane of DYK cube[Fig. 1 (b)] is formulated as the next set of differentialequations5˙ x = − x ( b + a c + b ) + x a p ++ x b + x a c, (A1)˙ x = − x ( b + a c + a c ) + x a p ++ x b + x b , (A2)˙ x = − x ( b + a c + a c ) + x a c ++ x a p + x b , (A3)˙ x = − x ( b + b + b ) + x a c ++ x a p + x a c, (A4)where x + x + x + x = 1 − y , y = x + x + x + x . The fractions of subunits in state ijk aredenoted as x ijk .In the case of high [IP ] the IP binding sites are sat-urated. Therefore, we can apply the condition of de-tailed balance between upper and lower planes of DYK: a i px jk = b i x jk for four sets of i , j , and k such as i = 1, j = 0 , k = 0 and j = 3, j = 0 , k = 1. From theseconditions we derive x = x b a p , (A5) x = x b a p , (A6) x = x b a p , (A7) x = x b a p . (A8)We substitute eqs. (A5)–(A8) into the system ofeqs. (A1)–(A4). The resulting four state model appearsas ˙ x = − x ( a c + b ) + x b + x a c, (A9)˙ x = − x ( a c + a c ) + x b + x b , (A10)˙ x = − x ( b + a c ) + x a c + x b , (A11)˙ x = − x ( b + b ) + x a c + x a c, (A12)where x + x + x + x = 1 − y due to conservationof probability, y = x + x + x + x .Now we introduce a reduction of DYK similarly to Liand Rinzel [8]. However, here we separate the upper andlower planes of the DYK cube [Fig. 1 (a)]. Consideringthe lower plane of the DYK cube in a steady state weobtain x = d d y ( c + d )( c + d ) , (A13) x = d cy ( c + d )( c + d ) , (A14) x = d cy ( c + d )( c + d ) , (A15) x = c y ( c + d )( c + d ) , (A16) where d i = b i /a i , i = 1 , x = d d (1 − y )( c + d )( c + d ) , (A17) x = d c (1 − y )( c + d )( c + d ) , (A18) x = d c (1 − y )( c + d )( c + d ) , (A19) x = c (1 − y )( c + d )( c + d ) . (A20)To calculate the transition rates between these twogroups we use steady state equations for upper plane(A14)–(A16) and lower plane (A18)–(A20) K − = a p ( x + x ) + a p ( x + x ) == ( a c + a d ) pyc + d , (A21) K + = b ( x + x ) + b ( x + x ) == ( b c + b d )(1 − y ) c + d , (A22)here we define transition rates between upper andlower planes of the DYK cube as C + = b c + b d c + d , C − = ( a c + a d ) pc + d , (A23)which is used in (14), the corresponding transition ratesare given in Tab. III.The results presented by the system of eqs. (A9)–(A12)are applicable to a subunit of a channel. In order to de-scribe the behavior of a cluster of channels, similarly to[23], we introduce variables a , g , h ′ and z that representTABLE III: The parameters of the IP dependentmodel adapted from [24, 25]. Parameter ValueFirst order rates, µ M − × s − a a a − b × − b b . × − Dissociation constants d i = b i /a i , µ M d . d d . dependence the sys-tem (A10)–(A12) transforms tod a d t = k + a cz − k − a a + k − i g − ˜ k + i a, (A24)d g d t = k + a ch ′ − k − a g − k − i g + ˜ k + i a, (A25)d h ′ d t = k − a g − k + a ch ′ + k + i cz − k − i h ′ , (A26)d y d t = C + (1 − y ) − C − ch y, (A27)where z = 1 − a − g − h ′ − y as sum of the fractionsof all states to unity, y defines the fraction of IP R un-bound channels (occupying the lower plane of the DYKcube), C + and C − ch are given by eqs. (A23) and (A29),respectively.The kinetics of these variables is governed by rates k ± a , k ± i and concentrations c and c s , where we substitute twoparameters k + i and c s as ˜ k + i = k + i c s . These transitionrates are obtained from the original DYK rates. Earliersimulations [23] showed that most of the reactions ( k − a , k + i , and k − i ) in the four-state model are obtained from b , a , and b , respectively. In contrast, k + a rate has beenconstructed from the condition of channel opening (atleast three subunits should be active). The probabilityof the single subunit activated is P act = a c/ ( a c + b ).According to the binomial distributions, the probabilityof none of the subunits is active is P = (1 − P act ) . If onlyone of four subunits is active and all remaining ones areinactive P = 4 P act (1 − P act ) . There are C = =6 ways for two of four subunits to be activated, thus, P = 6 P (1 − P act ) . Therefore, the transition to thisstate is possible only after the activation of two subunits.The probability that two subunits are active under thecondition that the channel is closed is P (2 |{ , , } ) = P / ( P + P + P ).We calculate k + a as a transition rate from a closed chan-nel with two active subunits to an open channel withthree activated subunits. One of two remaining inactivesubunits might be activated, thus, we obtain k + a = 2 a P (2 |{ , , } ) . (A28)Similarly to k + a we derive C − ch which requires four sub-units to bind IP [26] C − ch = C − P IP (3 |{ , , , } ) , (A29)where C − is given by eq. (A23), P IP (3 |{ , , , } ) = P IP3 / ( P IP0 + P IP1 + P IP2 + P IP3 ), similarly P IPact = C − / ( C + + C − ), P IP0 = (1 − P IPact ) , P IP1 = 4 P IPact (1 − P IPact ) , P IP2 =6( P IPact ) (1 − P IPact ) , P IP3 = 4( P IPact ) (1 − P IPact ).In order to reduce the number of variables, similarlyto [23] we introduce a compound state h which contains g and h ′ [Fig. 1 (d)]. The effective rates of transitions tothis state are obtained from a detailed balance between g and h ′ . We introduce the fraction of channels in thecompound state as the sum of components h = g + h ′ ,where from detailed balance h = k − a k + a c g + g = k + a c + k − a k + a c g, (A30) h = h ′ + k − a k + a c h ′ = k + a c + k − a k − a h ′ . (A31)Considering that the rates of transition to compoundstate are k = k − i k + a ck + a c + k − a = k − i g , (A32) k = k − i k − a k + a c + k − a = k − i (1 − g ) , (A33)where g = k + a c/ ( k + a c + k − a ).The resulting model reads asd a d t = k + a ( c )(1 − a − h − y ) − k − a a + k ( c ) h − ˜ k + i a = △ f p ( a, h, y, c ) , (A34)d h d t = k + i c (1 − a − h − y ) − k ( c ) h − k ( c ) h + ˜ k + i a = △ g p ( a, h, y, c ) , (A35)which where f p and g p are used in eqs. (12), (13).The solution of eq. (A27) separately is y ( t ) = Ae − ( C + + C − ) t + C + C + + C − . (A36)It follows that the equilibrium solution for the systemis y eq = C + C + + C − = ( b c + b d )( c + d )( a d + a c )( c + d ) p + ( b d + b c )( c + d ) p ≫ −→ constp , (A37)in the case of IP saturation, we consider p = [IP ] to be high. Taking that into account and assuming that7waiting time is long enough t ≫ / ( C + + C − ) we obtain y →
0. Considering all the assumptions made before,we further operate with the three-variable system takinginto account the summing fractions to unity a + h + z = 1.Using the assumption of high [IP ] from (A34), (A35),(A27) we obtain the model illustrated in Fig. 1 (d) andgiven by eqs. (1), (2). Appendix B: Stochastic model
In this appendix, we introduce the Langevin modelof Ca release from a single cluster of channels. Byintroducing the noise terms into the model from eqs. (1),(2) we obtain d a = f ( a, h, c ) d t + As a d W, (B1)d h = g ( a, h, c ) d t + As h d W, (B2) where d W are Wiener increments (generated by theWiener process also known as Brownian motion [53]), s a and s h are the functions of the system state in thegeneral case of multiplicative noise, but here for simplic-ity, we consider s a = s h = 1 which corresponds to theadditive noise case, A = 0 .
37 is the amplitude of thenoise.We integrate the system (B1), (B2) under c = c d as-sumption using the Euler-Maruyama method [54] withparameters given in the last column of Tab. I and N = 7.The results presented in Fig. 10 are obtained by analyz-ing the intervals between 622 puffs using [Ca ] traces ofseveral thousands of seconds. [1] M. J. Berridge, Biochem Soc T , 297 (2012).[2] K. Mikoshiba, Biochem Soc Symp , 9 (2007).[3] Taufiq-Ur-Rahman, A. Skupin, M. Falcke, and C. W.Taylor, Nature , 655 (2009).[4] I. F. Smith, D. Swaminathan, G. D. Dickinson, andI. Parker, Biophys J , 834 (2014).[5] G. D. Dickinson, D. Swaminathan, and I. Parker,Biophys J , 1826 (2012).[6] M. J. Berridge, Biochem Soc Symp , 1 (2007).[7] G. W. De Young and J. Keizer,P Natl A Sci USA , 9895 (1992).[8] Y.-X. Li and J. Rinzel, J Theor Biol , 461 (1994).[9] A. Politi, L. D. Gaspers, A. P. Thomas, and T. H¨ofer,Biophys J , 3120 (2006).[10] E. Smedler and P. Uhl´en,BBA-Gen Subjects , 964 (2014).[11] M. Falcke, Biophys J , 42 (2003).[12] A. Skupin, H. Kettenmann, U. Winkler, M. Wartenberg,H. Sauer, S. C. Tovey, C. W. Taylor, and M. Falcke,Biophys J , 2404 (2008).[13] K. Thurley, I. F. Smith, S. C. Tovey, C. W. Taylor,I. Parker, and M. Falcke, Biophys J , 2638 (2011).[14] R. Kupferman, P. P. Mitra, P. C. Hohenberg, and S. S.-H. Wang, Biophys J , 2430 (1997).[15] A. Atri, J. Amundson, D. Clapham, and J. Sneyd,Biophys J , 1727 (1993).[16] J. Sneyd, P. D. Dale, and A. Duffy,SIAM J Appl Math , 1178 (1998).[17] G. Dupont and A. Goldbeter, Biophys J , 2191 (1994).[18] A. Calabrese, D. Fraiman, D. Zysman, and S. PonceDawson, Phys Rev E , 031910 (2010).[19] J. Keizer and G. D. Smith, Biophys Chem , 87 (1998).[20] J. Sneyd and J. Sherratt,SIAM J Appl Math , 73 (1997).[21] S. P. Dawson, J. Keizer, and J. E. Pearson,P Natl A Sci USA , 6060 (1999).[22] S. Coombes and Y. Timofeeva,Phys Rev E , 021915 (2003). [23] S. R¨udiger, Phys Rev E , 062717 (2014).[24] M. R¨uckl, I. Parker, J. S. Marchant, C. Na-gaiah, F. W. Johenning, and S. R¨udiger,PLOS Comput Biol , e1003965 (2015).[25] M. R¨uckl and S. R¨udiger, Eur Phys J E , 1 (2016).[26] C. W. Taylor and V. Konieczny, Sci Signal , 2 (2016).[27] Z. Qu, M. B. Liu, and M. Nivala,Sci Rep , 35625 (2016).[28] J. R. Groff and G. D. Smith, Biophys J , 135 (2008).[29] R. Hinch, Biophys J , 1293 (2004).[30] S. M. Wiltgen, G. D. Dickinson, D. Swaminathan, andI. Parker, Cell Calcium , 157 (2014).[31] D. R. Laver, C. H. T. Kong, M. S. Imtiaz, and M. B.Cannell, J Mol Cell Cardiol , 98 (2013).[32] S. R¨udiger, P. Jung, and J.-W. Shuai,PLOS Comput Biol , e1002485 (2012).[33] I. Bezprozvanny, J. Watras, and B. E. Ehrlich,Nature , 751 (1991).[34] S. R¨udiger, J. W. Shuai, and I. M. Sokolov,Phys Rev Lett , 048103 (2010).[35] J. Shuai, J. E. Pearson, and I. Parker,Biophys J , 3738 (2008).[36] A. Dhooge, W. Govaerts, Y. A. Kuznetsov, H. G. Meijer,and Sautois. B., Math Comp Model Dyn , 147 (2008).[37] Y. A. Kuznetsov, Elements of Applied Bifurcation Theory ,edited by J. E. Madsen and L. Sirovich (Springer, 2000)pp. 91–112.[38] N. Callamaras, J. S. Marchant, X.-P. Sun, and I. Parker,J Physiol , 81 (1998).[39] N. Allbritton, T. Meyer, and L. Stryer,Science , 1812 (1992).[40] S. R¨udiger, Phys Rep , 39 (2014).[41] G. D. Dickinson, K. L. Ellefsen, S. P. Dawson, J. E. Pear-son, and I. Parker, Sci Signal , ra108 (2016).[42] J. Marchant, N. Callamaras, and I. Parker,EMBO J , 5285 (1999).[43] R. FitzHugh, Biophys J , 445 (1961).[44] C. K. R. T. Jones, T Am Math Soc , 431 (1984). [45] L. F. Lopez and S. P. Dawson,Phys Biol , 036006 (2016).[46] Y. Tojyo, T. Morita, A. Nezu, and A. Tanimura,J Pharmacol Sci , 138 (2008).[47] I. F. Smith, S. M. Wiltgen, J. Shuai, and I. Parker,Sci Signal , ra77 (2009).[48] S. O. Marx, J. Gaburjakova, M. Gaburjakova,C. Henrikson, K. Ondrias, and A. R. Marks,Circu Res , 1151 (2001).[49] G. Ullah, I. Parker, D.-O. D. Mak, and J. E. Pearson, Cell Calcium , 152 (2012).[50] M. D. Bootman, M. J. Berridge, and P. Lipp,Cell , 367 (1997).[51] J. S. Marchant and I. Parker, EMBO J , 65 (2001).[52] I. Parker and Y. Yao, J Physiol , 663 (1996).[53] B. Øksendal, Stochastic Differential Equations , 6th ed.,Universitext (Springer Berlin Heidelberg, Berlin, Heidel-berg, 2003) pp. 109–112.[54] A. D. Horchler, “Matlab Toolbox for the Numerical Solution of Stochastic Differential Equations (Version 1.2) [Source code].”(2013). r X i v : . [ q - b i o . S C ] S e p Supporting Material for “A phenomenological cluster-basedmodel of Ca waves and oscillations for Inositol1,4,5-trisphosphate receptor (IP R) channels”
Svitlana Braichenko ∗ and Atul Bhaskar Faculty of Engineering and the Environment,University of Southampton, SO17 1BJ, Southampton, UK
Srinandan Dasmahapatra
School of Electronics and Computer Science,University of Southampton, SO17 1BJ, Southampton, UK (Dated: September 14, 2018) ∗ Corresponding author; [email protected]
1. BIFURCATION ANALYSIS
In the Subsec. III A in the main text, we obtain qualitatively different behavior of acluster of channels compared to R¨udiger [1] except for transition from the excitable tothe bistable state of a cluster with an increase in N . Here we propose the more detailedexplanation of such a transition by performing a bifurcation analysis of the two-variablesystem of Eqs. (1)–(2) assuming that c = c d according to Eq. (4). Thus, the system readsas d a d t = f ( a, h, c d ) , (S1)d h d t = g ( a, h, c d ) , (S2) c d = c + 12 αN a ( (cid:20) N a − ǫ (cid:21)) , (S3)The continuation analysis of the system has been performed in order to study the stabilityof a cluster depending on the number of channels it comprises. We plot the bifurcationdiagram of the system in Fig. S1. Three distinctive stability branches are shown in thefigure. For the clusters with N = 5 and fewer channels only one stable branch exists. Thepoint marked by the cross for N = 5 channels corresponds to the stable fixed point inFig. 3(a) in the main text. According to the lower stable ( blue solid ) branch, the positionof the lower fixed point does not depend on N and the eigenvalues of Jacobian are alwaysnegative. With an increase in N saddle-node bifurcation (or limit point marked by LP)occurs. Thus, for N = 6 we observe three fixed points as shown in Fig. 3 (b). Notably, theeigenvalues are negative on stable ( yellow solid ) and have opposite signs on the unstablesaddle ( red dashed ) branch. Also, subcritical Hopf bifurcation in unstable limit cycle occurs.Therefore, we conclude that the transition from the excitable to the bistable behavior ina cluster occurs as the result of saddle-node bifurcation (limit point LP) between N = 5and N = 6 channels. S2. ALTERNATE OSCILLATIONS
Here we present the extension of the results in Subsec. III C. The light spots in the mapshown in the Fig. 7 (a) and higher values of the periods in the Fig. 7 (c) correspond to the2 (cid:212)(cid:213)(cid:214)(cid:215)(cid:216)
Lower stable branchUpper stable branchUnstable branch
LP H λ =-26, λ =-2λ =-1.4, λ =518 λ =-8, λ =-143 FIG. S1. (Color online) The bifurcation diagram of the three-state system depending on thenumber of channels in the cluster N . The lower stable branch ( blue solid line) corresponds to thelower stable fixed point in Fig. 3. The eigenvalues of the Jacobian on this branch are always realand negative ( λ = − λ = − yellow solid line) and unstable ( red dashed line) branches occur as the consequence of the saddle-node (limit point) bifurcation marked by LP.Also, subcritical Hopf bifurcation occurs on a stable branch (marked by H). The positions of fixedpoints for N = 6 correspond to the fixed points in Fig. 3(b). The stability of the branches aredefined by the eigenvalues (for N = 6 eigenvalues λ = − λ = −
143 correspond to stable fixedpoint and λ = − . λ = 518 correspond to unstable saddle point). irregular behavior in Ca oscillations. After initial excitation passes the membrane, theoscillations emerge alternatively from the raised concentration at the bistable cluster firstlytowards x = 0 µ m, then to the other side towards x = 70 µ m. Such behavior causes theperiods of oscillations almost twice higher than in regular oscillations.3 IG. S2. (Color online) The alternate oscillations emerging from the bistable cluster at D =38 . µ m / s, λ = 170 s − . S3. IP3 CHANGE
We explain the termination of Ca release by IP unbinding in Subsec. III D. In this casetwo clusters are linked by the diffusion-like term as shown in eqs. (10), (11). The change ofthe [IP ] is also governed by the diffusion. To estimate the timescales of the IP change wesolve the diffusion equation in the infinite 2D domain (Fig. S3). The initial shape of the IP excitation ( blue line in Fig. S3) is assumed as Gaussian: p ( x, y,
0) = p a e − x − y σ , (S4)the amplitude is p a = 4 . µ M and variance σ = 1 µ m .The solution of the diffusion equation in infinite 2D domain is p ( x, y, t ) = Z + ∞−∞ Z + ∞−∞ p ( ξ, ζ , G ( x − ξ, y − ζ , t ) d ξ d ζ, (S5)where the Green’s function is G ( x − ξ, y − ζ , t ) = 14 πD p t e − ( x − ξ )2+( y − ζ )24 Dpt . (S6)4 p (cid:217)(cid:218)(cid:219) (cid:220)(cid:221) cl (cid:222)(cid:223)(cid:224) Æ(cid:226)ª(cid:228) (cid:229) r cl (cid:230)(cid:231)ŁØŒ º p p a p (cid:236) p (cid:237) FIG. S3. (Color online) The schematic representation of the IP diffusion in the infinite 2Ddomain containing two channels. By integrating Eq. (S5) we obtain the bottom case of the Eq. (15) in the main text p ( x, y, t ) = p a σ σ + 2 D p t e − x − y Dpt +2 σ , (S7)where D p = 10 µ m / s. The IP traces at each cluster (Fig. S3 at x = ± ( ℓ + r cl ) / y = 0)are shown as red dashed lines in Figs. 8 (d) and (e). [1] S. R¨udiger, Phys Rev E , 062717 (2014)., 062717 (2014).