Recurrence quantification analysis for fine-scale characterisation of arrhythmic patterns in cardiac tissue

1 Vol.:(0123456789) Scientific Reports | (2023) 13:11828 | https://doi.org/10.1038/s41598-023-38256-w www.nature.com/scientificreports Recurrence quantification analysis for fine‑scale characterisation of arrhythmic patterns in cardiac tissue Radek Halfar 1*, Brodie A. J. Lawson 2,3, Rodrigo Weber dos Santos 4 & Kevin Burrage 2,5 This paper uses recurrence quantification analysis (RQA) combined with entropy measures and organization indices to characterize arrhythmic patterns and dynamics in computer simulations of cardiac tissue. We performed different simulations of cardiac tissues of sizes comparable to the human heart atrium. In these simulations, we observed four classic arrhythmic patterns: a spiral wave anchored to a highly fibrotic region resulting in sustained re‑entry, a meandering spiral wave, fibrillation, and a spiral wave anchored to a scar region that breaks up into wavelets away from the main rotor. A detailed analysis revealed that, within the same simulation, maps of RQA metrics could differentiate regions with regular AP propagation from ones with chaotic activity. In particular, the combination of two RQA metrics, the length of the longest diagonal string of recurrence points and the mean length of diagonal lines, was able to identify the location of rotor tips, which are the active elements that maintain spiral waves and fibrillation. By proposing low‑dimensional models based on the mean value and spatial correlation of metrics calculated from membrane potential time series, we identify RQA‑based metrics that successfully separate the four different types of cardiac arrhythmia into distinct regions of the feature space, and thus might be used for automatic classification, in particular distinguishing between fibrillation driven by self‑sustaining chaos and that created by a persistent rotor and wavebreak. We also discuss the practical applicability of such an approach. Cardiac arrhythmias are a prevalent disease worldwide and a leading cause of death. The onset of these arrhythmias is complex and remains difficult to anticipate. In 1952, Hodgin and Huxley presented the first mathematical model of cell action potential (AP)1, which is a critical component of furthering our understanding of arrhythmia in a mechanistic way2,3. Some more modern examples include the pairing ofsimulations of mathematical models of cardiac activity with data to inform the study, and treatment, of atrial flutter4,5 and atrial fibrillation5–7. A particular topic of interest is ablative surgery, in which small scars on the heart muscle are created, blocking off regions thought to be producing or sustaining arrhythmogenic dynamics4,8. Given this treatment’s inconsistent success rate and invasive nature, there has been much interest in learning how to identify regions for ablation. One common idea is to target complex fractionated atrial electrograms (CFAEs), especially as a supplementary treatment if the pulmonary veins (a common trigger zone for atrial fibrillation) have already been electrically isolated9. Precise definitions of CFAEs vary, but the goal is to identify regions of conduction slowing/block, anchor points for re-entrant activation, sites of wavebreak or other pro-arrhythmic phenomena from electrograms. The issue is that simple criteria for identifying CFAEs are also prone to incorrectly labelling large regions of tissue not critical for atrial fibrillation maintenance as dangerous10. A natural path forward, then, is to improve the identification of ablation targets via more sophisticated analysis of electrograms. In this work, we use recurrence quantification analysis (RQA) combined with entropy measures and organization indices to characterize arrhythmic patterns in computer simulations of cardiac tissue. RQA has emerged as a tool with great potential in electrophysiology, owing to its established history in capturing the regular movements of nonlinear dynamical systems, identifying transitions between regular and OPEN 1IT4Innovations, VSB - Technical University of Ostrava, 708 00 Ostrava, Czech Republic. 2ARC Centre of Excellence for Plant Success in Nature and Agriculture, Queensland University of Technology, Brisbane 4000, Australia. 3Centre for Data Science, Queensland Univeristy of Technology, Brisbane 4000, Australia. 4Graduate Program in Computational Modeling, Universidade Federal de Juiz de Fora, Juiz de Fora 36036-330, Brazil. 5Department of Computer Science, University of Oxford, Oxford, UK. *email: [email protected] 2 Vol:.(1234567890) Scientific Reports | (2023) 13:11828 | https://doi.org/10.1038/s41598-023-38256-w www.nature.com/scientificreports/ chaotic states, and enabling the study of unstable periodic orbits11,12. In the past, RQA has been used to analyze the relationship between atrial rate and spectral centre frequency13, to quantitatively analyse CFAEs14,15, to classify atrial electrograms as normal, fractionated or temporally unstable16–18, and to examine dynamics before and after ablative treatment19. This is possible even for short electrogram sequences17,20. Most similar to our study, RQA has been used, together with sample entropy, to distinguish between areas with or without rotors21, by applying RQA to bipolar electrogram readings that can identify wavefront direction. Sites behaving actively or passively in the context of fibrillation have been identified using 15 s22 or 5 s23 signals. In this work, we examine a variety of RQA and related metrics applied to much shorter (1s) AP traces. These traces form cardiac activation maps, which are an essential tool for the treatment of cardiac arrhythmias24. Furthermore, the method’s advantages for finding residual conduction after ablation of the mitral isthmus and thus ensuring its complete blockage are shown25 and thanks to the constant technological development of this method, new possibilities for its use are constantly being found26. We demonstrate that by appropriately combining RQA metrics, along with their spatial correlation, we can identify the changing nature of electrical signalling. In particular, we show how samples of single-site AP traces can be monitored for the presence of a stable rotor. We also demonstrate the use of RQA-based features to differentiate regular activation due to spiral waves (including separation of anchored and meandering spirals) and fibrillation (including separation of self-sustained and rotor-driven fibrillation) when high-resolution spatial data is available. These techniques also could serve as powerful tools for automated analysis of in silico simulation data, where full spatial AP data is available. Materials and methods Electrophysiological scenario simulation. The propagation of the action potential (AP) in cardiac tissue is modeled using the monodomain equation in terms of a dimensionless transmembrane potential, u, Membrane potential diffuses with isotropic conductivity D, representing the flow of ions between myocytes through their gap junctions. Membrane potential also changes with the flow of ions in and out of cells, Jion , and any externally provided stimulus current, Jstim . Here, Jion is represented by the model of Fenton and Karma27. This is a phenomenological model that simplifies ion transport to single currents representing the fast inward, slow inward, and slow outward flow of positive ions through the cell membrane ( Jfi , Jsi and Jso , respectively), Further information on the definition of this model, including parameter values used, can be found in the Supplementary material (TableS1). We consider two different scenarios relevant to cardiac arrhythmia. In the first scenario, we simulate a twodimensional slice of tissue of dimension 12.5 ×12.5 cm , this large size is chosen to allow single simulations to contain multiple dynamical behaviours in accordance with other similar studies23. No-flux boundary conditions are applied to domain edges. To produce the type of chaotic dynamics observed during fibrillation, we provide a rapid pacing stimulus, with a period of 31ms. Numerical discretisation parameters for this scenario are grid spacing x=0.25 mm and timestep t=0.1 ms . In the second scenario, two-dimensional tissue slices of dimension 4×4 cm are used, a size corresponding to the heart atrium. In this scenario, scar regions are added to make the tissue heterogeneous, creating the potential for spiral waves to be stabilised by the obstructed region28. The scar is created by placing two concentric circles, with a fixed probability of cell damage for sites within the inner circle and a linear decrease in probability to zero along the radius of the outer ring. This creates a border zone, known to be an important component of a scar regions’ arrhythmogenic impact29, as well as being a strong predictor of post-infarct mortality30. Damaged tissue is modelled as non-conductive (applying no-flux boundary conditions to boundaries between damaged and healthy tissue). In order to generate many different types of activation dynamics, a wide range of magnitudes of tissue damage are simulated, from no damage to complete destruction of tissue in the scar region. This scenario used x=0.2 mm and t=0.05 ms . In the simulations, the diffusion was Dxx = Dyy = 0.25 cm2 / s . If there was a scar in a given site, the diffusion coefficient was set to D=0 cm2/s . Signal analysis. Analysis of membrane potential time series is here approached using RQA11,31, an important method for nonlinear data analysis. RQA is based on the analysis of recurrence in a dynamic system (events of a dynamic system returning to a previously visited area in the phase space), as summarised by the recurrence plot (RP). The RP is a two-dimensional array of the values zero and one, with a value of one in position (i,j) corresponding to the system being in the same place in phase space at times ti and tj . That is, the non-zero elements in the array mark events of recurrence. The RP is calculated using the equation Here xi denotes the values of the system’s m dependent variables at the moment of its i-th snapshot, with N the total number of snapshots. The Heaviside function θ(x) , (1) ∂ u ∂ t =D  ∂ 2 u ∂ x 2+ ∂ 2 u ∂ y 2  −Jion −Jstim . Jion =−Jfi−Jso −Jsi. (2) Ri,j =θ(ǫ − ||x i −x j ||),x i ∈R m ,i,j=1...N . (3) θ( x):=  1, x>0 0, x≤0 , 3 Vol.:(0123456789) Scientific Reports | (2023) 13:11828 | https://doi.org/10.1038/s41598-023-38256-w www.nature.com/scientificreports/ specifies that recurrence occurs when the difference falls under the threshold ǫ for some choice of norm �·� . Here, we use the Euclidean norm. Once the RP is obtained, a number of different informative quantities can be calculated from the arrangement of its zero and non-zero elements. Many of these quantities consider the occurrence of diagonal strings of non-zeroes in the RP, as these correspond to system dynamics recurring throughout a window of time, instead of just coinciding at a single moment. In this paper, the following RQA measures will be used: RQA measure Description REC Percentage of recurrence points in a recurrence plot DET Percentage of recurrence points that form diagonal lines RATIO Ratio between DET and REC Lmax Length of the longest diagonal string of recurrence points Lmean Mean length of the diagonal lines DIV The inverse of Lmax ENTR The Shannon entropy of the diagonal line lengths distribution LAM Percentage of recurrence points that form vertical lines Vmax Length of the longest vertical line Vmean Mean length of vertical lines In addition to these RQA metrics, we also consider the organization indices (OI) for these time series. Widely used in cardiac electrophysiology studies32–39, organization indices are a measure of order (or disorder) in a time series, calculated using its representation in the frequency domain. Specifically, the nth organization index is defined as the proportion of power in the power spectrum that is contained within the n highest-power peaks. A high organization index, with large proportions of the power in the power spectrum concentrated in a small number of peaks, corresponds to regular activation of a location at some fixed frequency(ies). Calculations of the power contained within a peak will sometimes also include the power associated with its harmonic frequencies, but here these are not included. Finally, we analyze time series in terms of their entropy. Accurate entropy calculation requires large amounts of data and is very sensitive to the noise of the system under investigation, and so is not suitable here. To overcome this limitation, Pincus et al.40 developed an approximate entropy measure (ApEnt), based on searching for similar subsequences in the analyzed time series. In order to avoid the occurrence of natural logarithms of zero in the ApEnt calculation, each sequence is counted as similar to itself, although this introduces a bias41. Richman and Moorman41 proposed sample entropy (SampEnt) as a modified entropy measure that is independent of the length of the data, is bias-free, and also requires fewer operations to calculate42. Further information about the calculation and use of these entropy metrics is offered by Delgado-Bonal and Marshak43. We consider both ApEnt and SamptEnt in this work. All metrics are calculated from time series of length 1000ms, with a sampling frequency of 1000Hz. Spatial analysis. As a means of automatically identifying different types of arrhythmic dynamics, we take together the values of the time series metrics discussed in the previous section at different spatial locations. The approach is based on the intuition that overall measures of order/disorder in the tissue distinguish between tachycardic and fibrillatory patterns of activation, and the level of spatial similarity in these metrics then provides information towards the phenomena driving those dynamics. Initially, to demonstrate the potential in spatially analysing RQA (or other) time series metrics, we choose a fine spatial resolution with which to carry out this analysis. By selecting one in four node points from the discretisation used for numerical simulation, this results in a spacing of 0.4mm between the locations for which a time series is considered available. We later discuss the effect of coarser spatial resolutions. To quantify spatial similarity, we use the spatial autocorrelation measure known as Moran’s I44. Given a weight matrix, W, that defines how different spatial locations relate to one another, this measure is defined Here N is the total number of spatial sites, and x is the quantity for which the spatial correlation is being calculated. ¯ x = 1 Ni x i is the mean value of this quantity, and wij are the elements of the matrix W each describing the level of relationship between sites i and j. So long as the weight matrix is chosen such that sites close together are weighted strongly and distant sites are weighted weakly, a Moran’s I value close to the maximal value of +1 corresponds to large regions of similar values of the quantity of interest, while a value close to the minimum value of −1 corresponds to sharp spatial variation in the quantity (negative correlation). A value of zero represents independence between sites. In this work, the weighting matrix is defined by the Moore neighbourhood (surrounding eight points including diagonals), with these neighbours being given full weight and other points given zero weight, I =N ij wij  i  jwij(xi−¯x)  xj−¯x  i (xi−¯x)2 . w ij =  1j∈Moore neighbourhood of i 0 otherwise . . 4 Vol:.(1234567890) Scientific Reports | (2023) 13:11828 | https://doi.org/10.1038/s41598-023-38256-w www.nature.com/scientificreports/ As each site is not its own neighbour, the diagonal elements of the weight matrix wii =0 . As an example, given a set of values x for sites arranged in a 3 × 3 grid, we would have a 9 × 9 matrix given by Results RQA measures encode key facets of localised activation patterns. We first analyse how the range of RQA measures we consider respond to different patterns of activation. For this purpose, we simulate the complex response of cardiac tissue to very fast pacing (further details in Methods), generating regions that experience regular activation, regions that experience irregular activation, and yet more intricate dynamics where these two phenomena interact (Fig.1a). Regular waves propagate outwards from the centre of the tissue (the location of the pacing stimulus), and continue undisturbed out towards the upper left corner. In the lower right corner, a stable rotor develops and persists throughout the simulation period. Two more transient rotors appear during the course of the simulation, x 1x2x3 x 4x5x6 x 7x8x9 ,W=             010110000 101111000 010111000 111000110 111001111 011010111 000111010 000111101 000011010             . Figure1. (a) Membrane potential snapshots, spaced by 143ms. Regular fast pacing in the centre of a large ( 12×12 cm ) slice of tissue creates regular waves radiating out toward the upper left corner. A spiral persists in the bottom right corner, and transient rotors appear in the lower left and upper right corners. (b) Tissue maps for the RQA measures from left to right, Lmax , RATIO, DET, and LAM. Lmax takes on high values where activation is strictly regular, in the upper left corner and where rotor cores persist. RATIO highlights similar regions (via lower values), but identifies both rotor tips and the regular activation dynamics surrounding them. Measures DET and LAM highlight (via lower values) the regions where activation dynamics are most chaotic, characterised by wavelets that form where the main propagating waves interact. DET emerges as most appropriate for this purpose, as LAM also reaches low values around the rotor tip in the lower right corner. (c) Blurring the original maps via 2D convolution (see main text) makes the regions identified by the different RQA measures more visually distinguishable. In particular, the boundaries between regions of regular and fibrillatory activation detected by measure DET become much more pronounced. 5 Vol.:(0123456789) Scientific Reports | (2023) 13:11828 | https://doi.org/10.1038/s41598-023-38256-w www.nature.com/scientificreports/ in the lower left and upper right corners. The remainder of the tissue experiences the interactions of the different waves or wavelets created by these different rotors and the pacing site, exhibiting irregular patterns of activation and AP propagation characteristic of fibrillation. We analyze the RQA measures by calculating their values for the time series of the (dimensionless) membrane potential, u, at many sites throughout the tissue, producing spatial maps of these measures. These maps, some of which are displayed in Fig.1b, are then compared with the spatially varying activation dynamics through the tissue. We also consider smoothed versions of these spatial maps, obtained by averaging the value of each pixel and its immediate neighbours (the 2D convolution principle). This reduces noise and better highlights differences in value between larger-scale regions of the tissue (Fig.1c). The RQA measure Lmax best detects the regular waves propagating from the centre to the upper left corner, and the rotor tips at the bottom of the simulation domain. The measure RATIO acts somewhat oppositely to Lmax , taking on higher values where the activation pattern is disordered. In contrast with Lmax , however, regions where activation is quite regular but does not occur at a single, fixed frequency are still labelled as ordered (such as the lower right corner of the domain). As such, RATIO is less sensitive to regularity violations than Lmax . The measures DET and LAM are the least sensitive to regularity violations. DET takes on high values (greater than 0.9) throughout the entirety of the tissue, but after spatial smoothing, it becomes visually clear that this measure takes lower values where the dynamics are the least predictable. That is, the portions of the tissue where the activation waves from the different rotors collide and interact with each other. LAM performs a similar role, largely matching DET in terms of where it takes on a relatively low or high value. These different sensitivities help to guide appropriate selection of RQA measure(s), depending on the investigation being conducted. For example, DET and LAM are best suited for identifying downstream regions of wavebreak, unlikely to house any persistent drivers of arrhythmia. On the other hand, in the following section, we demonstrate the use of the most sensitive measure Lmax , along with a second RQA measure, to very accurately identify locations of stable functional re-entry. Detection of rotor tips. As Fig.1 demonstrates, persistent functional re-entry can manifest in more chaotic, fibrillatory dynamics away from the rotor location. Such persistent re-entries are often termed “mother rotors”, and are one of the prominent theories of arrhythmias not driven by some external circumstance (for example structural heterogeneity due to scarring)45. Others argue that all rotors are liable to appear and disappear, and that the dynamics behind these arrhythmias are better interpreted as a stochastic process of creation and annihilation46. Regardless, should such stable, self-sustained spiral waves exist, their accurate detection is critical for informing interventions such as ablative surgery. Automatic detection of rotors using spatial information remains a topic of interest, although has concentrated more on identification of any spiral wave tips (including both sustained rotors and transient tips appearing in fibrillation)47. Here we consider the use of RQA metrics specifically for the detection of sustained, spatially-fixed rotors using only local time series information. As discussed previously, Lmax serves as a very sensitive measure of regularity in a site’s membrane potential time series, and is thus an important component in the detection of spiral waves. However, to assess the fixation of these spiral waves, we also involve the RQA measure Lmean . By indicating the average length of diagonal lines in the RP, Lmean is related to the mean prediction time48 and helps to distinguish between persistent spiral waves and other regular dynamics, such as persistent stimulation due to waves of activation emanating from a site of regular pacing. To enable automatic detection using these measures, we select restrictive thresholds for their values that pick out only the most evident examples of the types of dynamics they characterise. The regions identified using threshold values for these two RQA measures are shown inFig.2, together with the intersection of these two sets of regions that serves as a detector of stable rotor cores. Here, this identifies the rotor in the bottom right, the only one that persists throughout the pictured time window in the simulation. The transient rotors that form, but are then annihilated, fail to reach the critical Lmax value. Meanwhile, regions of highly regular dynamics not reflective of a rotor are removed by failing to meet the threshold value of Lmean . Notably, as the value of Lmax represents the longest period of time the system remains in a periodic pattern, the threshold value selected for this RQA measure also reflects how long a rotor must persist before it will be labelled as stable by the detection process. This provides a degree of control over the detection process to the user. As we demonstrate in FigureS1, reducing the threshold value of Lmax presents a means of also identifying regions that are occupied by a stable rotor for a reduced amount of time. Considering the spatial maps seen in Fig.2, a further potential role for RQA in the detection of sustained rotors emerges when there is even limited access to spatially-distributed time series data. As discussed in the previous section, RATIO serves as a measure of irregularity, that is less sensitive to minor disruptions of a periodic activation pattern. This allows it to identify (with low values, indicating regularity), the region immediately surrounding the stable rotor (second column in Fig.2b,c). Distant from the rotor (top left), the RATIO value is similarly low, as this region is also activated in a regular fashion. The key feature distinguishing a stable rotor, at least in the presence of fibrillation, is a region of low RATIO surrounded by considerably higher values, and this could be detected by RQA applied to only a few time series collected from locations with the appropriate spatial separation. Use of RPs of AP traces to classify arrhythmia dynamics. We next consider a series of simulations using smaller slices of tissue (comparable to the size of the human heart atrium), each including a region of scarring of varying size and severity. In contrast to the study described in previous sections, characterised by a range of different dynamics and their interactions, these simulations using a smaller domain produce activation dynamics corresponding to tachycardia or fibrillation that can be classified on the tissue scale (Fig.3). 6 Vol:.(1234567890) Scientific Reports | (2023) 13:11828 | https://doi.org/10.1038/s41598-023-38256-w www.nature.com/scientificreports/ For each simulation, RPs are generated from the membrane potential time series at different locations using the approach described in the Methods. The time lag and embedded dimension are found using the nonlinearTseries package for the R programming language49. These parameters are used for phase space reconstruction. According to Takens theorem50, it is possible to use the time series that comes from a dynamic system to reconstruct a trajectory using time lag and embedding. Reconstructed trajectories will have the same dynamic properties as the original dynamic system. An RP is then constructed from these trajectories. For example, the time series X1,X2,X3,...,X10 with time delay 3 and embedding dimension 2 will be reconstructed into a twodimensional phase space as (X1,X4),(X2,X5),(X3,X6),...,(X7,X10) . The nonlinearTseries package computes the time lag using the first zero crossing of the autocorrelation function, and the embedding dimension using L. Cao’s algorithm51. The threshold ǫ for recurrence detection is determined to be 3% of the range of the analyzed time series. The resultant RPs are pictured in Fig.4, where they are seen to capture successfully some of the features of the different types of activation pattern. Tachycardia driven by a stable rotor. Where the re-entry anchors to the region of fibrosis, and the wave does not exhibit breakup away from the rotor core, this produces a regular but tachycardic pattern of activation, see Fig.4a. AP morphology is consistent throughout the tissue (except the damaged site), similar to sinus rhythm, except for the pacing frequency. The RPs in this case consist of diagonal lines that reveal the periodic nature of Figure2. Snapshots (spaced by 143ms) of AP propagation, with regions for which the membrane potential timeseries produces an RQA measure ( Lmax , Lmean , or both at the same time) exceeding a threshold value highlighted in bright green. These measures correspond to the duration of patterns of periodic activation (indicated in an RP by the length of diagonal lines), and as such the value of the Lmax threshold selects how long a regular activation pattern must persist to be identified. In this case, the threshold value Lmax ≥800 detects only the tip of the spatially pinned rotor, when combined also with a threshold Lmean > 15 (fifth and sixth rows, highlighted region in the lower right corner within the light blue circle). 7 Vol.:(0123456789) Scientific Reports | (2023) 13:11828 | https://doi.org/10.1038/s41598-023-38256-w www.nature.com/scientificreports/ activation. Patterns are superimposed on these lines at regular intervals, created by the portions of the time series where the tissue is at rest. Tachycardia driven by a wandering rotor. Re-entries that do not fix to a heterogeneity in the tissue and do not produce wavebreak also result in tachycardic activation of the tissue. The path of the rotor tip follows a twisting pattern, resulting in inconsistent lengths of time between activation events and as a result, differences in morphology between individual APs, see Fig.4b. RPs associated with these dynamics do not show the many regularly spaced diagonal lines indicating periodicity. However, as each completion of the rotor’s twisting path does represent a single period, diagonal lines may still appear far away from the main diagonal, with consistent spacing. Fibrillation driven by wavelets. Fibrillation may also be sustained by many individual wavelets, created and annihilated by their unpredictable interactions with each other, see Fig.4c. In these circumstances, membrane potential time series show highly varied AP morphology and activation timing, resulting in RPs consisting only of lines along the main diagonal, possibly with distinct blocks corresponding to prolonged periods of rest (example in Fig.4). Fibrillation driven by a stable rotor. Wave break can occur distant from an anchored spiral wave re-entry, resulting in disordered activation of the surrounded tissue, see Fig.4d. This results in very distinct membrane potential time series depending on the region of the tissue from which readings are taken. RPs distant from the anchored re-entry are chaotic, consisting mostly of a line along the main diagonal. In contrast, RPs at the location of the spiral wave show many evenly-spaced diagonal lines, characteristic of consistent, regular activation. Owing to minor disruptions by the surrounding chaotic dynamics, however, these lines are interrupted in a fashion not seen for tachycardia driven by a stable rotor. Using the membrane potential time series from just two locations, each different class of arrhythmic dynamics can be identified using the RPs leveraging the summarised distinguishing features. However, this relies upon qualitative observation and an appropriate selection of the locations to ensure capture of both regularly activated and irregularly activated tissue, in the case both are present. Such important locations will not generally be known prior to data collection. Consequently, we now examine whether these different arrhythmia dynamics Figure3. Snapshots (spaced by 100ms) of each group of AP propagation modelled on a smaller tissue square (with the size of 4×4 cm ). At each tissue site, the time series is analyzed using RQA, OI, RP eigenvalues and entropy. Examples of individual time series at the points in red and magenta are given in Fig.4. The dynamic parameters of these time series are then used to classify individual propagations. (A) A spiral wave anchored to a highly fibrotic region resulting in sustained re-entry. (B) Meandering spiral wave. (C) Tissue experiencing fibrillatory activation. (D) Spiral wave anchored to a scar region (bottom right), that breaks up into wavelets away from the main rotor. 8 Vol:.(1234567890) Scientific Reports | (2023) 13:11828 | https://doi.org/10.1038/s41598-023-38256-w www.nature.com/scientificreports/ on the tissue scale can be identified automatically, using the membrane potential time series recorded at many points throughout the tissue. Spatial behaviour of metrics for different arrhythmia mechanisms. To achieve our goal of distinguishing these four important types of arrhythmia in an automated fashion, we first consider how a large range of RQA metrics, along with OI and entropy as other measures of time series regularity, vary spatially throughout the tissue. We rely on the sophistication of these metrics to capture the intricacies of the different activation patterns and use relatively simple means to characterise this spatial variation. Specifically, we consider the mean value of each metric across all measurement locations in the tissue and the spatial correlation of each metric as evaluated by Moran’s I. This provides an interpretable and lower-dimensional model for the automatic classification we will suggest. Figure4. Membrane potential traces from two distinct locations (depicted in Fig.3) in the tissue (left side) and its recurrence plots (right side). (A) A spiral wave anchored to a highly fibrotic region resulting in sustained re-entry. The spiral wave regularly stimulates the tissue at both sites. This phenomenon is manifested in RPs by long uninterupted diagonal lines. On these diagonals there are small squares created by the plateau phase. (B) The regular meandering of the spiral wave manifests by repeating patterns in the time series. Due to the long period, uninterrupted lines are not visible in the RP as in the case of an anchored spiral wave. (C) Tissue experiencing fibrillatory activation. The time series created by this activation does not contain any regular patterns. Therefore long uninterrupted lines are not visible in the RP. The large black square visible in the right RP is due to the membrane potential remaining in the resting state between 200 and 400 ms (in the right time series). (D) Spiral wave anchored to a scar region (bottom right), that breaks up into wavelets away from the main rotor. This propagation manifests irregular patterns (left time series) and regular AP waveforms (right time series). However, the surrounding fibrillatory activation disrupts the regularity of these waveforms. As a result, the diagonal lines are interrupted. 9 Vol.:(0123456789) Scientific Reports | (2023) 13:11828 | https://doi.org/10.1038/s41598-023-38256-w www.nature.com/scientificreports/ In Fig.5, spatial maps of some of the different RQA measures are displayed, demonstrating how each can be used to identify or distinguish the features of the different types of arrhythmia. Across all the metrics, tachycardia driven by localised re-entry results in values of the metrics that are spatially homogeneous, suggesting that very high values of spatial correlation may identify this pattern of behaviour. In the case of a rotor in a wandering orbit around a location, this region is highlighted by both RATIO and LAM. Some metrics also produce multiple teardrop shapes around this region, indicating the points where the rotor tip turns around to continue its orbit. In fibrillatory arrhythmias, the spatiotemporal chaos results in far more irregular spatial maps of the different RQA metrics. When fibrillation is driven by an anchored rotor, larger contiguous regions of similar values appear, suggesting the possibility that spatial correlation might allow the two different causes of fibrillation considered here to be distinguished. To explore this, we now consider how our summarisation measures for the spatial RQA maps vary within a larger set of simulations, thus demonstrating the potential for automated classification. We give a detailed characterisation of a number of our measures in Table1. Figure5. Example of calculated tissue maps (with the size of 4×4 cm ) depicting the spatial distribution of selected RQA measures. Rows: A) anchored rotation, B) traveling rotor, C) fibrillation, D) fibrillation with anchored rotation. Columns: RQA measures 1) RATIO, 2) LAM, 3) ENTR, and 4) DET. A pattern indicating the type of analyzed propagation can be found in all maps. In anchored rotation (A), a spiral created by the rotational movement of the AP can be seen. Three spirals formed by the meandering rotor can be found (B) (in each spiral, the AP re-entries the tissue). The region around which the AP rotates is highlighted in the middle of these spirals. Fibrilatory activation (C) is manifested in maps by irregular, chaotically distributed shapes. Fibrillation associated with anchored rotation (D) in the maps (especially in measure RATIO) shows larger contiguous regions created by regular propagation made by anchored rotation. If we compare the mean values of individual RQA measures, we can see a lower mean value of RATIO (Ratio between DET and REC) for anchored rotation (see boxplots in Fig.6). This phenomenon is caused by the higher value of the RR measure due to diagonal lines manifesting periodic motion in RP (see Fig.4). There is also a high mean value for measure LAM (percentage of vertical lines in RP). A regular resting phase causes this phenomenon during anchored rotation (manifesting as black squares in RP). This phase can also be partially found in the meandering rotor, as its average value in the tissue is the second highest (see Fig.6). The average value of the ENTR measure is the highest for anchored rotation. ENTR achieves high values, as RPs of this type of propagation are the most complex (due to diagonal lines and their connected patterns). The Mean value of the DET measure is highest for anchored rotation due to the repetitive nature of this propagation. The irregularity of fibrillatory activation is reflected in the lowest value of this measure for all the types of analyzed propagation. 16 Vol:.(1234567890) Scientific Reports | (2023) 13:11828 | https://doi.org/10.1038/s41598-023-38256-w www.nature.com/scientificreports/ 63. Antzelevitch, C. & Fish, J. Electrical heterogeneity within the ventricular wall. Basic Res. Cardiol. 96, 517–527 (2001). 64. Britton, O. J. et al. Experimentally calibrated population of models predicts and explains intersubject variability in cardiac cellular electrophysiology. Proc. Natl. Acad. Sci. 110, E2098–E2105 (2013). Acknowledgements This work was supported by the project “The Doctoral School for Education in Mathematical Methods and Tools in HPC” reg. nr.: CZ.02.2.69/0.0/0.0/16_018/0002713, by the Ministry of Education, Youth and Sports of the Czech Republic through the e-INFRA CZ (ID:90254). This work was partially supported by CNPq, FAPEMIG, CAPES, and by the ARC Centre of Excellence for Mathematical and Statistical Frontiers (CE140100049). Author contributions R.H. and B.L. created the simulation tools used. R.H. created code for signal analysis. All authors contributed to the analysis of results and subsequent development of the study, original study concept, and drafting of the manuscript. competing interests The authors declare no competing interests. Additional information Supplementary Information The online version contains supplementary material available at https:// doi. org/ 10. 1038/ s4159802338256-w. Correspondence and requests for materials should be addressed to R.H. Reprints and permissions information is available at www.nature.com/reprints. Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 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 http:// creat iveco mmons. org/ licen ses/ by/4. 0/. © The Author(s) 2023