Introduction

Biological systems like the heart typically display complex dynamical behaviour through interactions among components such as genes, proteins, or cells (Dana et al. (2008)). These systems exhibit non-linear dynamics, feedback loops, and adaptability, reflecting phenomena like homeostasis, evolution, and rhythm stability (Xiong and Garfinkel (2023)). Central to the framework of dynamical systems is the concept of the phase space, a multidimensional representation in which each point corresponds to a unique state of the system. In geography when a small ball rolls on a raised-relief map, a state would include the planar position, elevation, speed and acceleration of the ball.

Within the phase space, the long-term behaviour of a system can often be described in terms of attractors, which are subsets of the phase space that trajectories tend toward over time. Attractors can take various forms, including fixed points, periodic orbits, and chaotic sets, reflecting the rich diversity of behaviours possible in dynamical systems. By replacing the ball by local rain in the raised-relief map, rivers, lakes and other bodies of water are attractors where that rain is collected. A watershed, also known as a drainage basin, is an area of land where any rain that falls will drain into the same river, lake or body of water. Therefore they are basins of attraction. Mountain ridges are separatrices that divide the watersheds, or more formally theoretical surfaces that divide the phase space into different basins of attraction. A separatrix therefore denotes the transition between different regions of the phase space. The terminology introduced here comes back throughout our paper, with the addition of a “state”, i.e. one of the attractors in the system. When the system is not interacted with, it will stay at the same point(s) in phase space.

Some systems only have a single basin of attraction. However, in others the phenomenon of bi-stability emerges, where two distinct attractors coexist, each with its own basin of attraction. Bistability is a hallmark of many natural and engineered systems, from climate dynamics to neural networks, ecological models, and cardiac dynamics (Angeli et al. (2004)). The interplay between basins in such systems can give rise to complex behaviours, making the characterization of these basins a critical component of a system’s analysis.

In this manuscript, the focus is on the complex system of the heart. At rest, ventricular cardiomyocytes (i.e muscle cells of the lower heart chamber) have a membrane potential of −80 to −90 mV. This resting state is an attractor of the system. When pushed out of equilibrium (e.g. due to an influx of Na⁺ ions or electrical stimulation), cardiomyocytes first depolarize and subsequently get attracted again towards the resting state by an efflux of K⁺ ions. In phase space the trajectory that gets mapped out resembles a cycle, corresponding to one heartbeat. When multiple of these cells are connected to each other, emergent behaviour can arise. Within this dynamic system of cardiomyocytes, we investigated emergent bi-stability in cell monolayers under the influence of spatial depolarization patterns. By using a train of external (electrical) waves, the system could be actively moved from one basin of attraction to another.

Cardiac arrhythmias are typically initiated by specific triggers, with point-source waves, i.e. ectopic waves being the most prevalent. To study ectopy in laboratory settings, control needs to be exerted over the depolarisation of cardiomyocytes. This can be achieved by optogenetics, i.e. a biological technique to control the activity of cells with light. This approach relies on the introduction of recombinant genes encoding light-activatable proteins into cells to endow them with a new biological function (Tan et al. (2022)). One way to depolarise a cardiomyocyte is by using mini-singlet oxygen generator (miniSOG). Upon activation with blue light, miniSOG converts molecular oxygen from its ground state (³O₂) into highly reactive singlet oxygen (¹O₂). The resulting oxidative damage induces a state of prolonged depolarisation in cardiomyocytes following electrical stimulation (Jangsangthong et al. (2016)). In such a system, it was previously observed that spatiotemporal illumination can give rise to collective behaviour and ectopic waves (Teplenin et al. (2018)) originating from illuminated/depolarised regions (with high curvature). These waves resulted from the interplay between the diffusion current and the bi-stable state that was induced in the illuminated region.

Although miniSOG can create depolarised areas, it is difficult to control the process precisely since it depolarises cardiomyocytes through an indirect pathway. Therefore, in this manuscript, we used light-sensitive ion channels to obtain more refined control over cardiomyocyte depolarisation. These ion channels allow the cells to respond to specific wavelengths of light, facilitating either depolarization or hyperpolarisation (Ördög et al. (2021, 2023)) through a direct pathway. By inducing cardiomyocyte depolarisation or hyperpolarisation only in the illuminated areas, optogenetics enables precise spatiotemporal control of cardiac excitability, an attribute we exploit in this manuscript (Appendix 2 Figure 1). With the aid of the light-gated ion channel CheRiff (Hochbaum et al. (2014)), ectopic waves similar to those obtained in the miniSOG system could be elicited (Appendix 2 Figure 2). Due to the superior control over depolarisation imposed by CheRiff, we uncovered new effects and unobserved behaviour related to the transition from a mono-stable cardiac system to a bi-stable one that displays periodic oscillations. This transition resulted from the interaction of an illuminated region with depolarized cardiomyocytes and an external wave train. More specifically, we found that our system displayed frequency selectivity: it oscillated only as a reaction to a specific range of inter-pulse external stimuli times. In neuroscience, this behaviour is classified as “resonance”. However, to avoid confusion with classical resonance, where the strength of the response is dependent on the input, here the phenomenon is called “induced pacemaker activity”. Evidence for this induced pacemaker activity will be presented both in vitro and in silico.

Figure 1.

Figure 2.

Results

