A master stability function approach to cardiac alternans

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 ²⁺) (Bers 2002; Eisner et al. 2017). Essentially, electrical excitation leads to a transient rise of the intracellular Ca ²⁺ concentration, which in turn triggers contraction and subsequent relaxation of the cellular contractile machinery. Under healthy conditions, these cycles of electrical activity and Ca ²⁺ 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 ²⁺ 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 ²⁺ transients is located here (Bers 2002; 2008). CRUs are coupled via Ca ²⁺ 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 ²⁺ 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 ²⁺ concentrations: subsarcolemmal space (\(c_{\mathrm {s}}^{\mu }\)), bulk cytsosol (\(c_{\mathrm {i}}^{\mu }\)), unrecruited SR (\(c_{\mathrm {u}}^{\mu })\) and total Ca ²⁺ concentration in the SR (\(c_{\mathrm {j}}^{\mu }\)). The orange bidirectional arrows represent Ca ²⁺ diffusion though the bulk cytosol with time constant τ_c and through the SR with time constant τ_sr, respectively

Full size image

In the present study, we investigate the linear stability of the synchronous network state as the strength of Ca²⁺ diffusion is varied. Due to different intracellular morphologies and biochemical compositions, Ca ²⁺ 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 ²⁺ in the cytosol has been estimated to be 223 μm²s⁻¹ (Allbritton et al. 1992), Ca ²⁺ buffers reduce this value substantially (Smith et al. 1996; Wagner and Keizer 1994). The strength of Ca ²⁺ 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 ²⁺ alternans (Tian et al. 2012). Here, the global Ca ²⁺ signal, i.e. the Ca ²⁺ response averaged across an entire cardiac myocyte, looks healthy, while the dynamics of single CRUs is irregular. This interplay between macroscopic Ca ²⁺ signals that look physiologically healthy and pathological local Ca ²⁺ signals is interpreted as the earliest onset of Ca ²⁺ alternans and the first sign that the healthy synchronous network behaviour has lost stability.

We have recently shown that changing Ca²⁺ diffusion in a network of CRUs leads to two kinds of network instabilities (Veasy et al. 2019). If cytosolic Ca ²⁺ 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 ²⁺ amplitudes. At the same time, neighbouring CRUs are out-of-phase with one another: when the Ca ²⁺ transient is large at one CRU, the adjacent CRU displays a small amplitude Ca ²⁺ transient. On the other hand, if luminal Ca ²⁺ diffusion dominates, i.e. Ca ²⁺ 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 ²⁺ signal is almost identical for Ca ²⁺ alternans emerging through either bifurcation, but the local dynamics is distinct. Consequently, microscopic Ca ²⁺ alternans may possess a much richer pattern space than previously thought.

To further unravel the complexity of microscopic Ca²⁺ 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 ²⁺ 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 ²⁺ 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 ²⁺ 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 ²⁺ concentrations: the Ca ²⁺ concentration in the subsarcolemmal space (\(c_{\mathrm {s}}^{\mu }\)), the bulk Ca ²⁺ concentration (\(c_{\mathrm {i}}^{\mu }\)), the Ca ²⁺ concentration in the unrecruited SR (\(c_{\mathrm {u}}^{\mu })\) and the total Ca ²⁺ concentration (\(c_{\mathrm {j}}^{\mu }\)).

The ODEs for a single CRU read as

The last ODE captures the Ca²⁺ release current from the SR into the subscarolemmal space and depends on the Ca²⁺ 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

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 T_p. This reflects a common practice in cardiac research whereby cardiac myocytes are paced by an external stimulus with period T_p while recording the intracellular Ca ²⁺ concentration, see e.g. Shiferaw et al. (2003). At the network level with N nodes, Eq. (1) takes the compact form

where \(x=(x^{1}, x^{2}, \ldots, x^{N}) \in \mathbb {R}^{5N}\), \(A \in \mathbb {R}^{5N\times 5N}\) and F(t)=1_N⊗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 A_i. \(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, H_2,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)=x¹(t)=…=x^N(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 A_i. Therefore, linearising Eq. (3) and block-diagonalising it with the linear transformation ξ=(P⊗I₅)⁻¹δx, where P is the matrix of eigenvectors that diagonalises G and I_n is the n-dimensional identity matrix, we obtain

Here Λ is diagonal holding the eigenvalues λ_i of G, i.e. GP=PΛ. Because both I_N and Λ are diagonal, Eq. (4) is block diagonal in \(\xi ^{\mu } \in \mathbb {R}^{5}\) via

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 A_i, i.e. \(\Delta =\sum _{i} \Delta _{i}\), and ξ(0) some initial perturbation. Then the perturbation after one period is given by ξ^μ(Δ)=Γ(σλ_μ)ξ^μ(0) with

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, σ⁻¹=τ_c, \(\tau _{\text {sr}}^{-1}=0\) and the only non-vanishing component of H is H_2,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 ²⁺ 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 ²⁺ amplitudes. When comparing Fig. 2b and c, we find that when the Ca ²⁺ 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 q_c leaves the ellipse with cos(argq_c)=−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 T_p=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 ²⁺ 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

Full size image

We now contrast the results for cytosolic coupling with those for purely luminal coupling. Here, σ⁻¹=τ_sr, \(\tau _{\mathrm {c}}^{-1}=0\) and the only non-zero component of H is H_3,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 T_p. As we increase T_p, 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 T_p=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

Full size image

In Fig. 3b, we plotted the stability regions as a function of the general MSF parameter η and the pacing period T_p. 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 T_p. 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 ²⁺ 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 T_p, 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 T_p. When we fix T_p 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 ²⁺ 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 ²⁺ 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 T_p 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)

