This article was previously published in Scientific Reports 14, 28962 (2024). (Arno et al., 2024a) and is a chapter of my dissertation.

Authors:

• Louise Arno^1,3 (ORCID: 0000-0001-9384-519X)

• Desmond Kabus^1,2,3 (ORCID: 0000-0002-6965-5211)

• Hans Dierckx¹ (ORCID: 0000-0003-0899-8082)

Institutions:

1. Department of Mathematics, KU Leuven Campus Kortrijk (KULAK), Etienne Sabbelaan 53, 8500 Kortrijk, Belgium

2. Laboratory of Experimental Cardiology, Leiden University Medical Center (LUMC), Albinusdreef 2, 2333 ZA Leiden, The Netherlands

3. iSi Health, Institute of Physics-based Modeling for In Silico Health, KU Leuven, Oude Markt 13, 3000 Leuven, Belgium

Correspondence: [email protected]

DOI: 10.1038/s41598-024-73544-z

Abstract:

Many extended chemical and biological systems self-organise into complex patterns that drive the medium behaviour in a non-linear fashion. An important class of such systems are excitable media, including neural and cardiac tissues. In extended excitable media, wave breaks can form rotating patterns and turbulence. However, the onset, sustaining and elimination of such complex patterns is currently incompletely understood. The classical theory of phase singularities in excitable media was recently challenged, as extended lines of conduction block were identified as phase discontinuities. Here, we provide a theoretical framework that captures the rich dynamics in excitable systems in terms of three quasiparticles: heads, tails, and pivots. We propose to call these quasiparticles ‘cardions’. In simulations and experiments, we show that these basic building blocks combine into at least four different bound states. By representing their interactions similarly to Feynman diagrams in physics, the creation and annihilation of vortex pairs are shown to be sequences of dynamical creation, annihilation, and recombination of the identified quasiparticles. We draw such diagrams for numerical simulations, as well as optical voltage mapping experiments performed on cultured human atrial myocytes (hiAMs). Our results provide a new, unified language for a more detailed theory, analysis, and mechanistic insights of dynamical transitions in excitation patterns.

Some of the most intriguing waves and patterns in nature arise from the spatial coupling between excitable elements. Examples include the emergence of intelligence in neural tissue (Cannon et al., 2014), epidemic spread across a population (Smith et al., 2002), intracellular waves (Lechleiter et al., 1991), chemical oxidation waves (Kapral & Showalter, 1995; Rotermund et al., 1990) and the emergence of cardiac arrhythmias in the heart muscle (Allessie et al., 1973; Gray et al., 1998). Our work was motivated by the case of heart rhythm disorders, that had a global mortality rate of in 2019, with increasing numbers every year. Hence, cardiac arrhythmias remain one of the largest causes of death worldwide (WHO, 2020).

While the local excitation dynamics can be measured and modelled, a firm theoretical description at the mesoscale is still lacking. Hence, much effort is currently invested in numerical modelling (Clayton et al., 2011) and machine learning approaches (Trayanova et al., 2021). Our aim is to complement such approaches with semi-analytical methods that give more insight in why patterns develop over time in a certain manner.

In this paper, we revise the classical mesoscale theory of excitation that has been used for cardiac arrhythmia analysis over the past decades (Christoph et al., 2018; Clayton et al., 2006; Clayton et al., 2011; Gray et al., 1998). Similar to historical developments in particle physics, we coin that the classical phase singularities in excitation are formed of more elementary quasiparticles. While this finer structure remains hidden in simple dynamical regimes, we show below in simulation and experiment that it becomes relevant during the processes of wave block, vortex initiation and recombination. Such insight is relevant for excitation patterns, since precisely those events need to be thoroughly understood if one wants to predict or control the emergent dynamics in those systems. A graphical outline of the analysis method presented in this paper can be found in Fig. 1.

particle	notation	simplified notation	charge $Q$	charge $P$	composition	colour code
head	$\mathrm{h}^{\pm}$	$\mathrm{h}^{\pm}$	$\pm\frac{1}{2}$	0
tail	$\mathrm{t}^{\pm}$	$\mathrm{t}^{\pm}$	$\pm\frac{1}{2}$	0
pivot	$\mathrm{p}_{\pm}$	$\mathrm{p}^{\pm}$	$0$	$\pm\frac{1}{2}$
filament; tip	$\mathrm{f}^{\pm\pm}$	$\mathrm{f}^{\pm}$	$\pm 1$	$0$	h$^\pm$ + t$^\pm$
core; phase singularity	$\mathrm{c}^{\pm\pm}_{\pm\pm}$	$\mathrm{c}^{\pm}$	$\pm 1$	$\pm 1$	h$^\pm$ + t$^\pm$ + 2p$_\pm$
shrink site	$\mathrm{s}^{\pm}_{\pm}$	$\mathrm{s}^{\pm}$	$\pm\frac{1}{2}$	$\pm\frac{1}{2}$	p$_\pm$ + t$^\pm$
growth site	$\mathrm{g}^{\pm\pm}_{\pm}$	$\mathrm{g}^{\pm}$	$\pm 1$	$\pm\frac{1}{2}$	h$^\pm$ + t$^\pm$ + p$_\pm$

Table 1: Overview of quasiparticles (cardions) in excitable media, together with their charges, notation, composition, and colour code. In the extended notation, superscript denotes the head/tail charge $Q$ in units of $\frac{1}{2}$ and subscripts denote the pivot charge $P$. As only quasiparticles of the same charge bind, the simplified notation indicates this with a single superscript $\pm$. The list is non-exhaustive and reflects our current knowledge.

As we illustrate our general theoretical framework below with an analysis of experiments performed on cultured human heart cells (Harlaar et al., 2022), we now describe the current state-of-art and associated challenges in the description and understanding of cardiac excitation patterns.

The mechanical contraction of each myocyte in the heart is triggered by electrical depolarisation of the cells, which first propagates through a dedicated conduction system, but is thereafter passed on from cardiomyocyte to cardiomyocyte. As a result, the activation wave can re-excite itself, or a local focal source can start extra excitations. Such abnormal emerging behaviour manifests as a cardiac arrhythmia.

A simple arrhythmia in which the excitation wave follows a fixed path in time, or which originates from a region of ectopic impulse formation can be treated by ablation therapy, in which the trigger or path is physically destroyed (Haissaguerre et al., 1998). Even though efficiency has improved significantly throughout the years, the success rate is rather low for arrhythmias with incompletely understood spatiotemporal organisation (Cronin et al., 2019; Mujović et al., 2017). Such arrhythmias include ventricular fibrillation (Gray et al., 1998; Jalife et al., 1998), which is the most lethal one, and atrial fibrillation, which affects approximately of the world population, and is responsible for one third of stroke cases (Samol et al., 2016). In both atrial and ventricular fibrillation, as well as in ventricular tachycardia, rotating vortices of activation have been observed (Allessie et al., 1973), named spiral waves, scroll waves, or cardiac rotors. While computer simulations of the heart (Gray et al., 1995) and initial observations (Gray et al., 1998) often yield long-lived stable rotors, complex patterns detected in real hearts commonly show a more complex regime that exhibits conduction blocks and transient rotors which do not even complete a full turn. There is a consensus that interacting wave fragments, also called wavelets (Arno et al., 2021; Lee et al., 2020; Moe & Abildskov, 1959) play a role in atrial and ventricular fibrillation. However, the precise mechanisms of onset, sustainment, and termination of complex arrhythmias remain to be further elucidated (Aras et al., 2017; Shibata et al., 2022).

