A master stability function approach to cardiac alternans
Yi Ming Lai ORCID: orcid.org/0000-0002-5397-87531,
Joshua Veasy1,
Stephen Coombes1 &
Rüdiger Thul1
During a single heartbeat, muscle cells in the heart contract and relax. Under healthy conditions, the behaviour of these muscle cells is almost identical from one beat to the next. However, this regular rhythm can be disturbed giving rise to a variety of cardiac arrhythmias including cardiac alternans. Here, we focus on so-called microscopic calcium alternans and show how their complex spatial patterns can be understood with the help of the master stability function. Our work makes use of the fact that cardiac muscle cells can be conceptualised as a network of networks, and that calcium alternans correspond to an instability of the synchronous network state. In particular, we demonstrate how small changes in the coupling strength between network nodes can give rise to drastically different activity patterns in the network.
The heart consists of millions of muscle cells called cardiac myocytes (Bers 2002; Eisner et al. 2017). Upon electrical stimulation, cardiac myocytes first contract and then relax. What we perceive as a heartbeat is the coordinated contractile response of large numbers of cardiac myocytes, initiated by an electrical signal that travels across the heart muscle. The link between electrical stimulation and contraction lies in the dynamics of intracellular calcium (Ca 2+) (Bers 2002; Eisner et al. 2017). Essentially, electrical excitation leads to a transient rise of the intracellular Ca 2+ concentration, which in turn triggers contraction and subsequent relaxation of the cellular contractile machinery. Under healthy conditions, these cycles of electrical activity and Ca 2+ transients remain almost identical from heartbeat to heartbeat. However, molecular changes can induce irregular patterns (Qu and Weiss 2014; Qu et al. 2014; Karma 2013; Krogh-Madsen and Christini 2012; Landstrom et al. 2017). One of the earliest aberrations are so-called cardiac alternans, where the duration of the electrical signal and the maxima of the intracellular Ca 2+ concentration alternate in a large-small-large-small fashion (Alvarez-Lacalle et al. 2015; Weiss et al. 2006; Shiferaw et al. 2003; Cherry 2017; Tomek et al. 2018; Alvarez-Lacalle et al. 2013; Groenendaal et al. 2014; Shiferaw et al. 2005; Restrepo et al. 2008; Kanaporis and Blatter 2017; Edwards and Blatter 2014; Shkryl et al. 2012; Qu et al. 2016). While cardiac alternans are not life-threatening per se, they often form precursors to more severe if not fatal cardiac arrhythmias such as sudden cardiac death. Understanding the emergence and progression of cardiac alternans has hence been the focus of intense research.
Cardiac myocytes can be conceptualised as a network of networks (see Fig. 1). Each node in the network corresponds to a so-called calcium release unit (CRU). The majority of the molecular machinery that shapes the intracellular Ca 2+ transients is located here (Bers 2002; 2008). CRUs are coupled via Ca 2+ diffusion through the cytosol and the sarcoplasmic reticulum (SR), respectively, forming a large network. At the network level, the rhythm of a healthy cardiac myocyte corresponds to a synchronous network state, while Ca 2+ alternans emerge when a synchronous state loses stability.
Schematic of 3 coupled CRUs with network labels μ−1, μ and μ+1 showing the regions of the the different Ca 2+ concentrations: subsarcolemmal space (\(c_{\mathrm {s}}^{\mu }\)), bulk cytsosol (\(c_{\mathrm {i}}^{\mu }\)), unrecruited SR (\(c_{\mathrm {u}}^{\mu })\) and total Ca 2+ concentration in the SR (\(c_{\mathrm {j}}^{\mu }\)). The orange bidirectional arrows represent Ca 2+ diffusion though the bulk cytosol with time constant τc and through the SR with time constant τsr, respectively
In the present study, we investigate the linear stability of the synchronous network state as the strength of Ca2+ diffusion is varied. Due to different intracellular morphologies and biochemical compositions, Ca 2+ generally diffuses more quickly in the cytosol than in the SR (but see Petersen et al. (2017) for a different view). While diffusion of free Ca 2+ in the cytosol has been estimated to be 223 μm2s−1 (Allbritton et al. 1992), Ca 2+ buffers reduce this value substantially (Smith et al. 1996; Wagner and Keizer 1994). The strength of Ca 2+ diffusion in the SR has been controversial for more than a decade, and the verdict of whether it is fast or slow is still out (Swietach et al. 2008; 2010; Picht et al. 2011; Bers and Shannon 2013).
Our particular interest is in the emergent network patterns just after the onset of a synchronous instability. This corresponds to recently discovered microscopic Ca 2+ alternans (Tian et al. 2012). Here, the global Ca 2+ signal, i.e. the Ca 2+ response averaged across an entire cardiac myocyte, looks healthy, while the dynamics of single CRUs is irregular. This interplay between macroscopic Ca 2+ signals that look physiologically healthy and pathological local Ca 2+ signals is interpreted as the earliest onset of Ca 2+ alternans and the first sign that the healthy synchronous network behaviour has lost stability.
We have recently shown that changing Ca2+ diffusion in a network of CRUs leads to two kinds of network instabilities (Veasy et al. 2019). If cytosolic Ca 2+ diffusion dominates, the network undergoes the traditional period doubling bifurcation where each node follows a period-2 orbit with an alternating pattern of large-small-large peak Ca 2+ amplitudes. At the same time, neighbouring CRUs are out-of-phase with one another: when the Ca 2+ transient is large at one CRU, the adjacent CRU displays a small amplitude Ca 2+ transient. On the other hand, if luminal Ca 2+ diffusion dominates, i.e. Ca 2+ diffusion in the SR is faster than in the cytosol, we find a saddle-node bifurcation at the network level. In this case, each CRU follows a period-1 orbit, but the peak amplitudes of neighbouring CRUs alternate. This means that the global Ca 2+ signal is almost identical for Ca 2+ alternans emerging through either bifurcation, but the local dynamics is distinct. Consequently, microscopic Ca 2+ alternans may possess a much richer pattern space than previously thought.
To further unravel the complexity of microscopic Ca2+ alternans, we here compute the master stability function (MSF) for the network (Pecora and Carroll 1998). This approach has been instrumental in understanding instabilities of synchronous network states and has recently been generalised to more structured network dynamics such as cluster states and to the case of nearly-identical oscillators (Lai et al. 2018; Coombes et al. 2018; Coombes and Thul 2016; Ladenbauer et al. 2013; Sun et al. 2009; Pecora et al. 2000). A key input for the computation of the MSF is the synchronous network state and its Jacobian. Since for diffusively coupled nodes, the synchronous network state corresponds to the periodic solution of a single CRU, the mathematical tractability of the MSF significantly depends on the mathematical structure of the ordinary differential equations (ODEs) that describe the behaviour of a CRU. Traditionally, the dynamics of CRUs is governed by coupled nonlinear ODEs, which can only be solved numerically. This precludes any explicit construction of the MSF. To make progress here, we employ a piecewise linear (PWL) caricature (Thul and Coombes 2010) of a well established Ca 2+ cycling model (Shiferaw et al. 2003) for a single CRU. This allows for the explicit construction of the MSF, which is key for the results presented here. In particular, we employ the MSF to explain non-intuitive abrupt changes in the patterns of microscopic Ca 2+ alternans. The results from our theory are in excellent agreement with direct numerical simulations, illustrating the predictive power of our approach and the benefits of PWL models. Our findings also highlight that Ca 2+ diffusion exerts a different effect on the network dynamics of cardiac myocytes depending on whether it occurs predominantly in the cytosol or the SR.
Figure 1 shows a schematic of 3 nodes in the network. For each node, indexed by the label μ, we distinguish between 4 different Ca 2+ concentrations: the Ca 2+ concentration in the subsarcolemmal space (\(c_{\mathrm {s}}^{\mu }\)), the bulk Ca 2+ concentration (\(c_{\mathrm {i}}^{\mu }\)), the Ca 2+ concentration in the unrecruited SR (\(c_{\mathrm {u}}^{\mu })\) and the total Ca 2+ concentration (\(c_{\mathrm {j}}^{\mu }\)).
The ODEs for a single CRU read as
$$\begin{array}{*{20}l} \frac{\text{d}c_{\mathrm{s}}^{\mu}}{\text{d} t} &= \beta_{\mathrm{s}}\left[\frac{v_{\mathrm{i}}}{v_{\mathrm{s}}} \left(I_{\mathrm{r}}^{\mu}-\frac{c_{\mathrm{s}}^{\mu}-c_{\mathrm{i}}^{\mu}}{\tau_{\mathrm{s}}}-I_{\text{CaL}}^{\mu} \right) +I_{\text{NCX}}^{\mu}\right]\, {,} \end{array} $$
$$\begin{array}{*{20}l} ~ \frac{\text{d}c_{\mathrm{i}}^{\mu}}{\text{d}t} & = \beta_{\mathrm{i}} \left[\frac{c_{\mathrm{s}}^{\mu}-c_{\mathrm{i}}^{\mu}}{\tau_{\mathrm{s}}}-I_{\text{up}}^{\mu}\right]+\sum_{\epsilon\in{\mathcal{I}}_{n}} \frac{c_{\mathrm{i}}^{\epsilon} -c_{\mathrm{i}}^{\mu}}{\tau_{\mathrm{c}}}\, {,} \end{array} $$
$$\begin{array}{*{20}l} ~ \frac{\text{d}c_{\mathrm{j}}^{\mu}}{\text{d}t} & = -I_{\mathrm{r}}^{\mu}+I_{\text{up}}^{\mu}+\sum_{\epsilon\in{\mathcal{I}}_{n}}\frac{c_{\mathrm{j}}^{\epsilon}-c_{\mathrm{j}}^{\mu}}{\tau_{\text{sr}}}\, {,} \end{array} $$
$$\begin{array}{*{20}l} ~ \frac{\text{d}c_{\mathrm{u}}^{\mu}}{\text{d}t} & = \frac{c_{\mathrm{j}}^{\mu}-c_{\mathrm{u}}^{\mu}}{\tau_{\mathrm{a}}}\, {,} \end{array} $$
(1d)
$$\begin{array}{*{20}l} ~ \frac{\text{d}I_{\mathrm{r}}^{\mu}}{\text{d} t} & = -g I_{\text{CaL}}Q(c_{\mathrm{u}}^{\mu})-\frac{I_{\mathrm{r}}^{\mu}}{\tau_{\mathrm{r}}}\, {.} \end{array} $$
(1e)
The last ODE captures the Ca2+ release current from the SR into the subscarolemmal space and depends on the Ca2+ concentration in the unrecruited SR, \(c_{\mathrm {u}}^{\mu }\). We model the currents \(I_{\text {CaL}}^{\mu }\), \(I_{\text {NCX}}^{\mu }\) and \(I_{\text {up}}^{\mu }\) as PWL functions, which renders Eq. (1) a PWL model (see Appendix). For a detailed discussion of Eq. (1) together with its biological interpretation, we refer the reader to (Shiferaw et al. 2003; Thul and Coombes 2010). For the current study, the key terms are the coupling functions in Eqs. (1b) and (1c), which correspond to the sums over ε in each ODE. The set \({\mathcal {I}}_{n}\) indexes the nearest neighbours of the μth node. The linear differences in the coupling terms represent a discrete form of diffusion with time scales τc and τsr in the cytosol and SR, respectively.
For a single node, the PWL nature of Eq. (1) means that there are m switching manifolds, between which the dynamics can be written as
$$ \frac{\text{d} x^{\mu}}{\text{d} t}=A_{i} x^{\mu} +f(t)\, {,} $$
where \(x^{\mu }=(c_{\mathrm {s}}^{\mu }, c_{\mathrm {i}}^{\mu }, c_{\mathrm {j}}^{\mu }, c_{\mathrm {u}}^{\mu }, I_{\mathrm {r}}^{\mu })\), \(A_{i} \in \mathbb {R}^{5 \times 5}\), i=1,…,m, is constant and \(f(t) \in \mathbb {R}^{5}\) collects all explicitly time-dependent functions that describe the electrical activity in the model. For the analysis presented here, we assume that f is periodic with a period of Tp. This reflects a common practice in cardiac research whereby cardiac myocytes are paced by an external stimulus with period Tp while recording the intracellular Ca 2+ concentration, see e.g. Shiferaw et al. (2003). At the network level with N nodes, Eq. (1) takes the compact form
$$ \frac{\text{d}{x}}{\text{d}{t}} x t = A x + F(t) + \sigma G \otimes H x\, {,} $$
where \(x=(x^{1}, x^{2}, \ldots, x^{N}) \in \mathbb {R}^{5N}\), \(A \in \mathbb {R}^{5N\times 5N}\) and F(t)=1N⊗f(t). Here, ⊗ denotes the standard tensor product and \(1_{N} \in \mathbb {R}^{N}\) is a column vector containing only 1s. Note that A is always constant between switching events and is block diagonal with entries taken from the set of Ai. \(G \in \mathbb N^{N \times N}\) refers to the graph Laplacian of the network, and \(H \in \mathbb R^{5 \times 5}\) encodes through which variables the coupling occurs. For instance, for cytosolic coupling only, H2,2=1, while all other components of H vanish. The overall coupling strength is given by σ, which e.g. in the case of pure cytosolic coupling is \(\tau _{\mathrm {c}}^{-1}\).
To ascertain the linear stability of the synchronous network state s(t) where s(t)=x1(t)=…=xN(t), we introduce small network perturbations δx via \(\widetilde x(t)=s(t)+\delta x(t)\), where \(\widetilde x(t)\) corresponds to the perturbed network state. Since we perturb off the synchronous network state, we can assume that for the majority of time, all CRUs are described by the same matrix Ai. Therefore, linearising Eq. (3) and block-diagonalising it with the linear transformation ξ=(P⊗I5)−1δx, where P is the matrix of eigenvectors that diagonalises G and In is the n-dimensional identity matrix, we obtain
$$ \frac{d\xi}{dt} = \left [I_{N} \otimes A_{i} - \sigma \Lambda \otimes H \right ]\xi \, {.} $$
Here Λ is diagonal holding the eigenvalues λi of G, i.e. GP=PΛ. Because both IN and Λ are diagonal, Eq. (4) is block diagonal in \(\xi ^{\mu } \in \mathbb {R}^{5}\) via
$$ \frac{\text{d} \xi^{\mu}}{\text{d} t} = \left [A_{i} - \sigma \lambda_{\mu} H \right ]\xi^{\mu} \, {.} $$
In other words, the linear stability problem for the full 5N-dimensional network can be decomposed into N 5-dimensional problems parametrised by the eigenvalue λμ. If the Floquet multipliers for the solutions of each ξμ lie in the unit disk, the synchronous network state is stable, otherwise, synchrony is unstable.
Because the dynamics is PWL and the vector field is continuous at the switching manifolds, we can immediately solve Eq. (5) using matrix exponentials. Let Δ denote the period of the synchronous state, Δi the time-of-flight during which the dynamics is governed by Ai, i.e. \(\Delta =\sum _{i} \Delta _{i}\), and ξ(0) some initial perturbation. Then the perturbation after one period is given by ξμ(Δ)=Γ(σλμ)ξμ(0) with
$$ \Gamma(\sigma \lambda_{\mu})=\exp\left[(A_{N} - \sigma \lambda_{\mu} H)\Delta_{N} \right]\cdots \exp\left[(A_{1} - \sigma \lambda_{\mu} H)\Delta_{1} \right]\, {.} $$
Since the above derivation holds for any graph Laplacian, it is instructive to replace σλμ with \(\eta \in \mathbb {C}\) in Eq. (6), noting that the eigenvalues of the graph Laplacian are generally complex. The MSF is then the function that maps η to the largest real part of the Floquet exponents associated with Γ(η). Put differently, if q(η) is an eigenvalue of Γ(η), then the MSF returns κ(η)= maxqRe{log(q(η)}/Δ. If κ(σλi)<0 for all λi, then synchrony is stable, otherwise it is unstable. Therefore, the MSF can be computed independently of the choice of network and then used to assess linear stability of the synchronous network state for a particular network.
We first compute the MSF for pure cytosolic coupling. In this case, σ−1=τc, \(\tau _{\text {sr}}^{-1}=0\) and the only non-vanishing component of H is H2,2. Figure 2a shows the zero-contour of the MSF, where we observe that the MSF is negative inside the ellipse and positive in the remainder of the complex plane. Hence, if every σλi falls inside the ellipse, the synchronous network state is stable, otherwise, it is unstable. Because diffusive coupling is symmetric and the coupling strength is positive, the eigenvalues λi of the graph Laplacian are real and negative including a zero eigenvalue. The latter corresponds to the periodic orbit of an uncoupled node, which entails that a necessary condition for the existence of a stable synchronous network state is that the periodic orbit of an uncoupled node is linearly stable. By taking the limit σ→0, all σλi can be contained within the ellipse, indicating that the synchronous network state is linearly stable for weak coupling. By increasing the coupling strength, the eigenvalue with the most negative real part exits the ellipse on the left, rendering synchrony unstable. The emergent network state is depicted in Fig. 2b and c, which show the peak Ca 2+ concentration in the bulk cytosol on successive pacing periods. The alternations of yellow and blue in each figure indicate that neighbouring CRUs exhibit alternating values of the Ca 2+ amplitudes. When comparing Fig. 2b and c, we find that when the Ca 2+ transient is small during the first pacing period it is larger during the second pacing period, and vice versa. This pattern of network activity is consistent with the MSF where the critical eigenvalue qc leaves the ellipse with cos(argqc)=−1, indicating a period-doubling bifurcation at the network level. The eigenvector corresponding to the eigenvalue that has crossed the stability boundary is plotted in Fig. 2d and agrees very well with the results shown in Fig. 2b and c.
a Zero-contour of the MSF for cytosolic coupling and Tp=0.9. The MSF is negative inside the ellipse, labeled by S, and positive outside, denoted by U. The colour represents the value of cos(arg(q(η))). b, c Peak Ca 2+ concentration in the bulk on successive beats. d Eigenvector corresponding to the single eigenvalue λk for which η=σλk lies outside the ellipse. Here, σ=0.225
We now contrast the results for cytosolic coupling with those for purely luminal coupling. Here, σ−1=τsr, \(\tau _{\mathrm {c}}^{-1}=0\) and the only non-zero component of H is H3,3. As Fig. 3a illustrates, the topology of the MSF has changed significantly. There are now two regions of stability in the complex plane separated by a region of instability. Because instabilities can only occur when eigenvalues q move along the negative part of the real axis, we can characterise Fig. 3a by taking a cut along the negative real axis. Figure 3b summarises the resultant regions of stability and instability for different values of Tp. As we increase Tp, the region of instability shrinks, up to a point when synchrony is always stable. This is consistent with experimental findings that show that cardiac myocytes do not undergo instabilities when stimulated at sufficiently low frequencies. While the network instability for pure cytosolic coupling occurs via a period-doubling bifurcation, the colour map in Fig. 3a reveals that the network loses stability through a saddle-node bifurcation where the critical eigenvalue leaves the unit disk through +1.
a Zero-contours of the MSF for luminal coupling and Tp=0.9. The MSF is negative in regions denoted by S and positive in regions labeled U. The colour represents the value of cos(arg(q(η))). b Zero values of the MSF for real values of η as pacing periods vary
In Fig. 3b, we plotted the stability regions as a function of the general MSF parameter η and the pacing period Tp. In a practical application, where the network structure is fixed, it is more natural to examine stability as a function of the coupling strength σ and the the pacing period Tp. Figure 4 provides an illustration of this. The green line at bottom indicates the critical pacing period when the period-1 orbit of an isolated CRU goes unstable. For periods faster than this, an isolated CRU displays Ca 2+ alternans, which then feed forward to the network level. Since we are interested in how coupling between CRUs induces instabilities, we restrict our attention to pacing periods above the green line. Note that we plot \(\sigma _{\text {SR}}^{-1}\) on the x-axis, so that weak coupling corresponds to the left part and strong coupling to the right part of the figure, respectively. For constant Tp, the synchronous network state is stable for very small and large coupling, while it is unstable for intermediate coupling strengths. For fixed values of \(\sigma _{\text {SR}}^{-1}\), faster pacing periods generally destabilise solutions, which mirrors experimentally observed behaviour. The most prominent feature of Fig. 4a is a series of small bumps in the stability line for larger values of Tp. When we fix Tp and vary \(\sigma _{\text {SR}}^{-1}\) we observe that the spatial patterns of the emergent network solutions vary drastically as we cross from one 'instability bump' to the next, see Fig. 5. While for weaker coupling, i.e. smaller values of \(\sigma _{\text {SR}}^{-1}\), the Ca 2+ concentration exhibits a multi-modal distribution with peaks in the corner on one side and in the middle on the opposite side, respectively (Fig. 5a), stronger coupling leads to stripes of the intracellular Ca 2+ concentration (Fig. 5b)
a Stability region of the synchronous network state as a function of the luminal coupling strength \(\sigma _{\text {SR}}^{-1}\) and the pacing period Tp for a network of 4 CRUs. The synchronous state is unstable in the region labeled U and stable in the region denoted by S. The green line indicates the pacing period at which the period-1 orbit of a single CRU loses stability. b Stability regions of individual eigenvectors. The corresponding eigenvectors are shown in blue, and the corresponding eigenvalues are listed. The green line is the same as in (a)
Peak values of the bulk Ca2+ concentration during one pacing period for σSR=3.2 (a), σSR=4.1 (b) in a 10×25 network of CRUs at Tp=1.05
We can explain these sudden changes in the activity patterns of the intracellular Ca2+ concentration by starting with Fig. 4b. This figure shows that the central region of instability seen in Fig. 4a is in fact a superposition of instability regions associated with different eigenvalues. For illustrative purposes, we computed Fig. 4b for a network of only 4 nodes. However, the main features of the diagram remain unchanged as we increase the number of nodes in the network. Since one of the eigenvalues always vanishes (see above), the three non-trivial eigenvalues are shown. Associated with each of these is a distinct region of instability. For instance, the region furthest to the left belongs to λ1, while the region to the right is controlled by λ3.
To compute these regions, we make use of the MSF as shown in Fig. 3a with zooms provided in Fig. 6a and b. Note that these correspond to a fixed pacing period and hence map onto one horizontal line in Fig. 4b. Figure 6a shows the case when the instability is driven by a single eigenvalue, say λk. The corresponding argument for the MSF, i.e. η=σλk is indicated by the red circle. As we change σ, η traces the space between the two instability lines. By computing the critical values of σ such that η intersects with the instability lines, we determine the left and right boundaries of the instability regions for a fixed eigenvalue in Fig. 4b. Figure 6b illustrates that upon changing σ, a different eigenvalue compared to the one in Fig. 6a shapes the instability. Note that there is again only one eigenvalue that is responsible for the instability. This corresponds exactly to the case when we increase σ in Fig. 4b for larger values of Tp and move from the 'bump' for λ1 to the 'bump' for λ2. Because each eigenvalue is associated with a specific eigenvector, the abrupt changes in the emergent network patterns reflect the often considerable variations among eigenvectors. In Fig. 4b, this can be seen by inspecting the eigenvectors that correspond to the respective eigenvalues. For a much larger network, this is confirmed by comparing the eigenvectors plotted in Fig. 6c and d with the simulation results in Fig. 5. There is excellent agreement between them. In Fig. 6e, we summarise the mechanism that gives rise to the different Ca 2+ activity patterns in the network. For weak coupling, only one eigenvalue (λ1) drives the instability. Upon increasing the coupling, both σλ1 and σλ2 move towards the left (recall that all eigenvalues are negative). Therefore, for some values of σ, synchrony is stable again. However, a further increase in the coupling strength causes σλ2 to move into the instability region, giving rise to a different emergent network state compared to the one for weak coupling.
a,b Zoom of the MSF for dominant luminal coupling around the the real axis with values of η=σλk (crosses) superimposed for σSR=3.2 (a), σSR=4.1 (b) in a 10×25 network of CRUs paced at Tp=1.05. The value of η for which the MSF is positive is circled in red. (c,d) Eigenvectors corresponding to the critical values of η in (a) and (b), respectively. e Schematic of how changing the coupling strength σ can give rise to different patterns of the network activity. See text for details
Figure 3b already suggests that the regions where the MSF is negative strongly depends on the pacing period Tp. There, we focussed on the negative real axis for the MSF parameter η since η cannot be complex or positive for diffusive coupling. For a more detailed view, we now plot the zero-contours of the MSF in the complex plane as a function of Tp in Fig. 7a. This three-dimensional plot highlights that the topology of the zero-contour changes significantly as a function of Tp. For larger pacing periods (Fig. 7b), we find a connected region where the MSF is negative, which resembles half a bowtie on the right. As we decrease Tp, the narrow part of the bowtie contracts, until two disconnected regions emerge as exemplified by Fig. 3a. When we lower Tp even further, the stable central regions contracts along the real axis, but expands along the imaginary axis (see Fig. 7c). Here, the synchronous network state can lose stability via a saddle-node bifurcation, indicated by the yellow colour of one of the stability boundaries. For larger pacing periods as shown in Fig. 7b, we note that there is a period doubling bifurcation towards the right side of the central stable region. A similar line exists towards the right of the unstable region in Fig. 7c. However, in both cases, these period doubling bifurcations occur for Re(η)>0, which is not permissible for diffusive coupling.
a Zero-contour of the MSF for purely luminal coupling as a function of Tp. Zero-contour of the MSF for Tp=1.075s (b) and Tp=0.55s (c). The MSF is negative in regions labeled S and positive in regions denoted by U. The colour represents the value of cos(arg(q(η)))
So far, we have studied purely cytosolic and luminal coupling, respectively. Under physiologically realistic conditions, however, Ca 2+ diffuses through both the cytosol and the SR. We therefore computed the bifurcation lines in the (τc,τsr) plane. For this, we replace σH in Eq. (3) with a single matrix H whose entries are all zero except for \(H_{2,2}=\tau _{\mathrm {c}}^{-1}\) and \(H_{3,3}=\tau _{\text {sr}}^{-1}\). To compute the MSF, we further introduce the general MSF parameters ηc and ηsr, which run along the axis in Fig. 8. Note that we can restrict the values of ηc and ηsr to negative real values since diffusion does not lead to complex eigenvalues of the graph Laplacian. For the computation of the MSF, we replace the matrices (Ak−ηH), k=1,…,N, in Eq. (6) with (A−H(ηc,ηsr)). Here, H(ηc,ηsr) is a matrix whose non-zero elements are H2,2=ηc and H3,3=ηsr. As Fig. 8 shows, the MSF is negative for large ranges of ηc and ηsr. When luminal coupling is weak, i.e. when the absolute value of ηsr is small, we only observe a period-doubling bifurcation upon variation of ηc. This is indicated by the (−1) line. In contrast, when cytosolic coupling is negligible, the network instability occurs via a saddle-node bifurcation upon variation of ηsr, which we mark by the (+1) line. These findings illustrate that the network is generally stable when coupling is balanced, i.e. when ηc and ηsr are of similar magnitude. However, when one coupling dominates, we find either a period-doubling or a saddle-node bifurcation.
Zero-contour of the MSF in the (ηc,ηsr) plane for Tp=0.6. Note that both ηc and ηsr are real since we focus on diffusive couling. The MSF is negative in the region labeled S and positive in regions denoted by U. The green curve corresponds to saddle-node bifurcations (+1), while the blue curve refers to period-doubling bifurcations (−1)
Networks are ubiquitous in biology, and intracellular signalling cascades constitute a prime example. In the present study, we investigated the Ca 2+ dynamics in cardiac myocytes. On the one hand, this is highly relevant for our general well-being as disturbances in these networks are associated with numerous pathologies. On the other hand, Ca 2+ signalling in cardiac myocytes exemplifies a network of networks. Each node in the network corresponds to a CRU, whose dynamics in turn is governed by its own reaction network. A key aspect is that CRUs are coupled through two different channels: Ca 2+ diffusion in the cytosol and Ca 2+ diffusion in the SR.
Our particular interest is in understanding how different coupling strengths shape the synchronous network state. The reason for this is that loss of synchrony in a cardiac Ca 2+ network is associated with the emergence of cardiac Ca 2+ alternans. These constitute one of the earliest cardiac arrhythmias and act as precursors to more severe cardiac abnormalities including sudden cardiac death. Until recently, Ca 2+ alternans could only be observed at an advanced stage, i.e. when the Ca 2+ concentration oscillates out-of-phase in mesoscopic parts of a myocyte. However, improved imaging techniques now allow the recording of microscopic Ca 2+ alternans (Tian et al. 2012). Here, the Ca 2+ concentration averaged across the cell suggests a healthy cardiac myocyte, while in fact single CRUs may already display pathological Ca 2+ alternans. From a conceptual point of view, microscopic Ca 2+ alternans correspond to the pattern that emerges when the synchronous network state has just lost stability.
We recently reported that microscopic Ca2+ alternans can emerge via two different mechanisms: a traditional period-doubling bifurcation and a novel saddle-node bifurcation (Veasy et al. 2019). Strikingly, the emergent network patterns of Ca 2+ activity vary substantially as we move along the stability boundaries. To understand this, we here computed the MSF for the CRU network.
A comparison of Figs. 2a and 3a reveals that the MSF differs significantly between purely cytosolic and luminal coupling. While the zero-contour takes on the shape of an ellipse in the former, there are multiple zero-contours in the latter, delineating distinct regions where the MSF is negative. This has direct implications for the stability of the synchronous network state. As Fig. 6 illustrates, we can understand changes in stability when the MSF parameter σλ moves along the real line. By increasing σ, the value of σλk becomes more negative for all k. Hence, for purely cytosolic coupling, once one eigenvalue leads to a value of η outside the ellipse, synchrony is unstable. On the contrary, for luminal coupling, the small region where the MSF is positive means that increasing σ induces a sequence in which the synchronous network states alternates between stable and unstable. These alternations are responsible for the 'stability bumps' in Fig. 4.
Because each of these bumps are linked to a different eigenvalue — and correspondingly with a different eigenvector — we observe the abrupt changes in the emergent network activity when the synchronous network state loses stability. Figure 4b shows the explicit sequence of eigenvectors for a network of 4 nodes. Similar behaviour is seen in much larger networks, too. Figure 5 illustrates this for a bigger network. As the network size grows, the pattern space of the eigenvectors becomes richer, which also means that the patterns that can occur in the network exhibit more features. The above argument rests on the assumption that eigenvectors are good predictors of the emergent network state. This is generally true when only one eigenvector gives rise to the instability as is the case for Fig. 6a and b. Then, the associated eigenvectors (Fig. 6c and d) match very well with the results from direct numerical simulations shown in Fig. 5. Once σλ is positive for more than one eigenvalue, the direct predictive power of eigenvectors is reduced, as now the emergent network state is a linear superposition of several eigenvectors, but the weights are not known a priori.
As our results demonstrate, knowledge of the MSF is key for understanding the non-intuitive behaviour of the CRU network. Unfortunately, the MSF can often only be obtained numerically, which can be computationally expensive. We here make progress in this direction by employing a PWL model for the dynamics of a single CRU, which allows us to explicitly construct the synchronous network state and the map that propagates network perturbations. We can then compute the MSF in a semi-analytical manner, which makes it possible to produce three-dimensional visualisations as those depicted in Fig. 7. Based on plots like this, we can infer how the stability of synchrony in the CRU network varies as the pacing period is altered (which is a common experimental practice).
Since the computation of the MSF rests on the PWL nature of the CRU model, it is worth asking how well this approximation describes the dynamics of the full nonlinear model. As shown in Lai et al. (2019), the PWL model captures the core dynamics of the nonlinear model very well. Of course, there are parameter regimes of the nonlinear model that cannot be captured with the current parametrisation of the PWL model used here. For example, the L-type Ca 2+ current in the present study is either zero or takes on a constant non-zero value iCaL. If the nonlinear shape of the L-type Ca 2+ current is central to a study, it requires either the construction of a piecewise constant L-type Ca 2+ current with multiple levels, or one needs to resort to the nonlinear model. Similar considerations apply to the closure of the L-type Ca 2+ channel. We currently include voltage-dependent inactivation only. However, Ca 2+-dependent inactivation also exists (Josephson et al. 2010; Grandi et al. 2010). We can again amend the PWL model used here with an additional switch for Ca 2+-dependent inactivation, or if the interplay between the timescales of the inactivation processes becomes important, one might have to consider the nonlinear model. Furthermore, care needs to be taken when Ca 2+ buffers are the focus of attention. The current model treats buffer contributions as constant, but given the dynamic nature of the intracellular Ca 2+ concentration, the fraction of Ca 2+-bound buffers changes over time. In this case, the PWL model cannot be amended and the nonlinear model is the only choice. In case the PWL model can be tailored to the question at hand, the study in Rotstein et al. (2012) provides a conceptual blueprint for it. It shows an extension of the classical 3-piece approximation of the Fitzugh-Nagumo model (FitzHugh 1961; Nagumo et al. 1962) to investigate canard-like solutions.
As the focus of this study lies on Ca2+ alternans, we only consider diffusive, i.e. nearest-neighbour coupling. Consequently, the MSF parameter η is always real and negative. This only leads to period-doubling or saddle-node bifurcations, respectively, where cos(argq(η))=±1. However, the MSF provides information for arbitrary values of η and hence arbitrary network topologies. The line colours in Figs. 3a and 7c show that if η crosses the zero-contour of the MSF at positions different than the real axis, cos(argq(η))≠±1. This corresponds to a Neimark-Sacker bifurcation at the network level, and it will be interesting to explore the emergent network patterns.
In conclusion, a combination of PWL modelling and MSF techniques facilitated a detailed investigation of microscopic Ca2+ alternans in a network of CRUs. Crucially, our results explain the previously reported abrupt variations in network activity as the coupling strength in the network changes (Veasy et al. 2019). Moreover, our findings demonstrate that depending on whether Ca 2+ diffusion is stronger in the cytosol or in the SR, different microscopic Ca 2+ alternans emerge, with each mode of diffusion giving rise to distinct network patterns of the intracellular Ca 2+ concentration. While these findings may have implications for cardiac health, they also highlight on a more fundamental level that cell signalling more generally may be usefully conceptualised as a network of networks.
Here, we provide details of the currents, the load release function and the clamped voltage used in Eq. (1). For a further discussion, we refer the reader to Thul and Coombes (2010).
The L-type Ca2+ current is given by ICaL=Θ(V−VL)iCaL with the threshold voltage VL=Vmax−1 and a constant conductance
$$ i_{\text{CaL}}= -\,4\,\overline{i}_{\text{Ca}} P_{\text{Ca}} \frac{a_{\text{CaL}} F \gamma_{\mathrm{o}}\text{Ca}_{o}}{\text{exp}(2a_{\text{CaL}})-1}\, {,} $$
where aCaL=VmaxF/RT. The current through the NCX is modelled as \(I_{\text {NaCa}} = \phi (V)-\psi (V)c_{\mathrm {s}}^{\mu }\), where
$$\begin{array}{*{20}l} \phi &= \overline{I}_{\text{NaCa}} \frac{\eta_{q} \text{Na}_{\mathrm{i}}^{3} \text{Ca}_{\mathrm{o}}}{(K_{\text{mNa}}^{3}+\text{Na}_{\mathrm{o}}^{3}) (K_{\text{mCa}}+\text{Ca}_{\mathrm{o}})}\, {,} \end{array} $$
$$\begin{array}{*{20}l} \psi &= \overline{I}_{\text{NaCa}} \frac{\gamma_{\text{NaCa}} \text{Na}_{\mathrm{o}}^{3} \times 10^{-3}}{(K_{\text{mNa}}^{3}+\text{Na}_{\mathrm{o}}^{3}) (K_{\text{mCa}}+\text{Ca}_{\mathrm{o}})}\, {.} \end{array} $$
Here, we introduce the piecewise constant function
$$ \begin{aligned} \gamma_{\text{NaCa}}= \left\{\begin{array}{ll} 0.45\, {,} & V>V_{\text{NaCa}}\, {,}\\ 4\, {,} & V \leq V_{\text{NaCa}}\, {,} \\ \end{array}\right. \end{aligned} $$
which switches at a threshold voltage of VNaCa=−50 mV. The function ηq is given by ηq=0.0501α2+0.3816α+0.9182 with α=FV/RT. Ca2+ uptake takes the form \(I_{\text {up}}=v_{\text {up}} c_{\mathrm {i}}^{\mu }\). The load-release function \(Q=Q(c_{\mathrm {u}}^{\mu })\) is described by the PWL function
$$ \begin{aligned} Q(c_{\mathrm{u}}^{\mu})=10^{-3} \left\{\begin{array}{ll} 0\, {,} & 0 \leq c_{\mathrm{u}}^{\mu} < 50\, {,}\\ c_{\mathrm{u}}^{\mu}-50\, {,} & 50 \leq c_{\mathrm{u}}^{\mu} < 115\, {,}\\ k c_{\mathrm{u}}^{\mu}+s\, {,} & c_{\mathrm{u}}^{\mu} \geq 115\, {,} \end{array}\right. \end{aligned} $$
where k measures the steep nonlinear dependence of Ca2+ release on the SR Ca2+ concentration, and s is a constant that is chosen such that Q is continuous. We model the clamped voltage for a pacing period Tp as
$$ \begin{aligned} V(t)= \left\{\begin{array}{ll} V_{+}(t)\, {,} & k T_{p} \leq t \leq (k+x) T_{p}\, {,}\\ V_{\text{min}}\, {,} & (k+x) T_{p} \leq t < (k+1) T_{p}\, {,} \end{array}\right. \end{aligned} $$
where \(k \in \mathbb N\) counts the number of APs and x=ax/(ax+Tp) with ax=2/3. The resting potential is given by Vmin=−70mV, and V+(t) captures the shape of the clamped voltage, given by
$$ V_{+}(t)=V_{\text{min}}+(V_{\text{max}}-V_{\text{min}}) \sqrt{1-\left (\frac{t-kT_{p}}{xT_{p}} \right)^{2}}\, {,} $$
for kTp≤t≤(k+x)Tp, where the maximal AP is given by Vmax=30mV. Note that since V is an explicit function of time, all variables that only depend on V are explicitly time-dependent. We collect these time-dependent functions in the function f(t) used in Eq. (2). The standard set of parameter values used in this study are listed in Table 1
Table 1 Standard parameter values used in the study
Availability of data and material
Ca2+ :
CRU:
Calcium release unit
MSF:
Master stability function
ODE:
Ordinary differential equation
PWL:
Piecewise linear
SR:
Sacroplasmic reticulum
Allbritton, NL, Meyer T, Stryer L (1992) Range of messenger action of calcium ion and inositol 1,4,5-trisphosphate. Sci (New York, N.Y.) 258(5089):1812–1815.
Alvarez-Lacalle, E, Cantalapiedra IR, Peñaranda A, Cinca J, Hove-Madsen L, Echebarria B (2013) Dependency of calcium alternans on ryanodine receptor refractoriness. PLoS ONE 8(2):55042.
Alvarez-Lacalle, E, Echebarria B, Spalding J, Shiferaw Y (2015) Calcium alternans is due to an order-disorder phase transition in cardiac cells. Phys Rev Lett 114(10):108101.
Bers, DM (2002) Cardiac excitation-contraction coupling. Nature 415(6868):198–205.
Bers, DM (2008) Calcium cycling and signaling in cardiac myocytes. Annu Rev Physiol 70:23–49.
Bers, DM, Shannon TR (2013) Calcium movements inside the sarcoplasmic reticulum of cardiac myocytes. J Mol Cell Cardiol 58(1):59–66.
Cherry, EM (2017) Distinguishing mechanisms for alternans in cardiac cells using constant-diastolic-interval pacing. Chaos 27(9):093902.
Coombes, S, Lai YM, Şayli M, Thul R (2018) Networks of piecewise linear neural mass models. Eur J Appl Math 9:1–22.
Coombes, S, Thul R (2016) Synchrony in networks of coupled non-smooth dynamical systems: Extending the master stability function. Eur J Appl Math 27(6):904–922.
Edwards, JN, Blatter LA (2014) Cardiac alternans and intracellular calcium cycling. Clin Exp Pharmacol Physiol 41(7):524–532.
Eisner, DA, Caldwell JL, Kistamás K, Trafford AW (2017) Calcium and Excitation-Contraction Coupling in the Heart. Circ Res 121(2):181–195.
FitzHugh, R (1961) Impulses and physiological states in theoretical models of nerve membrane. Biophys J 1(6):445–466.
Grandi, E, Morotti S, Ginsburg KS, Severi S, Bers DM (2010) Interplay of voltage and Ca-dependent inactivation of L-type Ca current. Prog Biophys Mol Biol 103(1):44–50.
Groenendaal, W, Ortega FA, Krogh-Madsen T, Christini DJ (2014) Voltage and Calcium Dynamics Both Underlie Cellular Alternans in Cardiac Myocytes. Biophys J 106(10):2222–2232.
Josephson, IR, Guia A, Lakatta EG, Lederer WJ, Stern MD (2010) Ca 2+-dependent components of inactivation of unitary cardiac L-type Ca 2+ channels. J Physiol 588(Pt 1):213–223.
Kanaporis, G, Blatter LA (2017) Alternans in atria: Mechanisms and clinical relevance. Medicina (Kaunas, Lithuania) 53(3):139–149.
Karma, A (2013) Physics of Cardiac Arrhythmogenesis. Ann Rev Condens Matter Phys 4:313–337.
Krogh-Madsen, T, Christini DJ (2012) Nonlinear dynamics in cardiology. Annu Rev Biomed Eng 14:179–203.
Ladenbauer, J, Lehnert J, Rankoohi H, Dahms T, Schöll E, Obermayer K (2013) Adaptation controls synchrony and cluster states of coupled threshold-model neurons. Phys Rev E 88(4):042713.
Lai, YM, Coombes S, Thul R (2019) Calcium buffers and L-type calcium channels as modulators of cardiac subcellular alternans. http://biorxiv.org/cgi/content/short/683011v1.
Lai, YM, Thul R, Coombes S (2018) Analysis of networks where discontinuities and nonsmooth dynamics collide: understanding synchrony. Eur Phys J Spec Top 227(10-11):1251–1265.
Landstrom, AP, Dobrev D, Wehrens XHT (2017) Calcium Signaling and Cardiac Arrhythmias. Circ Res 120(12):1969–1993.
Nagumo, J, Arimoto S, Yoshizawa S (1962) An Active Pulse Transmission Line Simulating Nerve Axon In: Proceedings of the IRE, 2061–2070.
Pecora, L, Carroll T (1998) Master stability functions for synchronized coupled systems. Phys Rev Lett 80(10):2109–2112.
Pecora, L, Carroll T, Johnson G, Doug M, Fink KS (2000) Synchronization stability in coupled oscillator arrays: Solution for arbitrary configurations. Int J Bifurcation Chaos Applied Sci Eng 10(2):273–290.
Petersen, OH, Courjaret R, Machaca K (2017) Ca 2+ tunnelling through the ER lumen as a mechanism for delivering Ca 2+ entering via store-operated Ca 2+ channels to specific target sites. J Physiol 595(10):2999–3014.
Picht, E, Zima AV, Shannon TR, Duncan AM, Blatter LA, Bers DM (2011) Dynamic calcium movement inside cardiac sarcoplasmic reticulum during release. Circ Res 108(7):847–856.
Qu, Z, Hu G, Garfinkel A, Weiss JN (2014) Nonlinear and stochastic dynamics in the heart. Phys Rep 543(2):61–162.
Qu, Z, Liu MB, Nivala M (2016) A unified theory of calcium alternans in ventricular myocytes. Sci Rep 6:35625.
Qu, Z, Weiss JN (2014) Mechanisms of Ventricular Arrhythmias: From Molecular Fluctuations to Electrical Turbulence. Annu Rev Physiol 77:29–55.
Restrepo, JG, Weiss JN, Karma A (2008) Calsequestrin-Mediated Mechanism for Cellular Calcium Transient Alternans. Biophys J 95(8):23–23.
Rotstein, HG, Coombes S, Gheorghe AM (2012) Canard-Like Explosion of Limit Cycles in Two-Dimensional Piecewise-Linear Models of FitzHugh–Nagumo Type. SIAM J Appl Dyn Syst 11(1):135–180.
Shiferaw, Y, Sato D, Karma A (2005) Coupled dynamics of voltage and calcium in paced cardiac cells. Phys Rev E 71(2 Pt 1):021903.
Shiferaw, Y, Watanabe MA, Garfinkel A, Weiss JN, Karma A (2003) Model of intracellular calcium cycling in ventricular myocytes. Biophys J 85(6):3666–3686.
Shkryl, VM, Maxwell JT, Domeier TL, Blatter LA (2012) Refractoriness of sarcoplasmic reticulum Ca 2+ release determines Ca 2+ alternans in atrial myocytes. Am J Physiol Heart Circ Physiol 302(11):2310–20.
Smith, GD, Wagner J, Keizer J (1996) Validity of the rapid buffering approximation near a point source of calcium ions. Biophys J 70(6):2527–2539.
Sun, J, Bollt EM, Nishikawa T (2009) Master stability functions for coupled nearly identical dynamical systems. Europhys Lett 85(6):60011.
Swietach, P, Spitzer KW, Vaughan-Jones RD (2008) Ca 2+-mobility in the sarcoplasmic reticulum of ventricular myocytes is low. Biophys J 95(3):1412–1427.
Swietach, P, Spitzer KW, Vaughan-Jones RD (2010) Modeling calcium waves in cardiac myocytes: importance of calcium diffusion. Front Biosci 15:661–680.
Tomek, J, Tomková M, Zhou X, Bub G, Rodriguez B (2018) Modulation of Cardiac Alternans by Altered Sarcoplasmic Reticulum Calcium Release: A Simulation Study. Front Physiol 9:1306.
Thul, R, Coombes S (2010) Understanding cardiac alternans: a piecewise linear modeling framework. Chaos 20(4):045102.
Tian, Q, Kaestner L, Lipp P (2012) Noise-free visualization of microscopic calcium signaling by pixel-wise fitting. Circ Res 111(1):17–27.
Veasy, J, Lai YM, Coombes S, Thul R (2019) Complex patterns of subcellular cardiac alternans. J Theor Biol 478:102–114.
Wagner, J, Keizer J (1994) Effects of Rapid Buffers on Ca 2+ Diffusion and Ca 2+ Oscillations. Biophys J 67(1):447–456.
Weiss, JN, Karma A, Shiferaw Y, Chen P-S, Garfinkel A, Qu Z (2006) From pulsus to pulseless: the saga of cardiac alternans. Circ Res 98(10):1244–1253.
This work was supported by the Engineering and Physical Science Research Council (grant number EP/P007031/1).
School of Mathematical Sciences, University of Nottingham, Nottingham, NG7 2RD, UK
Yi Ming Lai, Joshua Veasy, Stephen Coombes & Rüdiger Thul
Yi Ming Lai
Joshua Veasy
Stephen Coombes
Rüdiger Thul
YML, JV, SC, and RT contributed to the design and implementation of the research, to the analysis of the results, and to the writing of the manuscript. All authors read and approved the final manuscript.
Correspondence to Yi Ming Lai.
Lai, Y.M., Veasy, J., Coombes, S. et al. A master stability function approach to cardiac alternans. Appl Netw Sci 4, 90 (2019). https://doi.org/10.1007/s41109-019-0199-z
DOI: https://doi.org/10.1007/s41109-019-0199-z
Network instability