Full size image

Peak values of the bulk Ca²⁺ concentration during one pacing period for σ_SR=3.2 (a), σ_SR=4.1 (b) in a 10×25 network of CRUs at T_p=1.05

Full size image

We can explain these sudden changes in the activity patterns of the intracellular Ca²⁺ 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 λ₁, while the region to the right is controlled by λ₃.

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 T_p and move from the ‘bump’ for λ₁ to the ‘bump’ for λ₂. 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 ²⁺ activity patterns in the network. For weak coupling, only one eigenvalue (λ₁) drives the instability. Upon increasing the coupling, both σλ₁ and σλ₂ 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 σλ₂ 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 T_p=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

Full size image

Figure 3b already suggests that the regions where the MSF is negative strongly depends on the pacing period T_p. 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 T_p in Fig. 7a. This three-dimensional plot highlights that the topology of the zero-contour changes significantly as a function of T_p. 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 T_p, the narrow part of the bowtie contracts, until two disconnected regions emerge as exemplified by Fig. 3a. When we lower T_p 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 T_p. Zero-contour of the MSF for T_p=1.075s (b) and T_p=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(η)))

Full size image

So far, we have studied purely cytosolic and luminal coupling, respectively. Under physiologically realistic conditions, however, Ca ²⁺ 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 (A_k−ηH), k=1,…,N, in Eq. (6) with (A−H(η_c,η_sr)). Here, H(η_c,η_sr) is a matrix whose non-zero elements are H_2,2=η_c and H_3,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 T_p=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)

Full size image

Networks are ubiquitous in biology, and intracellular signalling cascades constitute a prime example. In the present study, we investigated the Ca ²⁺ 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 ²⁺ 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 ²⁺ diffusion in the cytosol and Ca ²⁺ 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 ²⁺ network is associated with the emergence of cardiac Ca ²⁺ 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 ²⁺ alternans could only be observed at an advanced stage, i.e. when the Ca ²⁺ concentration oscillates out-of-phase in mesoscopic parts of a myocyte. However, improved imaging techniques now allow the recording of microscopic Ca ²⁺ alternans (Tian et al. 2012). Here, the Ca ²⁺ concentration averaged across the cell suggests a healthy cardiac myocyte, while in fact single CRUs may already display pathological Ca ²⁺ alternans. From a conceptual point of view, microscopic Ca ²⁺ alternans correspond to the pattern that emerges when the synchronous network state has just lost stability.

We recently reported that microscopic Ca²⁺ 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 ²⁺ 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 ²⁺ current in the present study is either zero or takes on a constant non-zero value i_CaL. If the nonlinear shape of the L-type Ca ²⁺ current is central to a study, it requires either the construction of a piecewise constant L-type Ca ²⁺ current with multiple levels, or one needs to resort to the nonlinear model. Similar considerations apply to the closure of the L-type Ca ²⁺ channel. We currently include voltage-dependent inactivation only. However, Ca ²⁺-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 ²⁺-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 ²⁺ buffers are the focus of attention. The current model treats buffer contributions as constant, but given the dynamic nature of the intracellular Ca ²⁺ concentration, the fraction of Ca ²⁺-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 Ca²⁺ 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 Ca²⁺ 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 ²⁺ diffusion is stronger in the cytosol or in the SR, different microscopic Ca ²⁺ alternans emerge, with each mode of diffusion giving rise to distinct network patterns of the intracellular Ca ²⁺ 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 Ca²⁺ current is given by I_CaL=Θ(V−V_L)i_CaL with the threshold voltage V_L=V_max−1 and a constant conductance

where a_CaL=V_maxF/RT. The current through the NCX is modelled as \(I_{\text {NaCa}} = \phi (V)-\psi (V)c_{\mathrm {s}}^{\mu }\), where

Here, we introduce the piecewise constant function

which switches at a threshold voltage of V_NaCa=−50 mV. The function η_q is given by η_q=0.0501α²+0.3816α+0.9182 with α=FV/RT. Ca²⁺ 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

where k measures the steep nonlinear dependence of Ca²⁺ release on the SR Ca²⁺ concentration, and s is a constant that is chosen such that Q is continuous. We model the clamped voltage for a pacing period T_p as

where \(k \in \mathbb N\) counts the number of APs and x=a_x/(a_x+T_p) with a_x=2/3. The resting potential is given by V_min=−70mV, and V₊(t) captures the shape of the clamped voltage, given by

for kT_p≤t≤(k+x)T_p, where the maximal AP is given by V_max=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

Not applicable.

Ca²⁺ :

Calcium

CRU:

Calcium release unit

MSF:

Master stability function

ODE:

Ordinary differential equation

PWL:

Piecewise linear

SR:

Sacroplasmic reticulum

Not applicable.

This work was supported by the Engineering and Physical Science Research Council (grant number EP/P007031/1).

Authors and Affiliations

1. School of Mathematical Sciences, University of Nottingham, Nottingham, NG7 2RD, UK

   Yi Ming Lai, Joshua Veasy, Stephen Coombes & Rüdiger Thul

Contributions

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.

Corresponding author

Correspondence to Yi Ming Lai.

Ethics approval and consent to participate

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Reprints and permissions