We now review the concepts of phase and phase singularities that are used to describe excitable media. Winfree (1973) already noted that many biological processes are cyclic in nature, and therefore could be described naturally by a ‘phase’. Near the centre of a rotating spiral wave, phase singularities may occur (Winfree & Strogatz, 1983), which form filament curves in three dimensions. Given two observed variables $V$ and $R$ that depend on the position ${{\bm{{r}}}}$ and time $t$ in a spatiotemporal pattern, it is common to call the polar angle in the state space spanned by $V$ and $R$ the ‘phase of activation’: $$\begin{aligned} \varphi({{\bm{{r}}}},t) &= \operatorname{atan2}( R({{\bm{{r}}}},t) - R_*, V({{\bm{{r}}}},t) - V_*) \end{aligned}\qquad{(1)}$$ Here, $\operatorname{atan2}$ is the two-argument arc-tangent function, delivering values in $(-\pi,\pi]$. In a seminal work, Gray et al. (1998) analysed ventricular fibrillation patterns in this manner (Gray et al., 1998), with $V({{\bm{{r}}}},t)$ the optical intensity representing normalised transmembrane potential, and $R({{\bm{{r}}}},t) = V({{\bm{{r}}}},t+\tau_0)$ its time-delayed version. The main idea here is that points that differ in phase by an integer multiple of $2\pi$ are effectively in the same state. After spatial filtering of the signal, distinct points occurred where all phases meet, called phase singularities (Winfree & Strogatz, 1983). Since then, analysis tools (Bray & Wikswo, 2003; Gray et al., 1998), ablation strategies (Narayan et al., 2012), and mathematical theory (Clayton et al., 2006; Panfilov & Dierckx, 2017) have been developed based on the assumed existence of phase singularities during arrhythmias. A phase singularity can also be regarded as the intersection of two contour lines of two different state variables, which also allows to locate them numerically (Fenton & Karma, 1998). An extension of this topological approach that focuses on the dynamical transitions between the contour intersections was described by Marcotte & Grigoriev (2017).

Characterising excitable media in terms of phase singularities was pioneered based on observations in oscillating chemical reactions (Winfree, 1973; Zhabotinsky & Zaikin, 1973). When a phase singularity is tracked in computational models of cardiac excitation (Fenton & Cherry, 2008), however, it is often seen to follow a star-like or zig-zag pattern (Krinsky et al., 1992). Such a ’linear rotor core’ also occurs at the interface between regions with different action potential duration, where it was historically called a ‘reverbator’ (Krinsky, 1966, 1968). Linear core rotors have for decades escaped a thorough theoretical description, in contrast to the better-behaved ‘rigidly rotating rotors’. The latter are now well understood in terms of response function theory (Biktashev et al., 1994; Henry & Hakim, 2000; Keener, 1988; Verschelde et al., 2007), which enables to forecast their motion if they are long-lived and weakly interacting with other rotors, boundaries or inhomogeneities. However, in arrhythmia management, the key processes to control are the transient processes causing the ‘birth’ of rotors from wave breaks during arrhythmogenesis, and the elimination of rotors under pharmacological treatment, scarification of the medium, or by mutual collisions. Even worse, atrial and ventricular fibrillation are characterised by an incessant interaction between short-lived wavebreaks, for which a theory is hitherto largely lacking.

In recent years, it was suggested independently by Tomii et al. (2021) and Arno et al. (2021) that the classical phase singularity concept may not be an accurate description for the structure of a linear rotor core. Here, the zig-zag pattern emerges since the classical spiral wave tip, i.e., the end point of the wave front, is moving along the border of a region that has not yet recovered from the previous activation. As soon as this region recovers, the spiral wave tip turns over about in a so-called pivot point and then travels in the opposite direction, after which this process is repeated. The mentioned authors (Arno et al., 2021; Tomii et al., 2021) remarked that the parts of the medium at both sides of the linear rotor core have a different phase, as these regions were excited at different times. Hence, Tomii et al. (2021) coined that across the conduction block line at the centre of a linear-core rotor, there is a phase discontinuity. We also noted this (Arno et al., 2021), but remarked that due to electrotonic effects, i.e., diffusion of transmembrane potential, a boundary layer of finite width will be formed, across which the phase appears to change rapidly in space. The situation is reminiscent of domain walls or other interfaces in physics, and therefore we refer to them as ‘phase defects’ (Arno et al., 2021). The diffusive effect also explains why the phase defects were not identified before: If you have a line along which phase is discontinuous and smoothen the variables slightly, the resulting phase function tends to become continuous except in isolated points that happen to have the values $V=V_*$, $R=R_*$, where all phases meet. These points are exactly the phase singularities that have been used in cardiac analysis for nearly three decades now (Gray et al., 1995). While many of these points will be true phase singularities, a diffusive smoothing applied to a conduction block may introduce pairs of false phase singularities, which in our opinion explains the non-robust detection of phase singularity pairs in conduction block regions (Arno et al., 2021; Rodrigo et al., 2017).

The concept of a ‘phase defect’ in excitable media is able to link different objects already known in the field and will reveal new ones, as can be seen below. From the explanation above, the phase defect can represent a conduction block without a rotor attached to it, and the conduction block that is present in a linear-core rotor. In hindsight, Krinsky’s reverbators (Krinsky, 1966, 1968) were also phase defects, as across the interface, the phase was also discontinuous. Note that conduction blocks are key features in activation maps (Janse et al., 1980). In clinical acquisitions, it is common to record the local activation times (LATs) and draw the isochrones. Where different isochrones coincide, a conduction block and hence a phase defect is present. This situation is similar to cartography, where coincident isolines mark a cliff, i.e., a discontinuity in altitude.

The aim of this paper is to provide a more detailed theoretical framework that reconciles the concepts of wave fronts, wave backs, conduction blocks, rotors with different core types, fragmented wave fronts (wavelets) etc. The key concept will be the phase defect mentioned above, and we will show that tracking how such phase defect is connected to wave fronts and backs opens up new ways for pattern analysis. Our exposition has been designed to unravel the critical processes that cause, sustain or terminate complex arrhythmias. We purposely use terminology from particle physics, as in this work we bring several analogies from that discipline to a cardiac context.

We build up the framework by defining and motivating subsequently: phase, phase defects, quasiparticles, charges and diagrams. A sketch illustrating the main quantities in the following definitions can be found in Fig. 2.

Winfree wrote the following motivation to describe biological media with a phase concept (Winfree & Strogatz, 1983): “After excitation, [the phase] initially grows linearly with time, but after a while, it lingers near a value of $2\pi$, which is equivalent to phase 0, the completely recovered state.” So, when looking at an action potential as a function of time, we require that $\varphi(t) = 0$ for $t<0$ and $\varphi(t) = 2\pi$ for large $t$. Phase definition Eq. 1 is one possibility. It is certainly not unique, since different variables $V$, $R$ can be chosen, as well as their threshold values $V_*$, $R_*$. Two other choices were given in by Kabus et al. (2022), designed to exhibit no sharp phase variation around the wave front.