Bi-stability in cardiomyocyte monolayers was realised both in vitro and in silico (Figure 1). In vitro this happened through patterned illumination applied to neonatal rat ventricular myocyte (NRVM) monolayers, which were optogenetically modified using the blue light-activatable cation channel CheRiff (Appendix 2 Figure 1). In silico, use was made of a slightly modified (see Methods and Materials) NRVM computational model (Majumder et al. (2016)) complemented with a model for a depolarising light-gated ion channel (Williams et al. (2013)). In both cases, the illumination pattern consisted of a rectangle of variable light intensities. After initiating a periodic wave train (of variable pulse number and interpulse interval/period) from a pacing electrode, two main outcomes and one special case were observed (Figure 1A). The first outcome was that the tissue stayed in the stationary rest state (STA). The second outcome was that the illuminated region became a pacemaker, which sent out pulses from its edges as a form of oscillatory behaviour (OSC). The special case consisted of an in-between case, in which the illuminated region only sends out a limited amount of pulses before returning to the stationary state thus showing transient oscillatory behaviour (TR. OSC). Except for the occasional mention, this last behaviour of transient oscillations will be reserved for more detailed discussion in the supplement. All light intensity-influenced monolayer behaviours are visualised together in Appendix 2 Figure 3. Here, we will focus on the two main behaviours after periodic wave trains: STA and OSC. As shown below, these two behaviours can exist within the same monolayer under bi-stability and transitions from one state to the other are possible using periodic wave trains. The occurrence of this phenomenon appears to be influenced by the light intensity applied to the illuminated region (Figure 1B). In a step-by-step manner, we unravelled the intricacies of this tissue level (i.e. collective) bi-stability between oscillatory and stationary behaviour, linking experiment to theory.

Figure 3.

In vitro realisation of collective light-induced bi-stability

To test the scope of collective bi-stability in vitro, multiple parameters were varied within our experimental set-up. The first one of these was the frequency of stimulation, visualised in Figure 2, where the number of electrical pulses and the light intensity (in the range 0.03125-0.25 mW ∕mm²) were kept fixed. In every NRVM monolayer (26 biological replicates/monolayers), three frequencies (stimulation periods 200, 600, and 1200 ms) were tested subsequently, with a total of three repetitions (technical replicates) of the full set, totalling 26 × 3 = 78 observations at each pacing frequency. Optical voltage signals were measured both at the center of the illuminated area (top traces) as well as in the bulk of the tissue (bottom traces), in order to visualise action potentials throughout the tissue. Underneath these optical voltage traces is a space-time plot of the experimental monolayer, the vertical axis showing space and the horizontal one showing time. At each moment in time, the middle line of the monolayer image was plotted. White shows depolarized tissue, while black signifies repolarised cells. From these plots, intermediate stimulation periods (600 ms) were observed to change the monolayer from a collective stationary state to a collective oscillatory state, while low (200 ms) or high (1200 ms) ones had no effect. It was observed that the oscillatory state consisted of ectopic pacemaker activity with a frequency (735 ms-period) that was different from the stimulus train frequency (600-ms period). Similar results were found for the transient oscillatory state (Appendix 2 Figure 4A). In short, we found that collective bi-stability becomes apparent in the mid-range of electrical pacing frequencies.

Figure 4.

A second parameter tested was the number of pulses, visualised in Figure 3, where the period of stimulation (600 ms) and the light intensity (in the range 0.03125-0.25 mW ∕mm²) were kept fixed. In every NRVM monolayer (7 biological replicates/monolayers), two different pulse numbers (2 and 4) were tested subsequently, with a total of three repetitions (technical replicates) of the full set, totalling 7 × 3 = 21 observations for each number of pulses. Once again optical voltage traces were taken from the center of the illuminated area and from the bulk of the tissue. Our observations show that a minimum number of pulses is needed to change the collective behaviour of the monolayer from stationary to oscillatory.

Taken together, we found a way to control transitions of collective behaviour of cardiomyocyte monolayers from stationary to oscillatory through two different parameters, i.e. the frequency and number of pulses.

In silico realisation of collective light-induced bi-stability

Using the power of state-of-the-art computational modelling, we performed similar experiments to the in vitro ones with a modified version of the Majumder et al. (2016) model for NRVMs combined with the Williams et al. (2013) model for the H134R mutant of Chlamydomonas reinhardtii channelrhodopsin-2. As this light-sensitive ion channel slightly differs from CheRiff in its biophysical properties, the in vitro findings will have a qualitative rather than a quantitative correspondence to the in silico results, which are visualised in Figure 4. Under fixed constant illumination (1.7124 mW ∕mm²) and a fixed number of 4 pulses, frequency dependency of collective bi-stability was reproduced in Figure 4A. Once again, only intermediate electrical pacing frequencies (500-ms period) caused transitions from collective stationary behaviour to collective oscillatory behaviour and ectopic pacemaker activity had periods (710 ms) that were different from the stimulation train period (500 ms). Figure 4B shows the accumulation of pulses necessary to invoke a transition from the collective stationary state to the collective oscillatory state at a fixed stimulation period (600 ms). Once more, ectopic pacemaker activity had periods (750 ms) that were different from the stimulation train period (600 ms). Also for the transient oscillatory state, the simulations show frequency selectivity (Appendix 2 Figure 4B).

Having established a qualitative correspondence between the in silico model and in vitro experiments, the unique properties of the computational model were exploited to assert control over all variables in order to determine the parameter ranges of the state transitions. We previously identified three important parameters influencing such transitions: illumination strength, number of pulses, and frequency of pulses. Scanning through all three variables revealed the specific conditions at which transition from the stationary state towards the oscillatory state occurred as shown in the rightmost panel of Figure 4C with slices for fixed illumination extracted from this 3D-plot shown at the left. With increasing light intensity, the parameter range (pulse and period) in which bi-stability can be detected increases (the total area of TR.OSC. + OSC. becomes larger), and transient oscillatory states transition towards oscillatory states.