To make the link with isochrone maps, we here define the local activation time (LAT) as the last time when the transmembrane potential $V$ rose above a threshold value $V_*$: $$\begin{aligned} t_\mathrm{LAT}({{\bm{{r}}}},t) = \operatorname{argmax}_{t_0 < t} \{ V({{\bm{{r}}}},t_0) = V_* \wedge \partial_t V({{\bm{{r}}}},t_0) > 0 \} \end{aligned}\qquad{(2)}$$

The elapsed time since last activation is then $\tau = t - t_\mathrm{LAT}$. We now propose a linear activation time-based phase: $$\begin{aligned} \varphi({{\bm{{r}}}}, t) = \begin{cases} 2\pi \frac{\tau({{\bm{{r}}}}, t)}{\tau_0} & \tau({{\bm{{r}}}}, t) < \tau_0 \\ 2 \pi & \tau({{\bm{{r}}}}, t) \geq \tau_0 \end{cases} \end{aligned}\qquad{(3)}$$ Hence, the resting state has phase 0, and after activation phase grows linearly with time over a time interval $\tau_0$, after which the tissue is assumed to be fully recovered. Then the tissue has phase $2\pi$, which is equivalent to phase $0$. This simple phase has its limitations: It is identical to zero at the foot of the upstroke, when $V$ is still lower than $V_*$. After a time interval $\tau_0$, that we choose equal to the local action potential duration in that point, the medium is assumed to be fully recovered. Still, it suffices for the formal introduction of the concepts below.

By construction, $\varphi$ is continuous over a wave front, since there the phase amounts to $0$. It is also continuous over a wave back, where $\varphi$ is slightly lower than $2\pi$. At conduction block sites, however, the impeding wave does not make $V$ cross the threshold value $V_*$ from below anymore, and hence $\varphi$ will be discontinuous across a conduction block. As a consequence of the phase definition Eq. 3, a true phase discontinuity forms, of co-dimension 1. In an inhomogeneous medium where action potential duration depends on position, the edge of the inhomogeneity can take any shape, such that the discontinuity of $\varphi$ needs not be a line but can take any shape: Closed or branching phase defect structures may also form.

To numerically detect the phase defect, we compute the phase defect density $\rho({{\bm{{r}}}},t)$ over all pixels of the grid (Kabus et al., 2022; Tomii et al., 2021): $$\begin{aligned} \rho({{\bm{{r}}}}_a,t) = \frac{1}{2 N_a} \sum_{b \in \mathcal{N}(a)} [1 - \cos( \varphi({{\bm{{r}}}}_a,t) - \varphi({{\bm{{r}}}}_b,t) )] \end{aligned}\qquad{(4)}$$ Here, $\mathcal{N}(a)$ is the set of nearest neighbours in the Cartesian grid next to the point ${{\bm{{r}}}}_a$. $N_a$ is the number of such neighbours; it is everywhere four except near the domain boundaries or inexcitable obstacles. When the phases locally differ much, $\rho$ will have a value close to $1$. If all excitable elements excite simultaneously, $\rho=0$. The set $\mathrm{Z}$ of points inside a phase defect can be identified numerically as the set where $\rho$ exceeds a chosen threshold value $\rho_*$ (Kabus et al., 2022): $$\begin{aligned} \mathrm{Z}(t) = \{ {{\bm{{r}}}} \;|\; \rho({{\bm{{r}}}}, t) > \rho_* \} \end{aligned}\qquad{(5)}$$ In this manner, the phase defect region $\mathrm{Z}$ forms an open subset of the domain; it is not a one-dimensional curve anymore. The contours $\mathcal{C}$ obeying $\rho = \rho_*$ encircling each phase defect are well-defined and will be used below.

It is also possible to start from the classical activation phase, rescale it via a sigmoid function and then compute the phase defect density (Kabus et al., 2022). In that case, the phase defect also becomes an open set $\mathrm{Z}$, surrounded by a continuous contour $\mathcal{C}$.

The essence of our analysis is to divide the medium in three regions and to identify the points where all these regions meet, see Fig. 1.

First, we compute the phase defect density $\rho$ via Eq. 4 to find all points within the phase defect $\mathrm{Z}$, see Eq. 5. These points will be colored in gray.

Let $0<\varphi_\mathrm{F}<\varphi_\mathrm{B}<2\pi$, where $\varphi_\mathrm{F}$ is a small phase threshold close to $0$, and $\varphi_\mathrm{B}$ a large phase threshold close to $2\pi$, which must be small enough that a jump from $\varphi_\mathrm{B}$ to $\varphi_\mathrm{F}$ is considered a phase defect following Eq. 4 and Eq. 5. The set of points where the phase is between these values and outside a phase defect is the excited region ($\mathrm{E}$): $$\begin{aligned} \mathrm{E}(t) = \{ {{\bm{{r}}}}\;|\; \rho({{\bm{{r}}}}, t) < \rho_* \wedge \varphi({{\bm{{r}}}}, t) \in (\varphi_\mathrm{F}, \varphi_\mathrm{B}) \} \end{aligned}\qquad{(6)}$$ and is depicted in yellow in Fig. 1. The set of points outside of that range and phase defects is called unexcited ($\mathrm{U}$) and coloured white: $$\begin{aligned} \mathrm{U}(t) = \{ {{\bm{{r}}}}\;|\; \rho({{\bm{{r}}}}, t) < \rho_* \wedge \varphi({{\bm{{r}}}}, t) \not\in [\varphi_\mathrm{F}, \varphi_\mathrm{B}] \} \end{aligned}\qquad{(7)}$$

Where the excited region meets the unexcited one, the wave front ($\mathrm{F}$) and wave back ($\mathrm{B}$) are found, denoted in cyan and magenta, respectively: $$\begin{aligned} \mathrm{F}(t) = \{ {{\bm{{r}}}}\;|\; \rho({{\bm{{r}}}}, t) < \rho_* \wedge \varphi({{\bm{{r}}}}, t) = \varphi_\mathrm{F} \} \\ \mathrm{B}(t) = \{ {{\bm{{r}}}}\;|\; \rho({{\bm{{r}}}}, t) < \rho_* \wedge \varphi({{\bm{{r}}}}, t) = \varphi_\mathrm{B} \} \end{aligned}\qquad{(8)}$$

For completeness, we remark that two other boundaries of co-dimension 1 exist: the border between $\mathrm{Z}$ and $\mathrm{U}$, and the border between $\mathrm{Z}$ and $\mathrm{E}$. Together they form the aforementioned contours $\mathcal{C}$.

In the classical theory of excitation, the wave front meets the wave back in a single point, and this point coincides with the classical phase singularity (Clayton et al., 2006; Gray et al., 1995; Zykov, 1987).

In our case, wave fronts and wave backs will always end on a phase defects and medium boundaries. We call the point where a wave front and a phase defect meet, a ‘head’ (denoted in blue) and a point where a wave back and a phase defect meet, a ’tail’ (denoted in red). Essentially, heads and tails are the loci where the three distinct zones ($\mathrm{E}$, $\mathrm{U}$ and $\mathrm{Z}$) meet. As such, any continuous deformation of the boundaries between these zones will not affect the existence of such special points. Therefore, we call heads and tails topologically protected. A similar reasoning holds under time evolution: If the excited, unexcited and phase defect regions change their boundaries in a continuous manner, the head point can move in time, but is only lost if it meets another head or a medium boundary. An analogous reasoning justifies to call the tail points to be topologically conserved.

From Fig. 1E, we see that heads can exist with two different chiralities: If one circumscribes the head point counterclockwise and meets the regions $\mathrm{Z}$, $\mathrm{E}$, $\mathrm{U}$ in that order, we call it a ‘positive’ head, denoted $\mathrm{h}^+$. If we find the order $\mathrm{Z}$, $\mathrm{U}$, $\mathrm{E}$, we call it a ’negative’ head, written $\mathrm{h}^-$. These intuitive notions will be formalized below, to grant them the meaning of a ‘charge’.

For now, let us return to the depiction of a linear-core rotor in Fig. 1D. We observe yet another type of special point, namely the end point of the phase defect line. Such point are already known in literature as the loci around which the wave turns suddenly, and hence we also call them ‘pivots’. Using the LAT-based phase definition Eq. 3, it follows that the phase defects are truly curves of discontinuous phase. Hence these lines can have distinct end points, the pivots. The phase defect lines may also branch spatially; in such case a phase defect ‘joint’ will emerge, which however falls outside the scope of our present exposition. Pivots also exist in two flavours, see Fig. 1E: If the wave made a sudden half-turn to the left at that pivot, we call the pivot point a $\mathrm{p}_+$, with the plus sign marking counterclockwise rotation. If the half-turn was to the right, we call it a $\mathrm{p}_-$. We will discuss the topological preservation of pivots below, after the formal introduction of pivot charge.

We will here demonstrate that heads and tails carry half-integer charges that are identical to the classical ‘charge’ of a phase singularity. Pivots carry a charge that is also half-integer, but of a different kind.

First, consider a portion of the contour $\mathcal{C}$ that surrounds a phase defect in a counterclockwise sense, i.e., with the phase defect to the left. See Fig. 2. With arc length parameter $\ell$, we evaluate the classical charge density: $$\begin{aligned} Q = - \frac{1}{2\pi} \int d \varphi = - \frac{1}{2\pi} \int {{\bm{{\nabla}}}} \varphi \cdot {{\bm{{d \ell}}}}. \end{aligned}\qquad{(9)}$$ Here, $\varphi$ can be any phase related to the local state of cells; both Eq. 1 and Eq. 3 are possible, just like other choices (Clayton et al., 2006; Gray et al., 1995; Tomii et al., 2021). The minus sign is included such that a counterclockwise rotating spiral has $Q=+1$. Authors usually evaluate Eq. 9 around a closed contour, to find that $Q$ always takes an integer value if $\varphi$ is well-defined over the entire contour. However, if one takes only a portion of the contour, one can find the charge density $dQ = - \frac{d\varphi}{2\pi}$. Part of the phase change occurs during excitation, and part during the biological recovery process. If one assumes these changes to take place at the wave front and wave back, a phase change of $\pm \pi$ will be found each time the contour encircling the phase defect meets a wave front or wave back. These values correspond to localized charges of $\pm \frac{1}{2}$. Therefore, we here choose to ascribe a topological charge of $\pm \frac{1}{2}$ to heads and tails, with the sign chosen as in Fig. 1E. In principle, this choice can be made rigorous by defining the phase in the region $\mathrm{U}$ to be 0 and the phase over the region $\mathrm{E}$ to be $\pi$. Elegantly, the sum of all head and tail charges when circumscribing a phase defect equals the classical topological charge $Q$. To distinguish it from the pivot charge below, we call this the $Q$-charge. From the above, it follows that the $Q$-charge in a way ‘measures’ the amount of rotation in state space.

Second, we introduce the pivot charge as follows. Define the angle of wave propagation as: $$\begin{aligned} \psi = \operatorname{atan2}( \partial_y t_\mathrm{LAT}, \partial_x t_\mathrm{LAT}). \end{aligned}\qquad{(10)}$$ where the LAT is continuous and differentiable. $\psi$ is the angle with the positive $x$-direction under which the last arrived activation wave travelled. Now, in full analogy with Eq. 9, step along the contour $\mathcal{C}$ while keeping the phase defect to the left: $$\begin{aligned} P = \frac{1}{2\pi} \int d \psi = \int {{\bm{{\nabla}}}}\psi \cdot {{\bm{{d \ell}}}} \end{aligned}\qquad{(11)}$$

Let us verify how the $P$-charge density $d P = \frac{d\psi}{2\pi}$ behaves in standard cases. First, take the neighbourhood of the positive pivot at the top-right panel of Fig. 1E. There, the wave turns over $+\pi$ near the pivot, such that $P=+\frac{1}{2}$. This corresponds to the intuitive notion of ‘half a turn to the left’ that we advocated above. Similarly, a negative pivot $\mathrm{p}_-$ can now be proven to have $P=-\frac{1}{2}$. Third, consider a classical spiral wave core denoted as $\mathrm{c}^+$ in the bottom row of Fig. 1E. Near the classical phase singularity, phase gradients will become arbitrarily large, such that the classical phase singularity is also a special case of a phase defect: It is a phase defect of small size (Arno et al., 2021), whose exact diameter depends on the chosen $\rho_*$. Evaluating $P$ around a closed contour that circumscribes the phase singularity counterclockwise delivers: $P=+1$. Noteworthily, even circular rotor cores form over time from conduction blocks: Initially, a finite conduction block line is created, e.g. by an S1-S2 stimulus protocol, with two pivots each having $P=+\frac{1}{2}$. During successive wave rotations, this phase defect shrinks to a point. At that time, the two pivots with $P=+\frac{1}{2}$ as well as the head and tail with each $Q=+\frac{1}{2}$ merge to form a stable state, with $P=+1$ and $Q=+1$, that we call a ‘core particle’ here.