Except for the three aforementioned parameters, a fourth one was found to have an influence on the state transitions: the size of the illuminated region. Also this parameter was first tested in vitro resulting in 5 biological replicates/monolayers showing transitions to the oscillatory state and 6 biological replicates/monolayers showing transient oscillations. The effect of the size of the illuminated area on the occurrence of bi-stability was qualitatively reproduced in silico (Appendix 2 Figure 5). Varying the number and frequency of pulses together with the illuminated area at a fixed light intensity of 1.7200 mW ∕mm² resulted in the dynamic behaviour shown in Figure 4D. For the conditions that resulted in stable oscillations, the green vertical lines in the middle and right slices represent the natural pacemaker frequency in the oscillatory state. After the transition from the stationary towards the oscillatory state, the illuminated region settles into this period of sending out pulses. The location of this line depends on the light intensity and the size of the illuminated region. However, it is independent of the number of pulses and the frequency of the stimulus train. Figures Figure 4C and Figure 4D also illustrate that no oscillatory state could be induced with a high number of pulses for stimulus trains of high frequency (i.e. low stimulation period).

Figure 5.

Termination of collective pacemaker activity

A closer look at the high frequency range showed that it not only prevented the initiation of ectopic pacemaker activity, but can also cause its termination (Figure 5A). In NRVM monolayers (3 biological replicates/monolayers), collective ectopic pacemaker activity was first initiated using a stimulus train of 3 pulses with a period of 500 ms. Once ongoing, 8 pulses in the high frequency range (period 200 ms) terminated the ectopic pacemaker activity and made the tissue go back to its resting state. When copying this experimental protocol in silico, qualitative correspondence was shown using 4 pulses with a period of 500 ms to initiate pacemaker activity, and 7 pulses with a period of 180 ms to terminate it. Note that both experimentally and computationally every pulse of the stimulus train is captured (orange traces), despite the fact that the voltage traces from points inside the illuminated region do not give that impression.

Once again, the qualitative in silico correspondence allows us to determine the specific conditions at which initiation and termination of collective ectopic pacemaker activity occur in more detail (Figure 5B). In the period vs. number of pulses plot, which is a copy of the rightmost slice in Figure 4C, the region showing termination (red-lined) did not overlap with the region showing initiation (magenta-lined).

The collective nature of wave train-induced pacemaker activity

After exploring collective bi-stability in practice, we aimed to gain a deeper understanding of this phenomenon. Bi-stability is a wide-spread phenomenon and can take place at the single cell level (Guevara (2003)), which would be the simplest assumption for our monolayer system as well. When this would be the case, all cells would become oscillatory under illumination. We tested this hypothesis by making use of one of the four parameters we investigated earlier: the illuminated area. This parameter shows a transition from the stationary state towards the oscillatory state (Figure 6A) just like the light intensity parameter (Figure 1B) resulting in spontaneous oscillatory behaviour with waves emanating from the center of the illuminated area for a large area size and at a constant low light intensity (Figure 4D). These circumstances allowed us to compare the voltage traces between the center of the irradiated area and the bulk of the tissue, both in vitro and in silico. While the cells outside of the illuminated area showed oscillating membrane potentials those in the center of the illuminated region, i.e. those least affected by the coupling current of neighbouring cells, displayed a stable elevation of membrane potential. Hence, this implies that not all illuminated cells became oscillatory. We can therefore reject the simplest hypothesis of bi-stability at the single cell level, but have to speak of a collective phenomenon. It is for this reason that throughout the manuscript the observed effects are referred to as “collective” pacemaker activity.

Figure 6.

Insight into pulse-dependent collective state transitions

Knowing that we are dealing with a collective phenomenon, we wanted to gain a deeper understanding in the observed pulse-dependent collective transitions from stationary towards oscillatory pacemaker activity. In order to achieve this goal, in silico analysis was performed (Figure 7) because it offers the possibility to visualise all parameters relevant for the system. In every monolayer, we selected a single point from the centre of the illuminated region of the tissue and looked at its dynamic behaviour in phase space.

Figure 7.

As was shown in all figures from Figure 2 to Figure 6, illuminated cells do not show dampening oscillations when reaching the stationary state after a transition from the oscillatory state. This behaviour is justified when looking at a projection of the phase space of the model where the voltage V and the gating variable x_r for the rapid delayed outward rectifier K⁺ current were chosen for the x- and y-axis, respectively (Figure 7A). The stable fixed point (STA) is not located inside the stable limit cycle (OSC) (right panel), hence showing no Hopf bifurcation (left panel) with its associated damped oscillations. The two attractor regions (STA and OSC) were also reconstructed from experimental data and showed two distinct trajectories (Appendix 2 Figure 6).