Let us now investigate when the $P$-charge is conserved. Hereto, we apply Stokes’ theorem to Eq. 11 on a closed contour $\mathcal{C}$ around a region $S$. $$\begin{aligned} P = \frac{1}{2\pi} \oint d \psi =\oint {{\bm{{\nabla}}}}\psi \cdot {{\bm{{d \ell}}}} = \iint {{\bm{{\nabla}}}} \times {{\bm{{\nabla}}}}\psi \cdot {{\bm{{d s}}}} =0 \end{aligned}\qquad{(12)}$$ The last equality follows from the properties of curl and gradient. The second-last equality requires that the conditions for Stokes’ theorem are fulfilled: That ${{\bm{{\nabla}}}}\psi$ is everywhere in $S$ defined, and that it has continuous first derivatives. Result Eq. 12 implies that $P$ will not change its value if the contour is continuously deformed outside regions with phase defects or phase singularities (where $\psi$ is undefined), or where $\psi$ suddenly changes its orientation. The latter case happens at excitation sources, or on the boundary at the position where it was hit by a wave front. Given these limitations, the $P$-charge is conserved. E.g. a closed contour may contain two opposite rotating linear-core rotors, i.e., with $2\mathrm{p}_+$ and $2\mathrm{p}_-$ in it. The total $P=\frac{1}{2}+\frac{1}{2}-\frac{1}{2}-\frac{1}{2}=0$. If these rotors merge to form a single conduction block line, a $\mathrm{p}_+$, $\mathrm{p}_-$ pair annihilates, leaving the second $\mathrm{p}_+, \mathrm{p}_-$ inside the contour, and still $P=\frac{1}{2}-\frac{1}{2}=0$.

In conclusion, the $P$-charge measures the amount of rotation in physical space. It is thus different from the $Q$-charge, carried by heads and tails. Therefore, when using the detailed notation for charges (see Table 1), we denote the $P$-charge in subscript: $\mathrm{p}_\pm$.

Until here, we have identified three different special points in excitable media, that can carry half-integer ‘charges’ and that seem to be preserved under time evolution (for now disregarding boundaries and obstacles). We propose to refer to these special points as quasiparticles and call them ‘cardions’. Examples of the seven different cardions that we know now of are given in Fig. 1E. While three cardions are elementary, heads, tails and pivots, four of them arise as bound states between the elementary ones.

Strikingly, our analysis of activation patterns in simulation and experiments revealed that the heads, tails and pivots tend to form bound states, see Fig. 1E and Table 1. The situation reminds us of the hydrogen atom consisting of a proton and an electron, or the proton itself, which is made up of quarks in particle physics.

The first bound state consists of a head and a tail of equal $Q$-charge. Already in Fig. 1D, it can be seen that while the tail traces a phase defect, a head can follow it directly, exciting the newly-recovered tissue created by the tail particle. The two special points propagate as a pair along the phase defect, with an excitable gap between them. In the classical viewpoint, this corresponds to a wave front and wave back ending in one point, the classical ‘spiral tip’. In three spatial dimensions, the tips form a filament curve (Clayton et al., 2006; Winfree et al., 1996) around which the wave rotates. Therefore, we annotate the bound state of head and tail with f and call it a filament point and annotate it in the colour purple. The convention of assigning half a $Q$-charge to heads and tails has the useful property that the sum of these charges will be equal to the classical topological charge $Q=\frac{1}{2}+ \frac{1}{2} = 1$ of the filament. As the absolute charge of a filament particle is always the same, we abbreviate the ‘full’ notation f$^{\pm \pm}$ to $\mathrm{f}^{\pm}$, see second column in Table 1. In this simpler notation, we write the charges of all particles with superscript, for example, $\mathrm{p}^+$, rather than $\mathrm{p}_+$.

The second bound state is the core particle, which we introduced above as an example of $P$-charge. This core particle occurs with rigidly rotating spiral waves, where the phase defect is small and head and tails are close to each other. Since these cores typically form from shrinking of a phase defect line containing two pivots of the same charge, the core particle consists of a head, a tail and two pivots (where the latter are no more distinguishable). This ‘core particle’, denoted $\mathrm{c}^{\pm\pm}_{\pm\pm}$ or shorter $\mathrm{c}^{\pm}$ and indicated by an orange circle, see the bottom panel of Fig. 1E.

The next two bound states have no classical counterpart to our knowledge. Consider a propagating wave front that hits a wave back. Such event will generally occur in a single point first, creating a small phase defect. The end point of one side of the wave front will then travel along the existing wave back, and in the process, the conduction block region grows to become a phase defect line. Hence, there exists a state of a head, tail, and pivot of the same chirality that all move together. We call this a ‘growth’ particle and denote it with a teal circle. It has $Q=\pm 1$ and $P=\pm \frac{1}{2}$, with the same sign chosen. Hence, we identify quasiparticles g$^{++}_+$ and $\mathrm{g}^{--}_{-}$, or $\mathrm{g}^+$, $\mathrm{g}^-$ in shortened notation (see Table 1).

Finally, we were ourselves struck by the observation that the shrinking and growing of phase defect lines occurs in an asymmetric manner. When a wave back travels along a phase defect such that the tissue has recovered at both sides of the defect line, the phase defect line will retract, see Fig. 1E, bottom right. Such retraction implies that a tail and a pivot of same chirality travel together. We denote this as a ‘shrink’ particle s$^\pm_\pm$, shortened s$^\pm$, and mark it in brown. Table 1 summarises the hitherto discovered quasiparticles and their bound states in excitable media.

Until here, we have disregarded inexcitable obstacles in the medium, which also includes the medium boundary itself.

Let us first consider how the $Q$-charge behaves in the presence of obstacles. In the example of Fig. 1D, the total $Q$-charge in the medium disregarding the medium boundaries would be $-\frac{1}{2}-\frac{1}{2} = -1$. However, we now argue that it might be useful to also regard medium boundaries as a kind of phase defect. If we do this, e.g. by colouring the medium boundaries also in grey, it follows that every wave front is either closed or has two endpoints (heads) with different chirality. A similar reasoning holds for the wave backs, which are either closed, or end in two tails of opposite charge. Then, by construction the total $Q$-charge for both all heads in $H$, the set of heads, and tails in $T$, the set of tails, must be zero: $$\begin{aligned} \sum_{h \in H} Q_h &= 0 \\ \sum_{t \in T} Q_t &= 0 \end{aligned}\qquad{(13)}$$ These relations only hold if medium boundaries are considered as phase defects.

Taking the sum of the previous two statements proves that the total $Q$-charge in a two-dimensional medium is also always zero: $$\begin{aligned} \sum_{j \in H \cup T} Q_j &= 0 \end{aligned}\qquad{(14)}$$ considering medium boundaries as phase defects. This result is stronger than the conservation of phase singularities, which only holds in the bulk of the medium.

In the viewpoint that includes medium boundaries and inexcitable obstacles, the classical annihilation of a rotor core at the medium boundary can be regarded as the head-tail pair becoming a pair of a head and tail which are travelling apart from each other on the boundary. Also, the generation of a single rotor via S1-S2 stimulation now includes the boundary: The S2 pulse introduces a pair of heads, one of which joins a tail particle of the S1 pulse to form a filament quasiparticle.

An additional reason for including the medium boundary in the description lies in the interpretation of electrical signals. In leading order, the unipolar electrical potential measured by an electrode is determined by the angle subtended by the wave fronts and wave backs, seen from the electrode (Holland & Arnsdorf, 1977). The measured signal is thus only determined by the boundaries of the wave fronts and backs, namely the heads and tail curves in three dimensions (Arno et al., 2024b). Therefore, the concepts of heads and tails may also prove useful when improving efficient inverse methods to recover spatiotemporal activation sequences from patient measurements (Ramanathan et al., 2004).

Second, we consider how the pivot charge behaves in the presence of obstacles. If a wave front travels clockwise along the edge of a square domain, it will turn 90 ° to the right when crossing the corner. Thus, according to Eq. 11, the corners of a square domain carry a pivot charge of $\pm\frac{1}{4}$. Since the edge of obstacles can turn over any angle, pivot charges are not always half-integer, but can take on the boundary of the medium arbitrary values in $(-1,1)$. Another example is a wave hitting an obstacle within the medium, as shown in Fig. 3. Suppose the wave encircling the obstacle counterclockwise travels over an angle $\psi_\mathrm{l} >0$, while the wave travelling in the other direction around the obstacle over an angle $\psi_\mathrm{r} < 0$, with $|\psi_\mathrm{l}| + |\psi_\mathrm{r}| = 2\pi$ . Then, the $P$-charge of this obstacle will become via Eq. 11: $P = \frac{\psi_\mathrm{l} + \psi_\mathrm{r}}{2\pi} \in [-1, 1]$. In this situation, the $P$-charge is not conserved. Nonetheless, the value of the integral Eq. 11 contains useful information: If it has $P = \pm 1$, a rotor is attached to the obstacle, while $P = 0$ indicates a symmetric conduction block. Intermediate values of $P$ indicate whether the wave-obstacle interaction is close to the formation of a rotor attached to the obstacle or not.

In summary, $P$-charge is not conserved in the presence of medium boundaries. This was also the case for the $Q$-charge with classical phase singularities, as they can simply vanish by collision with the medium boundary.

In the last step of our conceptual framework, we describe the evolution of cardions via diagrams similar to Feynman diagrams (Feynman, 1949). A first example is given in Fig. 1F. A $\mathrm{g}^+$, $\mathrm{g}^-$ pair is born, after which the growth particles each split into triples of head, tail, and pivot. In this graphical depiction, time passes from left to right and positions in space are represented along a single, vertical axis. The arrows represent the continuous existence of the quasiparticle in space and time, until an event happens in which the particle merges, splits or annihilates. The premise of this work is that these diagrams can offer insight in the pattern’s dynamics, see the examples below.

Fig. 4 shows a numerical simulation in the smoothed Karma model (Karma, 1993; Marcotte & Grigoriev, 2017) for atrial fibrillation, in which a single rotor spontaneously breaks up into multiple vortices. We analyze the creation of a vortex pair, known in cardiology literature as a figure-of-eight re-entry (Wu et al., 1994). In the top two rows, phase singularities detected by the Kuklik method (Kuklik et al., 2016) are annotated on top of the normalised transmembrane voltage and phase. The vortex pair with labels 2 and 3 is situated away from a central persisting rotor, labelled 1. In the third row of panels, heads, tails and pivots, as well as their bound states are annotated. When the wave front hits the wave back, a pair of $\mathrm{g}^+$, $\mathrm{g}^-$ emerges. Each particle lasts until the vortex makes a pivoting turn around the formed conduction block line, which marks the decay of the $\mathrm{g}$-particles into triples of head, tail, and pivot. Just after $t={465~\mathrm{{m}{s}}}$, the tail pair annihilates, as the wave back detaches from the conduction block line to form a single wave back again. The recovery of the tissue at the location of the initial collision marks the splitting of the phase defect line, which involves the birth of a shrink pair, just before $t={510~\mathrm{{m}{s}}}$. Thereafter, both phase defect lines shrink until the head, tail, and two pivots are so close to each other that they effectively form a core particle. We conclude that the birth of a rotor pair, which is a single event in the classical theory, in fact consists of seven subsequent events, revealed by the vertices in the diagram at the bottom of Fig. 4.

The creation of a figure-of-eight spiral pair in a monolayer of human atrial myocytes (Harlaar et al., 2022) is thoroughly analyzed via the quasiparticles and diagrams in Figs. 5-9, see also the Methods section below. By adding a voltage-sensitive dye to the cell culture, the local transmembrane voltage could be optically mapped (Salama et al., 1987), which is shown here after normalization to $[0,1]$. The culture was stimulated at a high frequency from the top left corner. After several transient wave breaks, i.e., conduction blocks, a pair of oppositely rotating spiral waves was formed. The experiment presents an in-vitro realization of the onset of an arrhythmia, an event to be prevented in patients.

Fig. 5 presents a diagrammatic overview of the figure-of-eight spiral pair formation, while Figs. 6-9 present parts of this process in more detail in the same style as Fig. 4.

As can be seen in Fig. 6, the initial conduction block produces a growth pair just before $t={4056~\mathrm{{m}{s}}}$, in the same way as in the in-silico case above, which then break apart into a head, tail and pivot each. The tails and heads annihilate and leave the persisting pair of pivots behind, marking the end points of the conduction block line. Later on, a secondary conduction block forms in the same way, also producing a pivot-pair, Fig. 5C.

The process causing the the mere conduction block to develop into a dangerous figure-of-eight reentry, which can be life-threatening in peoples’ hearts, can be seen in Fig. 7 After the wave back hits the phase defect region of the conduction block, a tail pair is produced. Each tail then meets the persisting pivot particles to produce shrink particles that begin to shorten the conduction block line at $t={4164~\mathrm{{m}{s}}}$. However, shortly afterwards at $t={4170~\mathrm{{m}{s}}}$, the wave front also hits the phase defect, leading to a pair of heads, that combine with the shrinks to growth particles. Driven by this pair, the conduction block line grows again, until the growth particles decay into a head, tail and pivot each. Just after $t={4194~\mathrm{{m}{s}}}$, the tails annihilate, as well as the heads just after $t={4212~\mathrm{{m}{s}}}$, while again two pivot particles persist. As the growing outweighs the shrinking in this case, in total, the phase defect length accumulated. The distance between the pivots increased from ${(2.85\ ±\ 0.50)~\mathrm{{m}{m}}}$ at $t={4146~\mathrm{{m}{s}}}$ to ${(4.65\ ±\ 0.50)~\mathrm{{m}{m}}}$ at $t={4212~\mathrm{{m}{s}}}$. This process repeats six times, each time making the phase defect longer.

At $t={4854~\mathrm{{m}{s}}}$, also the second phase defect has formed, and both will merge, as visualized in Fig. 8. The merging process starts when both defects grow with a pair of $\mathrm{g}^+$ and $\mathrm{g}^-$ each. Then, the growth particles decay in triplets of head, tail and pivot and the tail-pairs of both phase defect sites annihilate just as in the accumulation stage Fig. 7). Meanwhile, the $\mathrm{h}^-$ of the large phase defect and the $\mathrm{h}^+$ of the smaller phase defect annihilate with each other just after $t={4866~\mathrm{{m}{s}}}$, as well as their respective $\mathrm{p}^-$ and $\mathrm{p}^+$, merging the two phase defects into one larger U-shaped phase defect line. At $t={4890~\mathrm{{m}{s}}}$, a pair of pivots and heads each with opposing charges remain. The pivots now have a distance of ${(6.72\ ±\ 0.50)~\mathrm{{m}{m}}}$.