Since we are dealing with a high-dimensional phase space in the in silico model, the visualisation of trajectories in time and the corresponding dynamic behaviour is rather complex. However, from the two-dimensional projections onto the V and x_r variables, we can see what causes the frequency and pulse dependency (Figure 2, Figure 3 and Figure 4) of ectopic pacemaker activity initiation, as visualised in Figure 7C. Four different pairs of stimulation period and number of pulses were selected, corresponding to different behaviour areas in the parameter space slice for a fixed illumination intensity (1.72 mW ∕mm²) and size of the illuminated area (67 pixels edge length) (Figure 7B). The first point (⍺) comes from the region that shows no transition towards collective pace-maker activity for a low number of pulses (2) and intermediate stimulation periods of 800 ms. The phase space projection (Figure 7C) shows that the pulses cycle around the stable limit cycle, and return back to the stable fixed point, indicating that two pulses are not enough to move within the attractor zone of the stable limit cycle. When the number of pulses is increased to 4 (β), the attractor zone of the stable limit cycle is reached from the outside, and collective ectopic pacemaker activity is induced. After switching to a faster stimulation period of 250 ms and sticking to just 3 pulses (γ), the phase space trajectories end up inside the stable limit cycle, but still within the attractor zone. Hence, also this time collective pacemaker activity was induced. If we now increase the number of pulses again, this time to 8 (δ), we are inside the stable limit cycle, but move again out of the region of attraction. Instead of resulting in oscillations, the system moves back to the stationary state through other dimensions by the electrotonic current.

Thus, there is frequency selectivity due to the region of attraction of the stable limit cycle, where too low and too high frequencies result in stationary behaviour. Also pulse accumulation is shown, where more pulses move the system closer to the attractor region. However, we are not dealing with a Hopf bifurcation, because otherwise we would not get back to the stationary state when the number of pulses is increased too much.

Insight into the termination of ectopic pacemaker activity

Let us look once more at a cell inside the illuminated tissue undergoing a stimulus train with a short period of 180 ms from outside and extract all time-dependent variables. In response to a high number of stimuli (17 pulses), the time trace in the phase-space projection plot goes inside the limit cycle (Figure 8A) towards a quasi-steady-state, and returns back to the stationary state through a strangely curved line (black). This quasi-steady-state is hidden within the oscillatory trajectory and no pacemaker activity is initiated. When the number of pulses is increased from 17 to 37, the voltage trace shows exactly the same behaviour, indicating that the number of pulses has no influence on the mechanism that moves the tissue back to the stationary state (Figure 8B). However, repolarisation reserve does have an influence and makes the transition longer (Appendix 2 Figure 7). This can be observed by either moving closer to the boundary of the illuminated region where there is a stronger diffusion/electrotonic influence from the non-illuminated region, or by introducing ionic changes.

Figure 8.

To understand what is going on, we studied single cell dynamics under illumination. We already know from Figure 6 that we are dealing with a purely collective phenomenon. Therefore, the single cell dynamics itself showed nothing peculiar, i.e. a standard current-voltage relationship (IV) curve with a single stable equilibrium (not shown). To add a multicellular component to our single cell model we introduced a current that replicates the effect of cell coupling and its associated diffusion. This so-called bias current equals the average in time and space of the current inside the illuminated region using wave trains with a period of 180 ms (the current produced during the quasi-steady-state in Figure 8B). When this bias current is present, the IV curve of the single cell displays two stable equilibria and shows bi-stability (Figure 8C-1). The second stable equilibrium cannot be reached with a single pulse, but needs a stimulus train and can only be reached with short stimulation periods (180 ms) and intermediate stimulation amplitudes (−0 mV) (Figure 8C-2).

Building on this result, we next applied global stimuli with a period of 180 ms and amplitudes 58, −20, and 50 mV to a virtual cardiomyocyte monolayer with an illuminated square in the center (Figure 8D). These pulses were applied when the tissue was either in the oscillatory state (top two rows) or the stationary state (bottom two rows). The stimulus amplitude of −20 mV that was necessary to get the single cell with a bias current into the second stable equilibrium, also terminates oscillatory behaviour (Figure 8D-2) while the other two amplitudes of −58 and 50 mV (Figure 8D-1,3) don’t alter the behaviour of the tissue. However, when starting from the stationary state, the behaviour is opposite and the intermediate amplitude of −20 mV fails to induce collective oscillatory behaviour while the other two amplitudes do elicit this behaviour. We thus conclude that high stimulation frequencies and intermediate stimulus amplitudes are capable of terminating collective oscillatory behaviour by transitioning through a quasi-steady-state.

We call this kind of bi-stability capable of terminating collective oscillatory behaviour “hidden bistability” for two reasons: (1) Single cell bi-stability depends on electronic interactions. Hence it can be said that this bi-stability is “hidden” and only emerges in collective behaviour. (2) It’s possible to escape the region of attraction of the oscillatory trajectory through opposite paths.

It’s complex, but not as complex as it seems

While it may be believed that the observed phenomena and explanations are intrinsic to entangled interactions specific to the complexity of a biological sample (in vitro) and a realistic computational model (in silico), we show here that this is not the case. By reducing the model equations, we show that collective oscillatory behaviour is not limited to the specific case of NRVMs. We started from the initial NRVM equations and reduced them to contain only a single depolarising current (I_Na), two repolarising currents (I_Kr and the time-independent I_K1), and a spatially dependent optogenetic current (I_ChR2). Fast variables in these currents were assumed to keep their steady-state value, which caused them to be time-independent. The exact equations can be found in Appendix 1.

Using this simplified system of three variables (the voltage V, the closing gate of the Na⁺ channel h, and the opening gate of the K⁺ channel x_r), we were able to reproduce all key features (frequency selectivity, pulse accumulation and bi-stability) of the more complex systems (Figure 9).

Figure 9.