Finally, in Fig. 9, it can be seen how the recovery of the initial zone of contact creates a shrink pair, which causes the splitting of the conduction block line just before $t={4902~\mathrm{{m}{s}}}$. Shortly after that, each shrink particle decays into a pivot and tail. At $t={4932~\mathrm{{m}{s}}}$, the system has evolved to a state with two nearly parallel conduction block lines, i.e., phase defects that comprise an effective functional isthmus around which a figure-of-eight rotor pair can revolve. Each of the phase defects have a head, a tail, and two pivot particles, all with the same charge. The phase defect with a clockwise rotating rotor has $Q=-1$ and $P=-1$, while the counter-clockwise rotor has $Q=+1$ and $P=+1$. From $t={4956~\mathrm{{m}{s}}}$, the pivots remain in place and the head and tail propagate along the phase defect line. Zooming out, these particles behave like a pair of cores.

In brief, the quasiparticle framework here reveals that the establishment of the figure-of-eight re-entry is a multi-step process that takes over ${800~\mathrm{{m}{s}}}$ to take place. The process features tight interplay between dynamically changing conduction blocks and wave front ends (heads). The formalism unifies the concepts of rotors and conduction block lines, and reveals how an initially flat conduction block line first takes a U-shape ($t={4890~\mathrm{{m}{s}}}$) and then becomes a functional isthmus. The details of this interaction are not visible in the classical phase singularity analysis, which would only reveal the creation, vanishing, and recreation of a phase singularity pair.

The above analysis and examples show that complex excitation patterns feature certain preserved special points. Noteworthily, these are critical points in space, rather than critical points of evolution equations, which also satisfy topological rules (Glass, 1975). While tips, pivots and cores have been observed before, the growth and shrink quasiparticles emerge from the framework itself. This situation is to a certain degree reminiscent of the ‘Standard Model’ in physics, which is based on a small set of fundamental particles (Mann, 2010) and explains almost all experimental results in particle physics. The theoretical framework presented in this work gains its power from simple geometrical arguments, as the heads and tails are points where three zones meet: excited, unexcited, and phase defect. The concept of ‘head’ itself has been implicitly used before, e.g. in Krinsky’s quantitative theory of reverbators (Krinsky, 1968; Krinsky et al., 1992) for inhomogeneous domains. Our framework is much wider valid, as it can describe long and short-lived rotors, wavelets and conduction blocks in homogeneous and inhomogeneous excitable media.

In physics, it is customary to give localised and preserved structures a name with the suffix ‘-on’, with the first part of the name referring to the nature of the system of the (quasi-)particle (Merriam-Webster, n.d.). Examples include not only most particles in the Standard Model (proton, electron, baryon, fermion, photon) but also other moving structures such as phonons, solitons and excitons. In analogy, we here propose to name the here-identified quasiparticles as ‘cardions’, as they were first observed in cardiac excitation patterns.

The motivation for this study was to investigate the interactions between short-lived rotors and approaching wave fronts during arrhythmogenesis (Clayton et al., 2006). We go beyond the observation that apparent phase singularities co-locate with conduction blocks (Arno et al., 2021; Rodrigo et al., 2017), by looking at a finer spatial scale. We discovered multiple interactions of the cardions in the data and represented them schematically. Our initial analyses show new mechanistic insight, like the dynamical isthmus formation in Fig. 9. It can be expected that the automated construction of those diagrams for larger fibrillation domains (see Fig. 5) will allow an enhanced classification and characterisation of those patterns.

The quasiparticle viewpoint allows us to introduce more concepts from particle physics in the world of excitable media. First, note that the terms ‘creation’ and ‘annihilation’ also applied to the classical viewpoint (Clayton et al., 2006; Gray et al., 1998). There, a single quasiparticle was present, i.e., the spiral wave centre (phase singularity), bearing a positive or negative charge ($Q=\pm 1$) depending on its rotation sense. In our viewpoint, creation and annihilation also occurs with heads, tails and pivots. Since cardiac tissue is, in reality, non-homogeneous, a wave front will often create a head pair near obstacles, e.g. areas of fibrosis, which annihilates shortly afterwards. Second, the difference between linear-core and circular-core rotors follows from the mutual interaction between pivots, which can be expressed with spiral wave response functions (Biktasheva & Biktashev, 2003; Dierckx et al., 2017; Marcotte & Grigoriev, 2016). If nearby pivot charges of the same charge attract each other, they form a core particle. Conversely, when they repel each other, a linear-core rotor emerges. However, such a phase defect cannot grow larger than the width of a travelling pulse in the medium, since otherwise, the defect line will break by the creation of a new $\mathrm{p}_+$, $\mathrm{p}_-$ pair. Therefore, pivot particles are in our opinion subjected to confinement, analogous to quarks in quantum chromodynamics (K. G. Wilson, 1974). This confinement may partially explain why the here-identified quasiparticles remained unnoticed during simple dynamical regimes, i.e., at time intervals where no new vortices or conduction blocks form or disappear.

Finally, we note that in Fig. 1B, there are two apparent symmetries. First, if one exchanges excited with unexcited regions, wave fronts become wave backs, heads become tails and vice versa. This corresponds to flipping the voltage axis around the value $V_*$ from Eq. 1. Second, one can reverse the arrow of time, which preserves $\mathrm{U}$ and $\mathrm{E}$ regions, but interchanges wave backs with fronts and heads with tails. However, when looking at the dynamics, there is a clear distinction between heads and tails. As heads lie on a wave front, they move at the conduction velocity in the medium. Since the repolarisation of the tissue occurs at a fixed delay after the depolarisation, the wave back is a phase wave, and wave backs and tails can move at any speed. Also, the medium can remain in the unexcited state for a unlimited time but only a limited time in the excited state. This effective breaking of symmetry also manifests in the growth and decay states. In the cases we have yet observed in simulations and experiment, the growth of a phase defect line each time involved the three different quasiparticles. Namely, a conduction block arises when a wave front (h) hits a wave back (t) and in the process, a phase defect is formed (p). Conversely, the vanishing of a phase defect only requires that the phases at both sides return to full recovery. We took the threshold of the wave back to mark the state where a standard impeding wave front can re-excite the tissue. Then, only a tail and pivot particle are needed. In case of a heterogeneous medium, we anticipate that also other situations can occur, but these fall outside the scope of our present investigation. Noteworthily, a statistical analysis of occurring processes in the topological approach of Marcotte & Grigoriev (2017) also demonstrated asymmetry between favoured processes in the creation and annihilation processes of spiral cores.

The parallels drawn here between biological excitation and theoretical physics are only at the conceptual level. Feynman diagrams in physics are representations of a path integral, allowing to structure tedious calculations. To use the diagrams in excitable context for quantitative calculations and predictions could be a next research step. Even then, the mathematical implications will be different from the physics case. For this reason, we refer to the diagrams here consistently as ‘Feynman-like’. Likewise, the terms of confinement, quasiparticle interaction, standard model and symmetry breaking are powerful analogues when used in complex systems, but will at a more detailed level have different meaning than their original physics use.

The local state of a cell was here assumed to lie along a cycle, such that it can be labelled with a single phase variable $\varphi({{\bm{{r}}}},t)$. Capturing specific memory effects, such as alternans (Karma, 1994) or after-depolarisations in cardiac tissue (Weiss et al., 2010) will require an extension of the formalism, e.g. using phase-amplitude equations (D. Wilson, 2020).