Just like in Figure 7C, the simplified model showed no initiation of collective pacemaker activity at a low number of 2 pulses and a low stimulation frequency of 2250 ms or at a high number of 8 pulses and a high stimulation frequency of 540 ms (Figure 9A). However, pacemaker activity was induced when applying 3 pulses and a low stimulation frequency or 2 pulses and a high stimulation frequency. At a low frequency, the phase space projection (again projected onto the V and x_r variable plane) shows that the limit cycle is approached from the outside, while at high frequencies it is approached from the inside. Altogether, these results show frequency selectivity and pulse accumulation in the simplified model (Figure 9A). To explain the findings in the rightmost panel of Figure 9A, we had a closer look at the single cell dynamics under high stimulation frequencies. As shown in Figure 9B, 4 pulses at a low frequency (1000-ms period) resulted in normal behaviour, but high frequency stimulation (500-ms period) pushed the cell towards a stably elevated membrane potential (Figure 9B). Looking at the influence of stimulus amplitude on the ability to get to this elevated membrane potential, it appears that after 7 pulses 3 of the 4 traces move towards the elevated membrane potential of −20 mV (Figure 9C). Panels B and C combined give us the same qualitative behaviour as seen in Figure 8C. The qualitative behaviours of Figure 4A and Figure 8D were also reproduced and can be found in Appendix 2 Figure 8.

Discussion

In this study, in vitro observations were combined with in silico modelling to study and gain mechanistic insight into collective electrophysiological behaviour in NRVM monolayers.

• Firstly, we demonstrated the existence of transitions from stationary cardiac monolayers towards ones that show pacemaker activity coming from an illuminated region, which are induced by wave trains of external pulses (Figure 1) as a manifestation of collective bi-stability. This transition between the two states is dependent on four variables: the light intensity applied to the illuminated region, the size of the illuminated region, the number of pulses of an incoming wave-train, and the frequency of the pulses of the incoming wave train (Figure 2, Figure 3, Figure 5 and Appendix 2 Figure 4, Figure 5). In vitro experiments were quantitatively reproduced in silico, thereby establishing a link between model and experiment (Figure 4, Figure 5 and Appendix 2 Figure 4, Figure 5)).

• Secondly, we studied the mechanisms underlying the initiation and termination of pacemaker activity. In the presence of a diffusion current (or bias current in the case of a single cell), the illuminated cells transition to a second equilibrium under high frequency pacing. Because this second equilibrium only appears in a coupled system (or in the presence of a bias current), we can also speak of a hidden bi-stability. Hence, hidden bi-stability (high vs. low stationary resting membrane) is shown to underlie the collective phenomena (stationary vs. oscillatory behaviour) (Figure 6 Figure 7 Figure 8, Appendix 2 Figure 6 and Appendix 2 Figure 7).

Below, we compare these phenomena to different types of arrhythmic events and mechanisms of bi-stability in cardiac tissue, and extend the relevance of our mechanisms to spiking behaviour in neuronal and other biophysical systems.

Ectopic pacemaker activity in cardiology

The generation of ectopic beats is considered as one of the main mechanisms triggering and maintaining cardiac arrhythmias (Pogwizd et al. (1992)). In most cases, ectopic sources are believed to be a result of triggered activity, which can manifest itself as early afterdepolarizations (EADs) or delayed afterdepolarizations (DADs) (Qu et al. (2014)). Such ectopic sources are typically initiated by high frequency pacing resulting in a Ca²⁺ overload that triggers EADs and DADs (Lerman et al. (2024)). However, as we didn’t observe Ca²⁺overload, the phenomenon observed in this study resembles induced pacemaker activity instead. In contrast with ectopy, this pacemaker activity is initiated at low pacing frequencies. Therefore our findings might explain the phenomenon of bradycardia-induced triggered activity (Brachmann et al. (1983)). The initiation of oscillations in our monolayers is the result of the spatial dynamics of the system, which stands in contrast with the frequency dependency of EADs that can be explained on a single cell level (Tran et al. (2009); Qu et al. (2014)). Our optogenetic current is absent in normal cardiac tissue, but its associated depolarisation could be caused by other late activating inward currents, which are also essential in the induction of EADs or DADs (Song et al. (2008); Qu et al. (2014)): e.g. the late Na current (I_Na,L, Belardinelli et al. (2015)), reactivation of I_CaL (Tran et al. (2009)), or NCX activity (Song et al. (2015)). In addition, collective bi-stability might be observed in the presence of oxidative stress, ischemia or chronic heart failure (Belardinelli et al. (2015); Shryock et al. (2013); Undrovinas et al. (1999)), where there is depolarisation-induced automaticity (triggered activity) caused by a background current (e.g. I_Na,L).

The ectopic pacemaker activity in this study is a direct result of the combination of single cell depolarisation with spatial effects. Due to the additional inward current produced by the activation of CheRiff, a single uncoupled illuminated cell goes to a steady depolarised state (Figure 6B). However, in a multicellular environment there is a repolarising electrotonic coupling current working against this depolarising light-induced current. This repolarising current is a consequence of the interaction between depolarised and non-depolarised parts in the monolayer. By playing with the illumination intensity and the size of the illuminated region, we were able to balance these de- and re-polarising effects to create a bi-stable regime (STA vs. OSC). This also explains why there are no oscillations present in the center of a large illuminated region, since there the electrotonic current is lower than at its edges. (Figure 6C). On the other hand, if the electrotonic current is too strong, it will prevent depolarisation by the optogenetic current and, as a consequence, the system wil stay in the stationary state. The balancing of currents also explains why the oscillatory boundary moves closer to the edge of the illuminated region at higher light intensities, since under those circumstances the depolarising current will overpower the spatial effects. However, when the illumination is very strong, it is possible to have oscillations that don’t need to be induced (Figure 1B) and can even come from the corners of the illuminated region (Appendix 2 Figure 3, Teplenin et al. (2018)). This last-mentioned phenomenon appears because the curvature of the boundary between illuminated and not-illuminated regions influences the electrotonic current.

Interplay between microscopic (single cell) and macroscopic (monolayer) behaviour

Collective bi-stability, as described in the previous section, is a bi-stability between two distinct behaviours (stationary STA and oscillatory OSC behaviour) that occur at the macroscopic level (mono-layer). However, by looking more deeply at the dynamics that happen at the microscopic level (single cell), we can get more insight into the underlying mechanism(s). At the single cell level, the dynamic transition between a stationary monolayer and one showing induced pacemaker activity relies on two main processes:

1. For the initiation of pacemaker activity, pulse accumulation (at intermediate frequencies) is necessary to bring the system into the attracting region of the oscillatory state.

2. For the termination of pacemaker activity, an interaction takes place between “hidden” bistable dynamics and coupling currents at high frequencies.

Both these processes occur as a consequence of wave pulses outside of the illuminated region. At low frequencies, no oscillatory activity can be induced. At intermediate frequencies, it is possible to escape the attracting zone of the stationary state and cross the separatrix into the attracting zone of the oscillatory state. However, for even higher frequencies, the system gets pushed out of the attracting zone of the oscillatory state again and through hidden bi-stability it returns to the stationary state.

Our data on the termination mechanism of induced pacemaker activity (Figure 8) indicate that we are not dealing with classical STA-OSC bi-stability (Winfree (2001)). In systems displaying classical bi-stability between a limit cycle and the resting state, pacemaker annihilation can be achieved by applying sufficiently high frequency voltage perturbations. Through pulse accumulation it is then possible to exit the basin of attraction in the reverse direction from which it was entered. In contrast, our system enters (from outside the limit cycle) and leaves (towards the inside of the limit cycle) the attraction zone in opposite directions. This was concluded from the absence of a Hopf bifurcation in our system (Figure 7).

The seemingly fragile nature of hidden bi-stability might be a requirement for coexistence of excitable (non-illuminated) and bi-stable (illuminated) regimes. Otherwise, if the hills and dips of the N-shaped steady-state IV curve are large (Figure 8C-1), they would be in the same range as the big currents of fast ion channels, thereby preventing the subtle interaction between the large ionic cell currents and small repolarising coupling currents.

Although the termination mechanism was shown to be amplitude-dependent (Figure 8C-D), the appropriate perturbation amplitudes are “automatically” selected by a high frequency train of travelling waves of excitation without any a priori knowledge about the state of the system. This wave train terminated ectopic triggers (Figure 2), which is in accordance with previous publications (Moak and Rosen (1984); Glikson et al. (2021)). These publications also showed that uniform high voltage shocks are unable to terminate ectopic sustained activity, similar to our findings (Figure 8). Although its manifestations are clear in a monolayer system, hidden bi-stability might be challenging to be observed in single cell patch-clamp experiments. In Figure 8C-1, the height of dips and hills of the N-shaped steady-state IV curve under bias current are ≈ 0.2pA∕pF. While it is possible to perform patch-clamping on individual cells in a monolayer experiencing diffusion currents, measurements of small amplitudes can be hindered by the statistical and instrumental errors of voltage-clamp recordings. The voltage bi-stability we are looking for can, for example, easily be lost via a shift in current balance due to the non-zero resistance of the patch pipette in current clamp mode. This might explain why the phenomenon of hidden bi-stability has only rarely been reported, for example in an in silico study dedicated to myotony in muscle fibres (Cannon et al. (1993)).

Extrapolation and translation

All key aspects of the initiation and termination of induced pacemaker activity, i.e. the points itemized above, were also demonstrated in a simplified 3 variable reaction-diffusion model (Figure 9, Appendix 2 Figure 8), implicating the more general nature of this observed phenomenon. A rigorous mathematical analysis can be accomplished in future work, but might be challenging due to the lack of efficient analytical tools to describe fully developed limit cycle oscillations in reaction-diffusion systems (Sherratt and Smith (2008)). Based on the reproduction of our results in such a simplified model, we showed that we don’t need the full complexity of an ionic cardiomyocyte model, enabling extrapolation and translation of our results towards other scientific disciplines.

Frequency-dependent initiation and termination of oscillatory activity is not exclusive to cardiac systems, but is well known in a variety of non-linear biophysical systems. For example, the actomyosin cortex in the social amoeba Dictyostelium discoideum responds to short pulses of 3’,5’-cyclic adenosine monophosphate (cAMP) by damped oscillations in average actin content, thereby allowing fine time-selective responses to either internal (pseudopodia initiation) or external (cAMP release) perturbations (Westendorf et al. (2013)). Moreover, in synthetic gene regulatory networks, subthreshold stimuli might be used to discriminate environmental stimuli (e.g. growth factor release, heat shock) (Guantes and Poyatos (2006)). Also neurons can differentiate external stimuli depending on their frequency. As mentioned in the introduction, the frequency-dependent initiation and termination of oscillatory activity is known as ‘resonance’ in neuroscience. It is interesting to compare the resonance properties of our system with those actively studied in neurons.