We are also not certain if all possible cardions states have been described here. We already noted that phase defects can have side branches, and the ‘joints’ where this happens will also be cardions, to be described elsewhere. In Fig. 8, co-localisation of a $\mathrm{h}$ and a $\mathrm{p}$ particle is seen in the frames. However, in our experience this state is transient and therefore not considered a bound state here. Nonetheless, it is possible that other regimes exist where $\mathrm{h}$-$\mathrm{p}$ pairs may travel together as a bound state.

The concept of $P$-charge was introduced in this work, and our understanding of it is currently incomplete. The $P$-charge is currently undefined before the first wave crosses the medium. The precise conditions under which it is conserved or not need to be further elucidated. For isolated, thin, conduction block lines, $P$-charge seems to be conserved: Basic features in the patterns observed are linear rotor cores ($P=\pm 1$) or mere conduction blocks ($P=0$).

Finally, we have applied the framework here already on in-silico and in-vitro datasets. An obvious next step would be to also analyse surface patterns in cardiac tissues, and to post-process phase maps or LAT maps obtained in patients. We encourage the scientific community to work together on these steps and are providing basic numerical methods for this on GitLab, see below. We are currently working on automated algorithms to detect cardions in datasets.

The above examples revealed dynamics at a finer level than in the phase singularity picture. Conducting similar analyses on existing datasets during arrhythmias is likely to give new clues on underlying mechanisms. More specifically, counting the elementary interactions (Marcotte & Grigoriev, 2017), computing statistics of potentially large diagrams, and deriving interaction laws between the quasiparticles could prove useful. The value of the $P$-charge for obstacles in the medium could also serve as a ‘distance’ to the rotor regime, and thus quantify how far a regime is from arrhythmia, complementing other risk scores in literature (Arevalo et al., 2016).

The pivot sites, here identified as quasiparticles, are actively being investigated by clinicians as possible ablation targets (Seitz et al., 2017). We believe that the analysis of the clinical data and the pivot ablation within our framework could help to answer if, why, and when these points are appropriate ablation targets.

In this manuscript, we have exposed several new concepts, such as the pivot charge, different quasiparticles and a diagrammatic approach to better understand their interactions. Continuing research by ourselves and colleagues will be needed to elaborate these concepts.

Within the context of continuous excitable media, several extensions are possible. Since we do not suppose an underlying evolution model, the analysis could also be performed to the surface of three-dimensional media. A next step is to analyse optical voltage mapping surface recordings of the heart muscle. Such analysis can be applied to any excitable surface pattern that is sufficiently sampled in space and time. Our original motivation is to perform predictive calculations in those patterns: Can one find the critical size of a wave break, or the timescale at which a rotor produces multiple wavelets? In three spatial dimensions, our preliminary findings (Arno et al., 2024b) indicate that the cardions become string-like objects (Verschelde et al., 2007). As such, cardions in three dimensions will refine the concept of a rotor filament (Clayton et al., 2006; Winfree, 1973), and may be used in the future to further investigate the three-dimensional organisation of turbulence in complex systems.

We presented a conceptual framework revising the classical theory of excitable media, intentionally designed for complex regimes: short-lived pivoting motion and prominent conduction blocks. By classifying the endpoints of wave fronts, wave backs, and conduction blocks as charged quasiparticles, Feynman-like diagrams can be created from phase maps or local activation times. We believe this framework has the potential to become a useful analysis tool in excitable media, with applications within and beyond the cardiac electrophysiology context.

The synthetic data from Fig. 1A were obtained by Euler-forward stepping of the Bueno-Orovio-Cherry-Fenton model (Bueno-Orovio et al., 2008), with finite-differences grid size ${0.5~\mathrm{{m}{m}}}$ and time step ${0.1~\mathrm{{m}{s}}}$ in a biventricular heart geometry using an S1S2 stimulation protocol in the Ithildin framework (Kabus et al., 2024). The data from Fig. 1B are obtained in the same way but for a square piece of tissue with grid size ${0.3~\mathrm{{m}{m}}}$. For the fibrillation-like data from Fig. 4, the Karma model (Karma, 1993, 1994) modified according to Marcotte & Grigoriev (2017) was used, with step size ${1~\mathrm{{m}{m}}}$ and time step ${0.1~\mathrm{{m}{s}}}$.

Monolayers of fully functional human atrial myocytes were generated as described in detail by Harlaar et al. (2022). A voltage-sensitive dye was added to the culture, after which a real-time recording can be made of the intensity of emitted light, which is a measure of the local transmembrane potential (Salama et al., 1987). The recording has a resolution of $100 \times 100$ pixels with pixel size ${0.25~\mathrm{{m}{m}}}$. The sampling time between frames was ${6~\mathrm{{m}{s}}}$.

Local activation times were converted to phase following Eq. 3 with time constant $\tau_0 = {110~\mathrm{{m}{s}}}$ for the in vitro data and ${70~\mathrm{{m}{s}}}$ for the in-silico data. We use the phase values $\varphi_\mathrm{F} = {0~\mathrm{rad}}$ and $\varphi_\mathrm{B} = {5~\mathrm{rad}}$. Phase defects are identified following Eq. 5 with the threshold $\rho_* = 0.2$. The positions of tips, heads, tails, pivots and compound particles were drawn manually on the resulting figures.

The simulation output and pre-processed optical voltage mapping data we have used for this article can be found on Zenodo (DOI: 10.5281/zenodo.13379783). This repository also contains the scripts and Python modules to generate the figures found in this article. The source code for the Ithildin finite-differences solver (Kabus et al., 2024) (https://gitlab.com/heartkor/ithildin), its Python module (Kabus et al., 2022) (https://gitlab.com/heartkor/py_ithildin), and a Python module to create the Feynman-like diagrams (https://gitlab.com/heartkor/pdl-feynman) can also be found on GitLab.

Acknowledgments:

The authors are grateful to D.A. Pijnappels, A.A.F. de Vries and N. Harlaar for collecting and sharing the hiAM dataset. The authors thank J. Ector and T. De Coster for helpful discussions. The authors thank A.A.F. de Vries for useful comments on the manuscript. The authors thank A. Gobeyn for coining the term ‘cardion’ for the collection of quasiparticles.

Funding:

L. Arno is funded by a FWO-Flanders fellowship, grant 117702N. D. Kabus is supported by KU Leuven grant GPUL/20/012. H. Dierckx was supported by KU Leuven grant STG/019/007.

Competing interests:

The authors declare no competing interests.

Copyright:

Open Access: This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit https://creativecommons.org/licenses/by/4.0/.

Author contributions:

LA: Conceptualisation, Methodology, Software, Validation, Formal analysis, Investigation, Data curation, Writing – original draft, Writing – review & editing, Visualisation. DK: Conceptualisation, Methodology, Software, Validation, Formal analysis, Investigation, Writing – original draft, Writing – review & editing, Visualisation. HD: Conceptualisation, Methodology, Resources, Writing – original draft, Writing – review & editing, Supervision, Project administration, Funding acquisition.