Resonant neurons, also known as resonate-and-fire neurons or type II excitability neurons (Izhikevich (2000)), are governed by sub-threshold oscillations initiated due to proximity of the resting state to a Hopf bifurcation (Hutcheon and Yarom (2000); Roach et al. (2018)). They react to a stimulus train as follows: The first pulse initiates a damped sub-threshold oscillation of the membrane potential, while the effect of the second pulse depends on its timing. If timed correctly, it adds to the first pulse and either increases the amplitude of the sub-threshold oscillation or evokes repetitive spiking (transition to a limit cycle). Thus, the response and spiking of the resonator neuron depends on the frequency content of the input. These types of spiking neurons are usually associated with Bogdanov-Takens bifurcations (Izhikevich (2000)). However, another type of bifurcation, a saddle-node-loop (SNL) bifurcation, is also capable of generating spiking neurons and allows a bi-stability regime, in which the stable state lies outside the limit cycle similar to our case (Figure 7).

SNL bifurcations were previously identified in neurons (Hesse et al. (2017); Hesse and Schreiber (2019); Schleimer et al. (2021)), in a model of transcriptional regulation of the cell (Ciliberto et al. (2005)), and in the design of combined genetic switches (Perez-Carrasco et al. (2018)). Our result might represent the emergent analogue of a SNL bifurcation with resonant transitions on a neuronal/non-linear network level. SNL bifurcations have already been identified in the networks of coupled Kuramoto limit cycle oscillators both in a broad (Martens et al. (2009); Childs and Strogatz (2008)) and in a neuroscience context (Buendía et al. (2021)). In these studies, SNL bifurcations were identified using dimensionality reduction techniques while ignoring the whole wealth of transient and finite-size phenomena. Moreover the Kuramoto model is a model without many physical details, describing a weak forcing regime on an idealized limit cycle oscillator, in which oscillations cannot be destroyed but only phase shifted. In our case, we provide a more realistic biophysically motivated model, which allows a strong forcing regime and destruction of oscillations due to the phenomenon of hidden bi-stability. This gives us reason to believe that idealized Kuramoto models could possess resonant properties if augmented/modified with sufficient physical/biophysical details. SNL bifurcations and hidden bi-stability could also be studied in a range of other unexplored biophysical situations including high frequency stimulation for the termination of neuronal oscillations during Parkinson disease (Touboul et al. (2020)), multi compartment cable computational models of the neuron (Schwemmer and Lewis (2012)), and bi-stability in neurons in vitro (Le et al. (2006)). Additional studies can be carried out in neuroscience by local optogenetic activation of depolarizing ion channels, both in vivo Lu et al. (2015) and in silico in non-local models Selvaraj et al. (2015); Heitmann et al. (2017).

Conclusion

Using both in vitro and in silico models, we explored a mechanism for frequency-dependent initiation and termination of oscillatory activity (i.e. induced pacemaker activity) in non-homogeneous cardiac tissue. This mechanism could be extended to more general excitable systems and is in neuroscience also known under the term resonance. Our collective “resonance” mechanism could be used to study or control collective states in cardiac but also other excitable systems.

Materials and methods

Preparation of monolayers of NRVMs expressing the light-activatable cation channel CheRiff

All animal experiments were reviewed and approved by the Animal Experiments Committee of the Leiden University Medical Center and carried out in accordance with Directive 2010/63/EU of the European Union on the protection of animals used for scientific purposes.. Monolayers of NRVMs were established as previously described Engels et al. (2015). Briefly, the hearts were excised from anaesthetised 2-day-old Wistar rats. The ventricles were cut into small pieces and dissociated in a solution containing 450 U∕ml collagenase type I (Worthington, Lakewood, NJ) and 18,75 Kunitz/ml DNase I (Sigma-Aldrich, St. Louis, MO). The resulting cell suspension was enriched for cardiomyocytes by preplating for 120 minutes in a humidified incubator at 37^◦C and 5% CO₂ using Primaria culture dishes (Becton Dickinson, Breda, the Netherlands). Finally, the cells were seeded on round glass coverslips (d = 15 mm; Gerhard Menzel, Braunschweig, Germany) coated with bovine fibronectin (100 μg/ml; Sigma-Aldrich) to establish confluent monolayers. After incubation overnight in an atmosphere of humidified 95% air-5% CO₂ at 37^◦C, these monolayers were treated with Mitomycin-C (10 μg/ml; Sigma-Aldrich) for 2 hours to minimise proliferation of the non-cardiomyocytes. At day 4 of culture, the NRVM monolayers were incubated for 20-24 h with lentiviral vector particles encoding CheRiff (Hochbaum et al. (2014)), a light-gated depolarising ion channel of the channelrhodopsin family (Wang et al. (2015)), at a dose resulting in homogeneous transduction of nearly 100% of the cells (Feola et al. (2017); Ördög et al. (2023)). Next, the cultures were washed once with phosphate-buffer saline, given fresh culture medium and kept under culture conditions for 3-4 additional days.

Optical voltage mapping and patterned illumination of NRVM monolayers

After 8-10 days of culturing, NRVM monolayers were optically mapped using the voltage-sensitive dye di-4-ANBDQBS (52.5 μM final concentration; ITK diagnostics, Uithoorn, the Netherlands) as reported previously (Feola et al. (2017)). The mapping setup (Appendix 2 Figure 1) was based on a 100 × 100 pixel CMOS Ultima-L camera (Scimedia, Costa Mesa, CA). The field of view was 16 × 16 mm resulting in a spatial resolution of 160 μm/pixel. For targeted illumination of monolayers the setup was optically conjugated to a digitally controlled micro-mirror device (DMD), the Polygon 400 (Mightex Systems, Toronto, ON), with a high power blue (470 nm) LED (BLS-LCS-0470-50-22-H, Mightex Systems). Before starting the actual experiments, all monolayers were mapped during 1-Hz electrical point stimulation to check baseline conditions. Electrical stimulation was performed by applying 10-ms-long rectangular electrical pulses with an amplitude of 8 V to a bipolar platinum electrode with a spacing of 1.5 mm between anode and cathode. Only cultures with an action potential duration (APD) at 80% repolarization (APD₈₀) below 350 ms and a conduction velocity (CV) above 18 cm∕s were used for further experiments. Stimuli were applied in trains that varied in number and period. Whenever the outcome was reproducible three consecutive times in a single monolayer, it was counted as successful and included as one measurement in the total dataset. The constant light intensity was varied in the range of (0.03125-0.25 mW ∕mm²) to search the critical value for the bi-stable regime. For experiments in which the light intensity was varied the rectangular-shaped area of illumination was in the range of 4 × 4 to 6 × 6 mm. The highest achievable irradiation intensity (0.3125 mW ∕mm²) was used to perform size modulation experiments. The resulting electrical activity was recorded for 6-24 s at exposure times of 6 ms per frame.

Post-processing of optical mapping data

Data processing was performed using specialized BV Ana software (Scimedia), ImageJ (Schneider et al. (2012)) and custom-written scripts in Wolfram Mathematica (Wolfram Research, Hanborough, Oxfordshire, United Kingdom). APD and CV were calculated as described previously (Feola et al. (2017)). To prepare representative frames of wave propagation, optical mapping videos were filtered with a spatial averaging filter (3 × 3 stencil) and a derivative filter.

Attractor reconstruction

Multiple attractor reconstruction is based on Takens’s embedding theorem (Takens (1981)) by plotting phase space trajectories in (V (t), V (t+τ)) coordinates. The optimal time delay τ was determined by calculating the minimum of the autocorrelation function (Clemson and Stefanovska (2014)). Linear drift removal was applied before plotting phase space trajectories.

Detailed electrophysiological model

In this study, the Majumder-Korhonen model (Majumder et al. (2016)) of NRVM monolayers was used, whose monodomain reaction-diffusion equation is formulated as follows:

where I_ion is the sum of all ionic currents, 𝒟 is the diffusion coefficient equal to 0.00095 cm²/ms and I_ChR2(x, y) is the spatially controlled optogenetic current. In comparison to the original model, I_Na was increased 1.3-fold and, I_Ca and I_Kr were both reduced 0.7-fold. The steady-state voltage dependency for the inactivation variable h of the fast Na⁺ current was changed to , where V is the transmembrane potential. The resulting CV was 20 cm/s. The optogenetic current I_ChR2 was formulated based on the Boyle et al. (2013) I_ChR2 model with the introduction of rectifying conduction properties and voltage-dependent state transitions (Williams et al. (2013)). The irradiation intensity varied in the range of 0.165–0.175 mW/mm².

The forward Euler method was used to integrate the equations with a time step Δt = 0.02ms and a centered finite-differencing scheme to discretise the Laplacian with a space step of Δx = 0.0625mm. The total computational domain size was 256 × 256 grid points; while the centrally irradiated square with the optogenetic current consisted of 80 × 80 grid points. To create stationary initial conditions, the model was integrated for 2 minutes before a train of stimuli was applied from the upper border of the domain.

Acknowledgements

We would like to thank Desmond Kabus for useful discussions that improved the manuscript and Annemarie Kip and Cindy Bart for the production of the CheRiff-encoding lentiviral vector LV.GgTnnt2.CheRiff eYFP.WHVoPRE.

This work was supported by the European Research Council (ERC consolidator grant 101044831 to DAP) and the Netherlands Organization for Scientific Research (NWO Vidi grant 917143 to DAP and ZonMW Off Road grant 04510012110051 to TDC).

Additional information

Author contributions

Conceptualization: AST, AVP, DAP, TDC

Methodology: AST

Software: AST, RM, NNK

Formal analysis: AST, TDC

Investigation: AST Resources: AAFdV, DAP

Writing—original draft: TDC

Writing—review & editing: AST, AAFdV, AVP, DAP, TDC

Visualization: AST, TDC

Supervision: DAP and TDC

Project administration: DAP

Funding acquisition: DAP, TDC

Appendix 1 Supplementary Methods sample

Simplified reaction-diffusion model

The simplified reaction-diffusion model consists of one voltage variable, two gating variables for the depolarising and hyperpolarising currents, and an inhomogeneity parameter C(x, y)

where D is the diffusion coefficient equal to 0.0003 cm²∕ms. The I_K1 current formulation was copied from Majumder et al. (2016). The C(x, y) was set to 0.1518115 inside the illuminated square and to 0 outside this. are described as:

The forward Euler method was used to integrate the equations with a time step Δt = 0.01ms and a centered finite-differencing scheme to discretise the Laplacian with a space step of Δx = 0.0625mm. The total computational domain size was 256 × 256 grid points while the square with the “background” current consisted of 80 × 80 grid points. To create stationary initial conditions, the model was integrated for 2 minutes before a train of stimuli was applied from the upper border of the domain.

Appendix 2 Supplementary Figures sample

Figure 1.

Figure 2.

Figure 3.

Figure 4.

Figure 5.

Figure 6.

Figure 7.

Figure 8.